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


# オーディオ・音声分析への道5 高速フーリエ変換 FFT

今回は、vDSPにて高速フーリエ変換(FFT)のやり方を紹介します。
vDSPにはFFTを行う為のいくつかの関数が用意されています。

この記事では、FFTを行うためのプログラムを一つにまとめた、FFT()という関数を作るところから始めたいと思います。

FFTを行う手順として、

1. FFTのフレームサイズを指定する。
2. FFTフレームサイズ分のオーディオデータを取得する。
3. 2で取得したオーディオデータに窓関数を掛け合わせる。
4. 窓関数を掛け合わせたオーディオデータをDSPSplitComplex型のデータ領域に格納する。
5. FFTを行う。
6. スケーリングを行う。（vDSP独特の仕様）

とここまでがあります。
この中で **DSPSplitComplex型** というデータタイプが現れますが、これは以下のような構造体をなしていて、信号の実数部と虚数部を格納する事ができます。

```C:main.c
typedef struct DSPSplitComplex {
    float  *realp;
    float  *imagp;
} DSPSplitComplex;
```

**DSPSplitComplex型** は **float型** のデータを扱う時に使用しますが、 **double型** の場合は **DSPDoubleSplitComplex型** を使います。
また、窓関数に関して、vDSPはいくつかの窓関数を生成する関数を持っていて、主な物だけでも
hanning窓

```C:main.c
extern void vDSP_hann_window(
    float       *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
```

hamming窓

```C:main.c
extern void vDSP_hamm_window(
    float       *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
```

blackman窓


```C:main.c
extern void vDSP_blkman_window(
    float       *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
```

 等があります。
これらは全てfloat型に対応した関数で、double型で処理を行いたい場合は、関数の最後に大文字の **D** を加えます。hanning窓でいうならば、


```C:main.c
extern void vDSP_hann_windowD(
    double      *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);

```

##main()関数のコード


```C:main.c
#define PAI 3.14159265358979323846264338 // PAI
#define SAMPLE_RATE 44100 //sampling rate
#define LENGTH 4 // 4 seconds
#define AMPLITUDE 1.0// 16bit
#define FREQUENCY 440 // frequency Hz
#define FFTSIZE 2048

int main(int argc, const char * argv[])
{
    //Libsndfile
    SNDFILE *fp;
    SF_INFO sfinfo;
    float *buffer;
    buffer = malloc(FFTSIZE*sizeof(float));

    memset(&sfinfo,0,sizeof(sfinfo));

    sfinfo.samplerate = SAMPLE_RATE;
    sfinfo.frames = SAMPLE_RATE*FFTSIZE;
    sfinfo.channels = 1;
    sfinfo.format = (SF_FORMAT_WAV | SF_FORMAT_PCM_16);

    if(!(fp = sf_open("Sawtooth.wav", SFM_WRITE, &sfinfo))){
		printf("Error : Could not open output file .\n");

		return 1;

	} 

    int i;
    for(i=0;i<FFTSIZE;i++){
        buffer[i] = AMPLITUDE * sin((float)FREQUENCY / SAMPLE_RATE * 2 * PAI * i);
    }

    int ffthalfsize = FFTSIZE/2;
    DSPSplitComplex splitComplex;

    //allocate complex (fft half size)
    splitComplex.realp = malloc(sizeof(float)*ffthalfsize);
    splitComplex.imagp = malloc(sizeof(float)*ffthalfsize);

    float *magnitude = malloc(sizeof(float)*ffthalfsize);

    FFT(buffer,&splitComplex);

    vDSP_zvabs(&splitComplex, 1, magnitude, 1, ffthalfsize);

    float freq_bins = (float)SAMPLE_RATE/FFTSIZE;

    for(i=0;i<ffthalfsize;i++)
        printf("%0.2f Hz: magnitude = %f\n",i*freq_bins,magnitude[i]);
    FFT_free(&splitComplex);

    if(sf_write_float(fp, buffer,sfinfo.channels *FFTSIZE)!= sfinfo.channels *FFTSIZE)
		puts(sf_strerror (fp));

    free(buffer);	
    sf_close(fp);

    return 0;

}
```

 ここで注目するのはsin波生成コードの後


