STM32F411 BlackPill CMSIS-DSPによる楽器用チューナー（７）コード

計測用コードのみ記載
１．コード
１）計測
（１）<fft_peak.h>

/**
 Project Name:FFT_F4

 File Name:fft_peak.h

 Description:
 USB speaker with isochronous transfer
 Generation Information :
 Device			:  STM32F411CEU6
 Version			:  01
 Created on		: 2023/12/20
 */
/**
 Copyright (c) 2023 luke24e-hbid

 This software is released under the MIT License.
 http://opensource.org/licenses/mit-license.php
 */

/* Define to prevent recursive inclusion -------------------------------------*/
#ifndef INC_FFT_PEAK_H_
#define INC_FFT_PEAK_H_

/* Includes ------------------------------------------------------------------*/
#include "main.h"
#define ARM_MATH_CM4
#include "arm_math.h"
#include "arm_const_structs.h"
#include "arm_common_tables.h"

/* Private defines -----------------------------------------------------------*/

/* Function prototype declaration ------------------------------------------*/
//	Read FFT Data
void read_data_sin(float32_t fin);

// FFT peak detect
float32_t rfft_spline(uint32_t nfft);

#endif /* INC_FFT_PEAK_H_ */

（２）<fft_peak.c>

/**
 Project Name:FFT_F4

 File Name:fft_peak.c

 Description:
 USB speaker with isochronous transfer
 Generation Information :
 Device			:  STM32F411CEU6
 Version			: 01
 Created on		: 2023/12/20
 */
/**
 Copyright (c) 2023 luke24e-hbid

 This software is released under the MIT License.
 http://opensource.org/licenses/mit-license.php
 */
/* Includes ------------------------------------------------------------------*/
// #include "main.h"
#include "fft_peak.h"
#include <stdlib.h>

/* Private macro -------------------------------------------------------------*/
// RFFT macro
#define FS1 16000                                 // sampling frequency 1 : 16000[Hz]
#define FT 4                                      // sampling period [Hz]
#define NS1 (FS1 / FT)                            // number of read sample data : 4000
#define NFFT 4096                                 // number of RFFT data
#define NPOWER (NFFT / 2)                         // number of RFFT power data
#define FRESO1 ((float32_t)FS1 / (float32_t)NFFT) // frequency resolution 1 : 16000/4096[Hz]

// measurement macro
#define FKMIN 55.0   // minimum frequency
#define FKMAX 880.0  // maximun frequency
#define PKMIN 200000 // minimum power
#define PKK 2        // minimum power rate

// Spline macro
#define NSPS 12         // number of spline data
#define NSPH 100        // number of spline interpolation
#define PSPH (NSPH / 2) //  magnification of spline interpolation

/* Private variables ---------------------------------------------------------*/
// global variable
int16_t data_sample[NFFT * 2]; // sampling data

