vDSPで高速フーリエ変換を行う関数の使い方です。

vDSP One-Dimensional Fast Fourier Transforms Referenceというリファレンスを見ると、一次元のフーリエ変換だけで30個ほどの関数が用意されています。

大きく分けてRealとComplexがあり、そのそれぞれがさらに、In-Place(変換前のバッファにそのまま上書き)かOut-of-Place(変換前と後で別々のバッファ)かで分かれています。さらに細かく見るとfloat用とdouble用とか、一度にたくさん処理するものとか、いろいろあります。

今回はIn-PlaceでComplexのfft_zip()を使ってみます。

ところで、このfft_zip()関数などは、リファレンスには関数名にvDSP_がついていますが、コード補完では引数付きの状態では出てこなくて、vDSP_がつかないfft_zipで出てきたりします。

//fft_zip()関数

void vDSP_fft_zip (FFTSetup setup,	//FFTセットアップ
   DSPSplitComplex * ioData,		//複素数配列の構造体
   vDSP_Stride stride,			//ストライド
   vDSP_Length log2n,			//FFTサイズ
   FFTDirection direction);		//正変換か逆変換か

だいたいのvDSPの関数は、その関数のみで完結しますが、FFTを使うにはvDSP Libraryというリファレンスでもわざわざ使い方が解説されているように、ちょっとした準備と後始末が必要です。

まず、Accelerateフレームワークをインポートしておきます。

#import <Accelerate/Accelerate.h>

create_fftsetup()という関数で、FFTSetupというインスタンスのようなものを作成します。第一引数で指定するFFTのサイズは2の何乗か、になります。9なら512サンプル、10なら1024サンプルといった感じです。使う関数によって指定できるサイズの範囲が決まっています。

第二引数ではradixを指定します。今回はFFT_RADIX2を指定します。ほかにFFT_RADIX3とFFT_RADIX5もあり、FFTするサンプル数が合えば（２の累乗の３倍とか５倍とか）、こちらを使う事もできるようです。

//FFTするサイズを2の何乗かで指定する。この場合1024サンプル。
vDSP_Length log2n = 10;

//FFTセットアップを作成する
FFTSetup fftSetup = create_fftsetup(log2n, FFT_RADIX2);

複素数を実数と虚数を別々の配列でメンバに持つDSPSplitComplexという構造体を宣言し、realp(実数配列)とimagp(虚数配列)に、それぞれ変換するデータを入れたメモリ領域を割り当てます。

//FFTを行うデータを作成する
vDSP_Length fftSize = 1 << log2n;
DSPSplitComplex splitComplex;
splitComplex.realp = calloc(fftSize, sizeof(float));
splitComplex.imagp = calloc(fftSize, sizeof(float));
	
/*
ここでsplitComplex.realpとimagpに、変換するデータをコピーする
 */

変換する方向を指定します。正変換ならFFT_FORWARD、逆変換ならFFT_INVERSEです。ストライドも指定します。ここでは連続したデータと想定して1にしておきます。

//FFTか逆FFTを指定する。逆ならFFT_INVERSE
FFTDirection direction = FFT_FORWARD;
//ストライドを指定する
vDSP_Stride signalStride = 1;

fft_zip()関数で変換を行います。

//FFTを行う
fft_zip(fftSetup, &splitComplex, signalStride, log2n, direction);

なお、fft_zip()関数において逆変換をした場合には、FFTサイズ倍の大きさのレベルそのままでデータが返ってきていますので、FFTサイズで割ります。

//逆FFTの場合FFTサイズでシグナルを割る
if (direction == FFT_INVERSE) {
    float scale = 1.0 / fftSize;
    vsmul(splitComplex.realp, 1, &scale, splitComplex.realp, 1, fftSize);
    vsmul(splitComplex.imagp, 1, &scale, splitComplex.imagp, 1, fftSize);
}
	
/*
ここでsplitComplexからデータを取り出す
 */

このvsmul()という関数は、float配列の全要素に対して、同じ数をかけ算する関数です。このようにIn-Placeでも使えますし、アウト側に別の配列を渡してOut-of-Placeで使う事もできます。

最後にdestroy_fftsetup()でFFTセットアップを解放します。

//FFTセットアップを解放する
destroy_fftsetup(fftSetup);

free(splitComplex.realp);
free(splitComplex.imagp);

続けて変換を行うなら、fft_zipを繰り返して終わった後で解放すれば良いと思います。FFTをクラスとして作ってしまうなら、initでcreateして、deallocでdestroyという感じでしょうか。

と、単純にフーリエ変換はこんなところですが、フーリエ変換時によく使われる窓関数をかける関数もvDSPに用意されています。

vDSP_blkman_window　（ブラックマン窓）
vDSP_hamm_window　（ハミング窓）
vDSP_hann_window　（ハニング窓）

という３種類(double用も含めると6種類)です。