iOSで高速フーリエ変換を使う

画像処理や音声処理でお世話になるFFT（高速フーリエ変換）ですが、iOSの場合、AppleよりFFTを計算するライブラリが用意されています。
その使い方について。

準備

プロジェクトにAccelerate.frameworkを追加してください。
そして、FFTを用いたいソースコードに以下を追記してください。
例：

FFTSample.m

# import <Accelerate/Accelerate.h>

FFTを計算する

以下が実際の計算になります。

FFTSample.m

// data: 入力
// length: 入力のうち利用する要素数
- (void)calcFFT:(float*)data dataLength:(int)length
{
    // データ長を2のn乗の大きさにする
    unsigned int sizeLog2 = (int)(log(length)/log(2));
    unsigned int size = 1 << sizeLog2;
        
    // fftのセットアップ
    FFTSetup fftSetUp = vDSP_create_fftsetup(sizeLog2 + 1, FFT_RADIX2);
    
    // 窓関数の用意
    float* window = calloc(size, sizeof(float));
    float* windowedInput = calloc(size, sizeof(float));
    vDSP_hann_window(window, size, 0);
    
    // 窓関数を入力値に適用し、windewedInputへ
    vDSP_vmul(data, 1, window, 1, windowedInput, 1, size);
    
    // 入力を複素数にする
    DSPSplitComplex splitComplex;
    splitComplex.realp = calloc(size, sizeof(float));
    splitComplex.imagp = calloc(size, sizeof(float));
    
    for (int i = 0; i < size; i++) {
        splitComplex.realp[i] = windowedInput[i];
        splitComplex.imagp[i] = 0.0f;
    }

    // FFTを計算する
    vDSP_fft_zrip(fftSetUp, &splitComplex, 1, sizeLog2 + 1, FFT_FORWARD);
    
    // 結果を表示する
    // FFTの性質から半分のデータのみ利用する
    for (int i = 0; i <= size/2; i++) {
        float real = splitComplex.realp[i];
        float imag = splitComplex.imagp[i];
        float distance = sqrt(real*real + imag*imag);
        NSLog(@"[%d] %.2f", i, distance);
    }
    
    // メモリを開放する
    free(splitComplex.realp);
    free(splitComplex.imagp);
    free(window);
    free(windowedInput);
    vDSP_destroy_fftsetup(fftSetUp);
    
}

あとは実行してみましょう。
元になるデータですが、例えばsin波などをつかってみてください。

FFTSample.m

- (void)execute
{
    int dataLength = 512;
    float* data = calloc(dataLength, sizeof(float));
    for(int i = 0; i < dataLength; i++) {
        // dataLength中で30回振動させる
        data[i] = sin(i * 30.0f / (float)dataLength * 2.0f * M_PI);
    }
    [self calcFFT:data dataLength:dataLength];
}

恐らく以下のような出力を得られるかと思います。
30付近のみ掲載します。あとは0なので。

[28] 0.00
[29] 128.00
[30] 256.00
[31] 128.00
[32] 0.00

参考

Appleの公式ドキュメント
http://developer.apple.com/library/ios/#documentation/Performance/Conceptual/vDSP_Programming_Guide/UsingFourierTransforms/UsingFourierTransforms.html

vDSPを使う（高速フーリエ変換編 その6）：なんてこったいブログ
http://nantekottai.com/2010/10/04/fft-with-vdsp-6/