vDSPを使う（高速フーリエ変換編その2）

いよいよvDSPを使って高速フーリエ変換してみます。※この記事はObjective-Audioさんのこちらの記事を参考に書いています。

vDSPの設定

加減乗除の時同様、vDSPを使えるように設定します。設定の仕方は加減乗除の記事を参照してください。基本的にはAccelerate.frameworkをリンクして、<Accelerate/Accelerate.h>をインポートするだけです。

解析する信号の準備

フーリエ変換させたい信号を用意します。ファイルやマイクから読み込んでもいいですが、長くなるので今回はサインウェーブを適当に合成して作ります。
信号の長さは前回の記事で触れたように2のn乗になっている必要があるので、2の9乗(512)にします。

float inputWave[512];
float phase = 0.0f;
float freq = 30.0f;
for (int i=0; i<512; i++) {
	phase += freq / 512.0f * 2.0f * M_PI; // 2 * PIで一周になります
	inputWave[i] = sin(phase) + sin(phase * 2.0f) * 0.5f;
}

512という長さの中で30回振動するサイン波と、その倍音成分を半分の振幅で足しています。サンプリング周波数を512Hzと考えれば、30Hzのサイン波と倍音成分(60Hzのサイン波)が含まれる波形になります。
本題とは離れますが、離散フーリエ変換ではここらへんの単位の計算が意外とこんがらがります。実際の音声信号は大抵44100Hz等でサンプリングされているので、512サンプルは 512 / 44100 = 0.012秒となり、0.012秒で30回振動するということは実際には上記の波形は2500Hzのサイン波とその倍音が含まれていることになります。

複素数への変換

ここが最初にわかりにくい部分だと思うのですが、フーリエ変換させたい信号は複素数の配列として渡してあげなければいけません。一次元の信号である音声信号をどうやって複素数の信号に変換するのか、という話ですが、ただ実数部分に音声信号をいれるだけです。虚数部分は0をいれておきます。
vDSPでは複素数を扱うためにDSPSplitComplexという構造体が用意されているのでこれを使います。この構造体には実数部分を示すrealpというポインタと虚数部分を示すimagpというポインタがあり、そこに配列を入れられるようになっています。その分のメモリは自分で確保してあげなければいけないので、mallocやcallocを使って確保します。今回は512サンプル分確保します。虚数部分は代入することはありませんが、メモリは確保してあげなければいけません。メモリを確保したら代入します。

DSPSplitComplex splitComplex;
splitComplex.realp = calloc(512, sizeof(float));
splitComplex.imagp = calloc(512, sizeof(float));
for (int i=0; i<512; i++) {
	splitComplex.realp[i] = inputWave[i];
}

フーリエ変換の準備

vDSPでフーリエ変換を行うための設定をします。とはいっても、設定しなければいけないのは配列の大きさと基数と呼ばれるものだけです。
配列の大きさは上で決めた512ですが、ここでは2のn乗の「n」の部分を指定する形になるので「9」を指定します。基数はFFTのアルゴリズムを決定する上で使用されるものですが、ここでは一般的な2の基数のFFTであるFFT_RADIX2を指定しますが、通常ここは気にしなくていいです。
フーリエ変換の設定は、実際にフーリエ変換を行う前にvDSP_create_fftsetup関数を使って設定し、フーリエ変換が終わった後にvDSP_destroy_fftsetupを呼んで破棄します。

FFTSetup fftSetup = vDSP_create_fftsetup(9, FFT_RADIX2);
//TODO: ここでフーリエ変換するよ
vDSP_destroy_fftsetup(fftSetup);

フーリエ変換の実行

変換結果の返し方

vDSPにはプログラム側のニーズにあわせていくつかFFT用の関数が用意されています。それらの関数の大きな違いの一つとして、変換結果をどのように返すか、という違いがあります。
フーリエ変換の入力は複素数の配列でしたが、変換後の出力も同様に複素数の配列となります。新しく複素数の配列をつくってそこに結果を代入してもらうこともできます。ただ、毎回配列を作り直していたら計算量が増えてしまうのでもったいないです。そこで、入力信号につかった配列にそのまま変換結果をいれて返してもらう事もできます。今回は面倒なので後者を使います。

変換の方向

前の記事に書いた通り、フーリエ変換は音声信号（グラフにしたとき横軸が時間）⇔周波数特性を相互に変換できます。vDSPでは音声信号→周波数特性に変換する方向がFFT_FORWARD、反対がFFT_INVERSEと定義されています。今回は信号から周波数特性を得たいので前者です。

FFTSetup fftSetup = vDSP_create_fftsetup(9, FFT_RADIX2);
vDSP_fft_zrip(fftSetup, &splitComplex, 1, 9, FFT_FORWARD);
vDSP_destroy_fftsetup(fftSetup);

vDSPには入出力の複素数の配列を使い回してFFTをやってくれる関数がいくつかあるのですが、今回のように一次の信号(realしか使っていないない入力)で使うのはvDSP_fft_zripになります。引数で、

FFTの設定
配列のアドレス
ストライド(加減乗除の記事を参照)
大きさ(2のn乗のnの部分)
方向

を指定します。
最後に不要になった複素数の配列用のメモリを解放しておきます。

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

変換結果の利用

フーリエ変換自体は以上でできますが、この変換結果をどう利用するのかという部分については長くなったのでまた次の記事で書きます。

参考リンク

Pocket