```C:main.c
    int ffthalfsize = FFTSIZE/2;
    DSPSplitComplex splitComplex;

    //allocate complex (fft half size)
    splitComplex.realp = malloc(sizeof(float)*ffthalfsize);
    splitComplex.imagp = malloc(sizeof(float)*ffthalfsize);  

    float *magnitude = malloc(sizeof(float)*ffthalfsize);

    FFT(buffer,&splitComplex);

    vDSP_zvabs(&splitComplex, 1, magnitude, 1, ffthalfsize);

    float freq_bins = (float)SAMPLE_RATE/FFTSIZE;

    for(i=0;i<ffthalfsize;i++)
        printf("%0.2f Hz: magnitude = %f\n",i*freq_bins,magnitude[i]);

    FFT_free(&splitComplex);
```

の部分です。ちなみにFFT = 2048です。

FFTフレームサイズの半分のサイズを作ります。

```C:main.c
    int ffthalfsize = FFTSIZE/2;
```

 次に **DSPSplitComplex型** のデータを宣言し、実数部と虚数部のデータ領域をallocateします。サイズはFFTサイズの半分になります。

```C:main.c
    splitComplex.realp = malloc(sizeof(float)*ffthalfsize);
    splitComplex.imagp = malloc(sizeof(float)*ffthalfsize);
```

次に、取得するスペクトラルの強度を格納する変数を作ります。サイズはFFTサイズの半分になります。

```C:main.c
     float *magnitude = malloc(sizeof(float)*ffthalfsize);
```

そしていよいよFFTを実行します。この関数はmain()の説明の後に紹介します。
FFT()は引数として、オーディオデータが入った **buffer** と **splitComplex** のポインタをとります。

```C:main.c
FFT(buffer,&splitComplex);
```

**FFT()** を実行すると、 **splitComplex** の実数部と虚数部にはそれぞれFFT変換された値が入ります。
よって、この実数と虚数の値のパワーをとります。

```C:main.c
vDSP_zvabs(&splitComplex, 1, magnitude, 1, ffthalfsize);
```

**vDSP_zvabs()** は、sqrt(real * real + img * img)の処理を行い、結果を **magnitude** に格納します。

```C:main.c
extern void vDSP_zvabs(
    const DSPSplitComplex *__vDSP_A,
    vDSP_Stride            __vDSP_IA,
    float                 *__vDSP_C,
    vDSP_Stride            __vDSP_IC,
    vDSP_Length            __vDSP_N)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
```

**Stride** はそれぞれ1になります。
vDSP_Lengthは計算するデータ量なので、 **ffthalfsize** です。

FFTの各ビンがどの周波数に対応するかを計算します。

```C:main.c
float freq_bins = (float)SAMPLE_RATE/FFTSIZE;
```

結果を出力します。

```C:main.c
for(i=0;i<ffthalfsize;i++)
        printf("%0.2f Hz: magnitude = %f\n",i*freq_bins,magnitude[i]);

```

結果の一部を抜粋します。

```C:result
0.00 Hz: magnitude = 0.000018
21.53 Hz: magnitude = 0.000018
43.07 Hz: magnitude = 0.000019
64.60 Hz: magnitude = 0.000020
86.13 Hz: magnitude = 0.000023
107.67 Hz: magnitude = 0.000026
129.20 Hz: magnitude = 0.000030
150.73 Hz: magnitude = 0.000036
172.27 Hz: magnitude = 0.000044
193.80 Hz: magnitude = 0.000055
215.33 Hz: magnitude = 0.000072
236.87 Hz: magnitude = 0.000096
258.40 Hz: magnitude = 0.000134
279.93 Hz: magnitude = 0.000195
301.46 Hz: magnitude = 0.000301
323.00 Hz: magnitude = 0.000504
344.53 Hz: magnitude = 0.000943
366.06 Hz: magnitude = 0.002103
387.60 Hz: magnitude = 0.006500
409.13 Hz: magnitude = 0.051471
430.66 Hz: magnitude = 0.221126
452.20 Hz: magnitude = 0.202370
473.73 Hz: magnitude = 0.034186
495.26 Hz: magnitude = 0.005429
516.80 Hz: magnitude = 0.001862
538.33 Hz: magnitude = 0.000858
559.86 Hz: magnitude = 0.000466
581.40 Hz: magnitude = 0.000281
602.93 Hz: magnitude = 0.000182
624.46 Hz: magnitude = 0.000125

```

ここで、解析したデータは440Hzのsin波でした。
結果は、430hz付近と452hz付近のFFT binが、高い強度値を持つことが分かります。よって430hzと452hzの間の周波数つまり440hzに近い周波数にもっとも大きな強度がある事を示しています。スペクトラルをプロットすると以下の様な画像になっています。

