Swift vDSP で 実FFT
C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理
サイン波は次式のように表せられる。
x(n) = Amp * sin(2 * PI * f0 * n / fs), (0<=n<N)
ここで、Ampは振幅、f0は基本周波数、fsは標本化周波数、nは標本化数。
FFTを実行する音データとして、以下の値をとるサイン波を考える。
- 振幅(Amp) 0.25
- 基本周波数(f0) 250 Hz
- 標本化周波数(fs) 8k Hz
- 標本化数(n) 64
このとき、サイン波は次のようになる。
x(n) = 0.25 * sin(PI * n / 16.0), (0 <= n < 64) (式1)
環境
- Xcode 7.3.1
- Swift 2.2
準備
Accelerateフレームワークをインポート。
import Accelerate
サンプルデータの用意
振幅(Amp)を 0.25、 基本周波数(f0)を 250 Hz、 標本化周波数(fs)を 8k Hz、 標本化数を 64 としたときのサイン波（式１）のデータの配列 samples を用意する。
// 標本化数 let numSamples = 64 // サンプルデータ var samples = [Float](count: numSamples, repeatedValue: 0) for n in 0..<numSamples { samples[n] = 0.25 * sinf(Float(M_PI) * Float(n) / 16.0) }
サンプルデータをDSPSplitComplex型に変換する
まず、vDSPで実FFTを実行する前に、vDSP_ctoz関数を使用してDSPSplitComplex型にパック（変換）する必要がある。 vDS_ctozを 実数データの配列 A = {A[0], A[1], … A[n]} に適用すると even-odd 配列 AEvenOdd = {A[0], A[2], …, A[n-1], A[1], A[3], … A[n]} のように変換される。
var reals = [Float](count: numSamples/2, repeatedValue: 0) var imgs = [Float](count: numSamples/2, repeatedValue: 0) var splitComplex = DSPSplitComplex(realp: &reals, imagp: &imgs) let samplesPtr = UnsafePointer<DSPComplex>(samples) vDSP_ctoz(samplesPtr, 2, &splitComplex, 1, vDSP_Length(numSamples/2))
FFT ウェイト配列(FFT Weights Arrays)
FFTのウェイト配列は vDSP_create_fftsetup 関数を使用することによって生成する。 vDSP_create_fftsetupの第１引数 __Log2n には、底を２とする標本化数の対数(i.e. log2(number_of_samples))以上を指定する。第２引数 __Radix には、基数オプションを指定する。基数オプションには、ビット演算ORでFFT_RADIX2、FFT_RADIX3、FFT_RADIX5の任意の組み合わせが使用できる。
// Create FFT setup // __Log2nは log2(64) = 6 より、6 を指定 let setup = vDSP_create_fftsetup(6, FFTRadix(FFT_RADIX2))
FFTの実行・出力
vDSP_fft_zrip 関数を使用してFFTを実行する。
// Perform FFT vDSP_fft_zrip(setup, &splitComplex, 1, 6, FFTDirection(FFT_FORWARD))
vDSP Programming Guide に記載されている通り、C_imp = 2 * C_math となっている。C_math を求めるには 1/2倍する必要がある。vDSP_vsmul 関数を使用して各配列要素の値を1/2倍する。また、splitComplex.realp[0], splitComplex.imagp[0] には DC成分、ナイキスト成分 がそれぞれ入っている。音データの場合は無視できるので、インデックスは１から出力する。
// splitComplex.realp, splitComplex.imagpの各要素を1/2倍する var scale:Float = 1 / 2 vDSP_vsmul(splitComplex.realp, 1, &scale, splitComplex.realp, 1, vDSP_Length(numSamples/2)) vDSP_vsmul(splitComplex.imagp, 1, &scale, splitComplex.imagp, 1, vDSP_Length(numSamples/2)) // 複素数の実部と虚部を取得する let r = Array(UnsafeBufferPointer(start: splitComplex.realp, count: numSamples/2)) let i = Array(UnsafeBufferPointer(start: splitComplex.imagp, count: numSamples/2)) for n in 1..<numSamples/2 { let rel = r[n] let img = i[n] let mag = sqrtf(rel * rel + img * img) let log = "[%02d]: Mag: %5.2f, Rel: %5.2f, Img: %5.2f" print(String(format: log, n, mag, rel, img)) } // setupを解放 vDSP_destroy_fftsetup(setup)
出力結果は以下のようになる。
[01]: Mag: 0.00, Rel: -0.00, Img: -0.00 [02]: Mag: 8.00, Rel: 0.00, Img: -8.00 [03]: Mag: 0.00, Rel: 0.00, Img: 0.00 以下0が続く . . .
C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理