/* Private user code ---------------------------------------------------------*/
//
// FFT peak detect
float32_t rfft_spline(uint32_t nfft)
{
  // RFFT parameter
  arm_rfft_fast_instance_f32 fft;
  float32_t fft_input[NFFT];
  float32_t fft_output[NFFT];
  float32_t fft_power[NFFT / 2];
  uint8_t ifftFlag = 0;
  float32_t maxValue;
  uint32_t maxIndex;

  // arm_spline_init_f32  parameter
  arm_spline_instance_f32 S;
  arm_spline_type type = ARM_SPLINE_NATURAL;
  float32_t x[NSPS];
  float32_t y[NSPS];
  uint32_t n = NSPS;
  float32_t coeffs[3 * (NSPS - 1)];
  float32_t tempBuffer[2 * (NSPS - 1)];

  // arm_spline_f32   parameter
  float32_t xq[NSPH];
  float32_t pDst[NSPH];
  uint32_t blockSize = NSPH;
  float32_t maxValue_s;
  uint32_t maxIndex_s;

  // measurement
  float32_t avg_power; // average of FFT power
  float32_t fm1;       // dominant frequency 1
  uint32_t fm1Index;
  uint32_t fm1Index_01;
  uint32_t fm1Index_02;
  uint32_t fm1_multi;
  float32_t fm1_01;
  float32_t fm1_02;
  float32_t xn;
  float32_t fs2; // sampling frequency 2
  float32_t freso2;
  uint16_t nk2;
  float32_t fm2;   // dominant frequency 2
  float32_t fmspl; // spline interpolation frequency

  uint16_t i;

  // Set FFT data 1
  for (i = 0; i < NS1; i++)
  {
    fft_input[i] = (float32_t)data_sample[i + nfft];
  }
  for (i = NS1; i < NFFT; i++)
  {
    fft_input[i] = 0;
  }

  // analyze signal 1
  arm_rfft_fast_init_f32(&fft, NFFT);
  arm_rfft_fast_f32(&fft, fft_input, fft_output, ifftFlag);
  arm_cmplx_mag_f32(fft_output, fft_power, NFFT / 2);
  // fill fft power with zeros
  fft_power[0] = 0.0;
  arm_mean_f32(fft_power, NPOWER, &avg_power);

  for (i = 0; i < NPOWER; i++)
  {
    if (i < (uint16_t)(FKMIN / FRESO1))
      fft_power[i] = 0;
    else if (i > (uint16_t)(FKMAX / FRESO1))
      fft_power[i] = 0;
    else
    {
      if (fft_power[i] < PKMIN)
        fft_power[i] = 0;
      else if (fft_power[i] < (avg_power * PKK))
        fft_power[i] = 0;
    }
  }
  // find dominant frequency 1-1
  arm_max_f32(fft_power, NPOWER, &maxValue, &maxIndex);
  fm1Index_01 = maxIndex;

  if (fm1Index_01 == 0)
  {
    return (0.0);
  }

  fm1_01 = (float32_t)fm1Index_01 * FRESO1;
  // fill fft power with zeros
  i = (uint16_t)fm1Index_01;
  do
  {
    fft_power[i++] = 0;
  } while (fft_power[i] != 0);
  i = (uint16_t)fm1Index_01;
  do
  {
    fft_power[i--] = 0;
  } while (fft_power[i] != 0);
  for (i = (uint16_t)(fm1Index_01 * 2 + fm1Index_01 / 2); i < NPOWER; i++)
  {
    fft_power[i] = 0;
  }

  // find dominant frequency 1-2
  arm_max_f32(fft_power, NPOWER, &maxValue, &maxIndex);
  fm1Index_02 = maxIndex;
  fm1_02 = (float32_t)fm1Index_02 * FRESO1;

  // find dominate frequency 1
  if (fm1Index_01 > fm1Index_02)
  {
    uint32_t t = fm1Index_01;
    fm1Index_01 = fm1Index_02;
    fm1Index_02 = t;
  }
  fm1_multi = 1;
  if (fm1Index_01 == 0)
    fm1Index = fm1Index_02;
  else
  {
    if (abs((int32_t)fm1Index_01 * 2 - (int32_t)fm1Index_02) < 2)
      fm1Index = fm1Index_01;
    else
    {
      fm1Index = fm1Index_01;
      fm1_multi = (uint32_t)((float32_t)fm1Index_01 / (float32_t)(fm1Index_02 - fm1Index_01) + 0.5);
    }
  }
  fm1 = (float32_t)fm1Index * FRESO1;

  // calculate sampling frequency 2
  xn = (log(fm1 * 16.0) - log(125.0)) / log(2.0) + 0.5;
  fs2 = 125 * (1 << (uint16_t)xn);
  freso2 = fs2 / NFFT;

  // Set FFT data 2
  nk2 = FS1 / (uint16_t)fs2;
  for (i = 0; i < (uint16_t)(fs2 / 4); i++)
  {
    fft_input[i] = 0;
    for (uint16_t j = 0; j < nk2; j++)
    {
      fft_input[i] = fft_input[i] + (float32_t)data_sample[i * nk2 + j + nfft];
    }
  }
  for (i = (uint16_t)(fs2 / 4); i < NFFT; i++)
  {
    fft_input[i] = 0;
  }
  // analyze signal 2
  arm_rfft_fast_init_f32(&fft, NFFT);
  arm_rfft_fast_f32(&fft, fft_input, fft_output, ifftFlag);
  arm_cmplx_mag_f32(fft_output, fft_power, NFFT / 2);
  // find dominant frequency 2
  maxValue = 0;
  maxIndex = 0;
  for (i = (fm1Index - 2) * nk2; i < (fm1Index + 2) * nk2; i++)
  {
    if (fft_power[i] > maxValue)
    {
      maxValue = fft_power[i];
      maxIndex = i;
    }
  }
  fm2 = (float32_t)maxIndex * freso2;

  // spline
  for (i = 0; i < NSPS; i++)
  {
    x[i] = (float32_t)(i + maxIndex - (NSPS / 2));
    y[i] = fft_power[i + maxIndex - (NSPS / 2)];
  }
  for (i = 0; i < NSPH; i++)
  {
    xq[i] = (float32_t)i / (float32_t)PSPH + (float32_t)(maxIndex - 1);
  }
  arm_spline_init_f32(&S, type, x, y, n, coeffs, tempBuffer);
  arm_spline_f32(&S, xq, pDst, blockSize);
  // find spline dominant frequency
  pDst[0] = 0.0;
  arm_max_f32(pDst, 100, &maxValue_s, &maxIndex_s);
  fmspl = xq[maxIndex_s] * freso2 / (float32_t)fm1_multi;

  // printf("%.3f\t%.0f\t%.3f\t%.3f\r\n", (double)fm1, fs2,(double)fm2, (double)fmspl);

  return (fmspl);
}