![spectram_fft01.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/12cfccd9-d64c-4574-3304-c81291a59eec.jpeg)

##FFT()

次はFFT()の部分を解説します。
コードは以下の通りです。

```C: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);

}
```

**FFT()** は、FFTフレーム分の長さを持ったfloat型のオーディオデータと、 **DSPSplitComplex型** のデータをを引数としてとります。

FFTサイズの半分を作り、 **FFTSIZE** の2の自然対数をとります。

```C:main.c
int ffthalfsize = FFTSIZE/2;
vDSP_Length log2n = log2(FFTSIZE);
```

 次に、vDSPでFFTを行うためのsetupを行います。

```C:main.c
//setup fft
 FFTSetup fftsetup = vDSP_create_fftsetup(log2n, kFFTRadix2);
```

第一引数には **FFTSIZE** の2の自然対数を、第二引数には **kFFTRadix2 = 0** を代入し、モードを設定します。
次に、窓関数を生成します。今回はhanning窓を使用します。


```C:main.c
float *window = malloc(sizeof(float)*FFTSIZE);
vDSP_hann_window(window, FFTSIZE, 0); 
```

hanning窓をプロットすると以下のような形になります。

![spectram_fft02.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/8f1f9361-5336-79ef-62ad-35f80d4df986.jpeg)

生成したhanning窓をオーディオデータに掛け合わせます。

```C:main.c
     vDSP_vmul(wavedata, 1, window, 1, wavedata, 1, FFTSIZE);
```

こうすることで、窓がけが完了しました。
結果をプロットしてみると、窓関数掛け合わせ前のオーディオシグナルは、

![spectram_fft03.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/8a733e99-84d0-eb87-247e-a4e7af1e91e8.jpeg)

となり、掛け合わせ後は

![spectram_fft04.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/f993eb15-3776-65da-a4ee-9456c8e8a4cf.jpeg)

となります。

さて、いよいよFFTをしましょう。
まず、窓関数を掛け合わせたオーディオシグナルを、 **COMPLEX型** に変換します。
これは以下の関数にて行われます。

```C:main.c
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, splitComplex, 1, ffthalfsize );
```

この場合の **Stride** ですが、wave dataは **FFTSIZE**分あるのに対して、 **splitComplex** は **ffthalfsize** 分のデータ領域しかありません。これは、vDSPの仕様で、計算量、メモリ使用量を減らす事を目的にしています。
よって、処理するデータ量も **ffthalfsize** （第五引数）になります。この場合、wave dataの **Stride** は2、 **splitComplex** の **Stride** は1となります。
次にsplitComplexをFFTにかけます。

```C:main.c
     vDSP_fft_zrip(fftsetup, splitComplex, 1, log2n, FFT_FORWARD);
```

先ほど **vDSP_create_fftsetup()** にて生成したFFTのセットアップと、 **splitComplex** , **splitComplex** の **stride** 、**FFTSIZE** の2の自然対数、そして、FFTを行う場合は **FFT_FORWARD** をセットします。iFFTを行う場合は、ここに **FFT_INVERSE** をセットします。

これでFFTの処理は終わりますが、vDSPの仕様で、出力された **splitComplex** の実数部と虚数部の値をスケーリングしなければなりません。

```C:main.c
    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);
```

以上のように **1.0/(FFTSIZE*2)** 分を全ての値に掛け合わせます。

最後に、使い終わったFFTのセットアップと窓関数のメモリ領域を解放します。
 

```C:main.c
    //free 
    vDSP_destroy_fftsetup(fftsetup);
    free(window);
```

これでFFT()の解説は終わります。

##FFT_free()						

splitComplexのデータ領域を解放します。

```C:main.c
void FFT_free(DSPSplitComplex *splitComplex)
{
    free(splitComplex->realp);
    free(splitComplex->imagp);
}
```

##おまけ　矩形波	

おまけで矩形波をこのプログラムを使ってFFTにかけてみます。
矩形波の生成方法は[この記事](http://qiita.com/programanx1/items/212453f55e01b44e303f)を参照してください。

FFTをしてスペクトラルの強度をとったものをプロットすると、以下のようになりました。
	
![spectram_fft05.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/77336775-2682-9f5c-636b-04b9a23a8171.jpeg)
				
ちゃんとピークが等間隔に並んでいます。上手く動いているようです。

 
次回は、iFFTを行って、元のオーディオデータに戻したり、その間にフィルターなどをかまして音を変調したりしてみたいと思います。