実数列を半分サイズの複素FFTにつっこむイディオム

謎のコード

ESP32の音声処理のプログラムに次のようなコードがあった。

  dsps_fft4r_fc32(data, nfft / 2);
  dsps_bit_rev4r_fc32(data, nfft / 2);
  dsps_cplx2real_fc32(data, nfft / 2);

dataは音声のサンプル列つまり実数の配列であり、nfft はFFTサイズである。これは一見して奇妙である。dsps_fft4r_fc32 は複素FFTの関数であって、第一引数は実部と虚部が交互に並んだ配列を、第二引数にはFFTサイズを要求される。それなのに実数の配列と半分のFFTサイズを渡している。これはどうしたことか？
また、dsps_bit_rev4r_fc32 はバタフライ演算結果の並べ直しと思われるが、dsps_cplx2real_fc32 はいったい何なのか？ドキュメントを読んでもよく分からない。

Gemini先生の解説

Gemini先生によると、「N点の実数FFTを、N/2点の複素FFTとして計算することで、計算量とメモリを半分に節約する」というテクニックらしい。

N点の実数データを、「偶数番目を実部、奇数番目を虚部」と見なして、N/2点の複素データとして FFT につっこむ。(dsps_fft4r_fc32関数)
バタフライ演算結果を並べ直す。(dsps_bit_rev4r_fc32関数)
N/2点の複素FFT結果を、N点の実数FFT結果へと換算する。(dsps_cplx2real_fc32関数)

ということらしく、dsps_cplx2real_fc32関数がこのトリックの肝なのだが、計算式を見てもよく分からず、ホンマかいな？という気持ちが拭えない。

検証

PCで実数FFT

まず、PC上で実数FFTを計算してみる。FFTの計算にはKiss FFTを利用した。

サンプリング周波数 : 16kHz
入力波形 : 1kHzと2kHzと3kHzの正弦波の合成波形
FFT点数 : 512点

#include <stdio.h>
#include <math.h>
#include "kiss_fftr.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// 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;
}

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

    FILE* fp = fopen("fftout.txt", "w");

    const int nfft = 512;
    kiss_fftr_cfg cfg = kiss_fftr_alloc(nfft, 0, NULL, NULL);

    float in[nfft];                 // 入力は実数の配列
    kiss_fft_cpx out[nfft / 2 + 1]; // 出力は (nfft/2 + 1) 個の複素数

    for (int i = 0; i < 512; i++) {
        double x = my_wave(fs, f, i);
        in[i] = x;
    }
    kiss_fftr(cfg, in, out);
    
    kiss_fftr_free(cfg);

    for (int i = 0; i <= 256; i++) {
        fprintf(fp, "[%d] %.2f %.2f\n", i, out[i].r, out[i].i);
    }

    fclose(fp);

    return 0;
}

これを実行すると、[32]と[64]と[96]にピークが出る。Δf = 16kHz / 512 = 31.25Hz であるから、32×31.25Hz = 1kHzで計算が合う。cosが実部に、sinが虚部に出ていることにも注意。

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

ESP32の問題のコード

今度は同じ計算をESP32の問題のコードで実行してみる。実験にはM5Stack Core2を用いた。

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// 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;

  if (dsps_fft4r_initialized != 0)
  {
    Serial.printf("DSP is already initialized\n");
    return;
  }
  if (dsps_fft4r_init_fc32(NULL, nfft / 2) != ESP_OK)
  {
    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;
  }

  dsps_fft4r_fc32(data, nfft / 2);
  dsps_bit_rev4r_fc32(data, nfft / 2);
  dsps_cplx2real_fc32(data, nfft / 2);
    
  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]);
  }
}

これを実行すると、細かい誤差はあれ、PCでの実数FFTと同等の結果が得られた。

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

じつは Kiss FFT も！

Kiss FFTのソースを読むと、じつはKiss FFTも内部ではESP32と同じく、実数列を半分サイズの複素FFTにつっこんで結果を換算するという処理をしてるっぽい。数学的背景がいまひとつよく分からないが、とりあえず計算結果が正しいことは確認できたし、イディオムだと思っておくことにする。

若い頃に数学を真面目にやらなかった報いであるが、しかし今はもう泥縄的に数学を深掘りできるほど若くはないのである。。。