２）呼び出し例
　nfft : サンプリングデータ先頭
　fmspl : 計測値[Hz]

/* Private variables ---------------------------------------------------------*/
uint32_t g_nfft_adc = 0xffff;

/* Private user code ---------------------------------------------------------*/
// ADC DMA harf
void HAL_ADC_ConvHalfCpltCallback(ADC_HandleTypeDef *const hadc)
{
  g_nfft_adc = 0;
}
// ADC DMA full
void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *const hadc)
{
  g_nfft_adc = 4000;
}

//	RFFT, Display
//
void display_tune(uint32_t nfft)
{
  float32_t fmspl;  // display frequency
  int16_t oct_cent; // octave cent
  int16_t cent1200; // cent
  int16_t note;     // note number
  int16_t cent;     // display cent +/- 50

  fmspl = rfft_spline(nfft);
  oct_cent = (uint16_t)(1200.0 * log(fmspl / 27.5) / log(2.0) + 50.5);
  cent1200 = oct_cent % 1200;
  note = (cent1200) / 100;
  cent = (cent1200) % 100 - 50;
}

//
//  main

extern int16_t data_sample[]; // sampling data
HAL_ADC_Start_DMA(&hadc1, (uint32_t *)data_sample, 8000);
HAL_TIM_Base_Start(&htim3);

while (1)
{
  if (g_nfft_adc != 0xffff)
  {
    display_tune(g_nfft_adc);
    g_nfft_adc = 0xffff;
  }
}

２．補足解説
１）コードについて
・コードは FFT計算とスプライン補間による基音周波数計測と基音周波数のセント値変換関数を記載
・呼び出し例としてAD変換DMA転送によるデータを計測を記載
・DMA転送はダブルバッファで構成
・表示用コードは音名表示に使用するフォントデータ量が大きいため記載していない

２）表示イメージ
（１）リソースの想定

リソース	デバイス	用途
I2C Master	OLED	音名、周波数等表示
LED 2port	2色LED 赤・黄緑	チューニング状態表示 　緑色：±5セント内 　橙色：±5セント外 　赤色：ADC入力レベル小 　　　　周波数測定範囲外
PWM 2ch	アナログメータ（電流計 両振れ）	セント値表示

（２）表示イメージ：YAMAHA TD-90
　音名、周波数等はOLEDで表示

jp.yamaha.com

おしまい