ARMのCMSIS-DSPライブラリで実数FFT

使い方

#include <arm_math.h>

arm_rfft_fast_instance_f32 S;

arm_status status = arm_rfft_fast_init_f32(&S, nfft);

arm_rfft_fast_f32(&S, data, data, 0);

検証

前回と同じ計算をARMのCMSIS-DSPライブラリを使って実行してみる。実験にはSPRESENSEを用いた。

#include <Arduino.h>
#include <math.h>
#include <arm_math.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// FFT構造体
arm_rfft_fast_instance_f32 S;

// fs: サンプルレート [Hz]
// f: 波形の周波数 [Hz]
// i: サンプル番号
// 戻り値: サンプル値 (-1.0 〜 +1.0)
double my_wave(double fs, double f, int i)
{
    double t = (double)i / fs; // 時刻 [秒]
    double x = cos(2.0 * M_PI *     f * t) 
       + 0.5 * sin(2.0 * M_PI * 2 * f * t)
       + 0.25 * cos(2.0 * M_PI * 3 * f * t);
    return x;
}

void setup() {

  Serial.begin(115200);
  Serial.println("Hello!");

  const double fs = 16000.0; // サンプルレート[Hz]
  const double f = 1000.0;   // 波形の周波数[Hz]

  const int nfft = 512;
  
  arm_status status = arm_rfft_fast_init_f32(&S, nfft);
  if (status != ARM_MATH_SUCCESS) {
    Serial.printf("DSP init error\n");
    return;
  }

  float data[nfft];  // 入力は実数の配列

  for (int i = 0; i < 512; i++) {
    double x = my_wave(fs, f, i);
    data[i] = x;
  }

  arm_rfft_fast_f32(&S, data, data, 0);
    
  //dsps_fft4r_deinit_fc32();

  for (int i = 0; i < 256; i++) {
    Serial.printf("[%d] %.2f %.2f\n", i, data[2 * i], data[2 * i + 1]);
  }
}

void loop() {
  delay(1000);
}

これを実行すると、前回と同等の結果が得られた。

[32] 256.00 -0.00
[64] 0.00 -128.00
[96] 64.00 0.00