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種類)です。

vDSPの基本的な使い方は、こちらを参考に。

Vectorizing with vDSP and vecLib (Forever For Now)

英文や数式が苦手な僕ですが、こうやってコードっぽく書いてあると分かりやすいです。

Mac OS XのAccelerateフレームワークには、連続したデータの配列にデジタル信号処理をするときなどに便利な関数がいろいろ用意されています。Accelerateフレームワークを使えば１つの記述をするだけで、内部で勝手に判断してAltiVecやSSE等で素早く処理を行ってくれるようです。場合によってはアクセラレートしない事もありますが、いちいち自分でfor文を回したりしなくても良くて、処理が速くなる可能性があるのですから、使わない手はありません。

では、Accelerateフレームワークの中にあるvDSP_deq22()という関数でbiquadのフィルタをかけるクラスを作ってみます。vDSP_deq22()関数は以下のように宣言されています。

void vDSP_deq22 (
   float * A,		//インプットする配列
   vDSP_Stride I,	//インプットのストライド
   float * B,		//係数の配列
   float * C,		//アウトプットされる配列
   vDSP_Stride K,	//アウトプットのストライド
   vDSP_Length N	//処理するサイズ
);

これはfloat配列用ですが、double配列用にvDSP_deq22D()という関数も用意されています。基本的にvDSPの関数はDがついているのがdouble用で、ついてないのがfloat用という名前のつけ方になっているようです。

まず、Accelerateフレームワークを追加しインポートします。

#import <Accelerate/Accelerate.h>

オーディオデータを細切れに分けて処理をすることを考え、使い回すデータのインスタンス変数を用意します。

@interface IIRUnit : NSObject {
    float *gCoefBuffer;	//係数の配列
    float *gInputKeepBuffer;	//インプット2サンプル分のバッファ
    float *gOutputKeepBuffer;	//アウトプット2サンプル分のバッファ
}

- (id)init
{
    self = [super init];
    if (self != nil) {
        gInputKeepBuffer = calloc(2, sizeof(float));
        gOutputKeepBuffer = calloc(2, sizeof(float));
        gCoefBuffer = calloc(5, sizeof(float));
    }
    return self;
}

- (void)dealloc
{
    free(gInputKeepBuffer);
    free(gOutputKeepBuffer);
    free(gCoefBuffer);
    [super dealloc];
}

biquadフィルタでは処理するデータに対して２サンプル前と１サンプル前のデータが必要なので、その配列を用意します。インプット用２サンプル分とアウトプット用２サンプル分です。
次に５個の係数の配列を用意します。係数の順番はB0、B1、B2、A1、A2です。６個の係数を使用する場合は先に５つの係数をA0で割っておきます。フィルタやEQの係数の求め方はAudio-EQ-Cookbook等を参考にしてください。

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

それでは、処理を行うメソッドです。

- (void)processWithIoData:(float *)ioData frames:(NSUInteger)frames
{
    //処理用のバッファを用意する
    float *tInputBuffer = malloc((frames + 2) * sizeof(float));
    float *tOutputBuffer = malloc((frames + 2) * sizeof(float));
    
    //処理用のバッファにデータをコピー
    memcpy(tInputBuffer, gInputKeepBuffer, 2 * sizeof(float));
    memcpy(tOutputBuffer, gOutputKeepBuffer, 2 * sizeof(float));
    memcpy(&(tInputBuffer[2]), ioData, frames * sizeof(float));
    
    //処理を行う
    vDSP_deq22(tInputBuffer, 1, gCoefBuffer, tOutputBuffer, 1, frames);
    
    //処理後のデータをコピー
    memcpy(ioData, tOutputBuffer, frames * sizeof(float));
    memcpy(gInputKeepBuffer, &(tInputBuffer[frames]), 2 * sizeof(float));
    memcpy(gOutputKeepBuffer, &(tOutputBuffer[frames]), 2 * sizeof(float));
    
    free(tInputBuffer);
    free(tOutputBuffer);
}

処理したいシグナルデータの配列のサイズ＋２サンプル分のサイズの配列をインプットA用とアウトプットC用に用意します。それぞれ頭の２サンプルには、前回処理した最後の２サンプルをコピーしておきます。３サンプル目以降に元のデータをコピーします。

ストライドI、Kはデータをいくつ飛ばしで読み込みor書き込みをするかを指定します。シグナルデータがNotInterleavedなモノラルの連続したデータであれば、ストライドは1です。もし、StereoのInterleavedなデータ（LRLRLR…の順番）の片方のチャンネルを処理するのであれば、ストライドは2です。StereoのRchとかならストライドを2にしつつ、vDSP_deq22関数に渡す配列のポインタをRchの先頭位置に持っていかなければいけません。

処理サイズNは、２サンプル足さない元データのサイズ（バイト数ではなくサンプル数）を指定します。

vDSP_deq22()関数で処理したら、処理前のデータA、処理後のデータCのそれぞれ最後の２サンプルを次の処理用にインスタンス変数にコピーしておきます。

処理後のデータCの３サンプル目から必要に応じてコピーします。今回はもとのデータを上書きしています。