[vDSP][信号処理]オーディオ・音声分析への道5.1 Swiftで高速フーリエ変換 FFT

Swiftが発表されてから随分日にちが経ちました。
Objective-Cを結局覚えなかった自分が、これを機にSwiftを勉強してみようと思い立ち早二ヶ月。以前書いた記事の内容をSwiftで書き直してみたいという欲求が湧いてきました。
今回はSwiftでvDSPを取り扱う際に、色々とつまづいた点や解決策を、他の人たちが行った方法などを交えて紹介したいと思います。

(注意点)筆者は、Swift初心者です。しかも、Objective-Cを飛ばし、C言語からいきなりSwiftに入門した経歴を持っています。SwiftをC言語と比較して説明してしまう時に、適切でない表現を行ってしまう事もあるかと思います。ご了承くださいませ。また、誤字脱字などは毎度毎度ご容赦ください。

(注意点2 : 2016年12月23日追記)
Swift3での実装はポインタの扱いが異なります。以下の記事も参考にしてください。
http://qiita.com/programanx1/items/7237344dac049a4b95a4

今回はデータ型にまつわる、つまづき所を紹介したいと思います。

厳格な型指定

Swiftでは変数の型指定が非常に厳格になっており、例えばFloat型とDouble型を同時に扱う場合は、どちらかの型にキャストして、データ型を統一してやる必要があります。

例

sample.c

var a:Float = 10.0
var b:Double = 2.56

var c:Float = a + Float(b)

さて、vDSPのライブラリ関数は、独自のデータ型の引数を使用します。
vDSP_LengthやComplex型などなどです。
これを厳密に取り扱う必要があります。
また、vDSPではスカラー計算を行う事が多く、配列やポインタの取り扱いに十分注意する必要があります。

ArrayとUnsafePointer

Swiftにも配列とポインタの概念があります。

配列は[データ型]で表します。例えば、Float型であるならば、
[Float]
となります。
ポインタは可変可能であれば、
UnsafeMutablePointer
可変不可ならば (C言語でいうところの、const)
UnsafePointer
です。

Arrayの宣言、サイズの指定は以下のとおりです。

Array.c

var array = [Float](count:100,repeatedValue:0.0)

上の例では、aという100個のFloat型の配列を宣言すると同時に、100個データすべてに0.0を代入し初期化しています。

では、ポインタで同じようにデータをallocateするにはどうするか、というと、以下のようになります。

Pointer.c

var p = UnsafeMutablePointer<Float>.alloc(100)

また、既存の配列を割り当てる場合は

Pointer2.C

var array = [Float](count:100,repeatedValue:0.0)
var p = UnsafeMutablePointer<Float>(array)

上記の通りにします。これは、配列をポインタに変換する方法になります。
このテクニックは、DSPComplex型を取り扱う際に必須になります。

具体例

Swiftでは、データ型以外にも、ポインタの可変可能、不可の違いも明確に指定する必要があります。

vDSP_vmul、vDSP_vsmul

まず初めに、vDSP_vmulと、vDSP_vsmulを例にします。
vDSP_vmulは配列同士の掛け算を行います。
また、vDSP_vsmulはscalarでの掛け算を行います。
C言語では以下のように使いました。

vDSP_vmul.C

// 全てのデータに同じ数を掛け合わせる場合
int scale = 10;
vDSP_vsmul(src, 1, &scale, 1, dist, 1, 100);

//配列同士
vDSP_vmul(src1, 1, src2, 1, dist, 1, FFTData.fftsize);

Swiftでは以下の通りです。(Playgroundに書いたコードをそのまま載せます)

vDSP_vmul.C

//PlaygroundでのAccelerateのimport法
import Accelerate

var array1 = [Float](count:10,repeatedValue:1.0)
var array2 = [Float](count:10,repeatedValue:10.0)
var dist = [Float](count:10,repeatedValue:0.0)


// 全てのデータに同じ数を掛け合わせる場合(scalar)
let scale:Float = 100.0
vDSP_vsmul(array1, 1, [scale], &dist, 1, 10)

//配列同士
vDSP_vmul(array1,1,array2,1,&dist,1,10)

C言語との差異は、scaleをUnsafePointerに変換するため、[]で囲います。また、distをUnsafeMutablePointerとするため、&をつけます。

