LoginSignup
39
41

More than 5 years have passed since last update.

@Talokay(Keitaro Takahashi)

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

Last updated at Posted at 2014-06-13

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

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

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

FFTを行う手順として、

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

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

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

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

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窓

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窓

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窓でいうならば、

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()関数のコード

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波生成コードの後

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フレームサイズの半分のサイズを作ります。

main.c
    int ffthalfsize = FFTSIZE/2;

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

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

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

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

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

main.c
FFT(buffer,&splitComplex);

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

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

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

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の各ビンがどの周波数に対応するかを計算します。

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

結果を出力します。

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

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

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

FFT()

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

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の自然対数をとります。

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

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

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

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

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

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

spectram_fft02.jpg

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

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

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

spectram_fft03.jpg

となり、掛け合わせ後は

spectram_fft04.jpg

となります。

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

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にかけます。

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 の実数部と虚数部の値をスケーリングしなければなりません。

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のセットアップと窓関数のメモリ領域を解放します。

main.c
    //free 
    vDSP_destroy_fftsetup(fftsetup);
    free(window);

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

FFT_free()

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

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

おまけ　矩形波

おまけで矩形波をこのプログラムを使ってFFTにかけてみます。
矩形波の生成方法はこの記事を参照してください。

FFTをしてスペクトラルの強度をとったものをプロットすると、以下のようになりました。

spectram_fft05.jpg

ちゃんとピークが等間隔に並んでいます。上手く動いているようです。

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

39
41
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
39
41