タコさんブログ

プログラミングメモと小言

Swift vDSP で 実FFT

C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理 ページ33 で取り上げられているサイン波を vDSP(vectorized digital signal processing) ライブラリを使用して、FFTにより振幅スペクトルを求める(ここでは簡単のために窓関数を用いずにFFTを実行する)。

サイン波は次式のように表せられる。

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)

環境

準備

Accelerateフレームワークをインポート。

import Accelerate

サンプルデータの用意

振幅(Amp)を 0.25、 基本周波数(f0)を 250 Hz、 標本化周波数(fs)を 8k Hz、 標本化数を 64 としたときのサイン波(式1)のデータの配列 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)
}

サンプルデータをDSPSplitComple型に変換する

まず、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の第1引数 __Log2n には、底を2とする標本化数の対数(i.e. log2(number_of_samples))以上を指定する。第2引数 __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成分ナイキスト成分 がそれぞれ入っている。音データの場合は無視できるので、インデックスは1から出力する。

// 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言語ではじめる音のプログラミング―サウンドエフェクトの信号処理 に記載されている値が出力されていることがわかる。

参考URL