COMPLEX型

今回苦労したのが、このCOMPLEX型です。
COMPLEX型は、

complex.c

DSPComplex // float
DSPDoubleComplex // double

などがあり、それぞれのデータ構造は以下のとおりです。

complex.c

struct DSPComplex {
    var real: Float
    var imag: Float
}
struct DSPDoubleComplex {
    var real: Double
    var imag: Double
}

vDSPを使用したFFTの手順は、[vDSP][信号処理]オーディオ・音声分析への道5 高速フーリエ変換 FFTで紹介した通りです。
FFTのコードをここに再度載せます。

*注意 Swift 3 ではポインタの扱いが変更になりました。
こちらに随時まとめていきますので、参照ください。
[vDSP][Swift3] ポインタの扱い
http://qiita.com/programanx1/items/7237344dac049a4b95a4

main.c

void FFT(float *wavedata, DSPSplitComplex *splitComplex)

{   
    int ffthalfsize = FFTSIZE/2;
    vDSP_Length log2n = log2(FFTSIZE);

    //setup fft
    FFTSetup fftsetup = vDSP_create_fftsetup(log2n, kFFTRadix2);

    //make window (fft size)
    float *window = malloc(sizeof(float)*FFTSIZE);
    //hanning window
    vDSP_hann_window(window, FFTSIZE, 0);
    //windowing
    vDSP_vmul(wavedata, 1, window, 1, wavedata, 1, FFTSIZE);


    //transform to complex
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, splitComplex, 1, ffthalfsize );

    //FFT forward
    vDSP_fft_zrip(fftsetup, splitComplex, 1, log2n, FFT_FORWARD);

    //scaling
    float scale = 1.0/(FFTSIZE*2);
    vDSP_vsmul(splitComplex->realp, 1, &scale, splitComplex->realp, 1, ffthalfsize);
    vDSP_vsmul(splitComplex->imagp, 1, &scale, splitComplex->imagp, 1, ffthalfsize);

    //free
    vDSP_destroy_fftsetup(fftsetup);
    free(window);

}

ここで、

main.c

//transform to complex
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, splitComplex, 1, ffthalfsize );

この一行が曲者となります。
vDSP_ctozを使用せず、Complex型の実部にオーディオデータをそのまま代入すると、出力結果が上手くいかないのがvDSPの厄介なところです。色々と試しましたが、この関数を避けては通れないようです。

ちなみに、Mattt ThompsonさんのvDSPライブラリSurgeでは、この処理を行わず、以下のようにcomplex型の実部に直接オーディオデータを代入しています。

surge

// 実部 inputはオーディオデータなど
var real = [Float](input)
//虚部 虚部は0で初期化されている
var imaginary = [Float](count: input.count, repeatedValue: 0.0)
//DSPSplitComplexの宣言と初期化
var splitComplex = DSPSplitComplex(realp: &real, imagp: &imaginary)

筆者はこの動作は未確認です。が、C言語で筆者が四苦八苦していた時に行った実験では、以下の部分を

main

vDSP_Length log2n = log2(FFTSIZE);

main

vDSP_Length log2n = log2(FFTSIZE*2);

と、FFTSIZEの二倍分与えないと結果が正しくなりませんでした。
この際、iFFTを行う場合に都合が悪くなります。

ですので、

main.c

//transform to complex
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, splitComplex, 1, ffthalfsize );

としてオーディオデータをsplitComplexに代入しますが、
Swiftでは、このまま書いても、勿論上手くいきません。
この問題を解決するには、Christopher LiscioさんのSMUGMath in Swiftが大変参考になりました。

Christopher Liscioさんの手法を筆者のプログラムに取り入れると以下のようになります。

main.c

//wavedataをComplex*型に変換
var audioBufferComplex = UnsafePointer<DSPComplex>(wavedata)

//Complex型に変換したものをvDSP_ctozに使用する
vDSP_ctoz(audioBufferComplex, 2, &splitComplex, 1, FFTData.ffthalfsize)

//FFT
vDSP_fft_zrip(setup,&splitComplex,1,log2n,FFTDirection(FFT_FORWARD))

これで上手くいきました。
今回はここまで。
Swiftで書き直したFFT全体のコードはまた気が向いた時に載せます。
次回はiFFTなどのSwiftでの実装について。