LoginSignup
34
28

More than 5 years have passed since last update.

@Talokay(Keitaro Takahashi)

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

Last updated at Posted at 2014-06-14

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

今回はFFTによって得られたデータにFiltering処理を施して、iFFTを行い、結果をwavデータに出力したいと思います。

手順としては、

  1. wavデータを読み込む
  2. FFTフレーム分を取り出す
  3. FFT
  4. イコライジング
  5. iFFT
  6. wavデータとして書き出す

以上です。

C言語でのwavデータ取り扱いについては、Libsndfileを使います。
Xcodeへの導入と使い方に関しては、この記事を参照してください。

以前の記事では、wavデータの書き出しは紹介しましたが、読み込みには触れていませんでした。今回はそれらもまとめてやりたいと思います。

まず2秒間のホワイトノイズを録音したwavデータを用意します。
https://soundcloud.com/keitaro-takahashi/input

既存関数のライブラリ化

前回のFFTのソースコードのままだと、使い勝手が悪いので、FFT関数とLibsndfile関数を別ファイルにまとめ、ライブラリ化ました。

FFT.h
#ifndef vDSP_01_FFT_h
#define vDSP_01_FFT_h

#include<stdio.h>
#include<Accelerate/Accelerate.h>

typedef struct {
    FFTSetup fftsetup;
    vDSP_Length log2n;
    DSPSplitComplex splitComplex;

    int fftsize;
    int ffthalfsize;
}fft_str;

void FFT_init(fft_str *fft, int fftsize);
void FFT(float *wavedata, fft_str *fft);
void iFFT(float *wavedata, fft_str *fft);
void FFT_free(fft_str *fft);

#endif
FFT.c
#import "FFT.h"

void FFT_init(fft_str *fft, int fftsize){

    fft->fftsize = fftsize;
    fft->ffthalfsize = fftsize/2;
    fft->log2n = log2(fft->fftsize);

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

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

}

void FFT(float *wavedata, fft_str *fft)
{

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

    //transform to complex
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, &fft->splitComplex, 1, fft->ffthalfsize );
    //FFT forward
    vDSP_fft_zrip(fft->fftsetup, &fft->splitComplex, 1, fft->log2n, FFT_FORWARD);


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


    //free
    free(window);

}

void iFFT(float *wavedata, fft_str *fft)
{
    vDSP_fft_zrip(fft->fftsetup, &fft->splitComplex, 1, fft->log2n, FFT_INVERSE);

    //transform to real
    vDSP_ztoc( &fft->splitComplex, 1, ( COMPLEX * )wavedata, 2, fft->ffthalfsize );


}

void FFT_free(fft_str *fft)
{
    vDSP_destroy_fftsetup(fft->fftsetup);
    free(fft->splitComplex.realp);
    free(fft->splitComplex.imagp);
}
Audio.h
#ifndef vDSP_01_Audio_h
#define vDSP_01_Audio_h

#include<stdio.h>
#include<stdlib.h>
#import"sndfile.h"


float *AudioFileLoader(char *filename, SF_INFO *sfinfo, float *data);
int AudioFileWriter(char *filename, SF_INFO *sfinfo, float *data);
void AudioFree(float *data);

#endif
Audio.C
#import "Audio.h"

float *AudioFileLoader(char *filename, SF_INFO *sfinfo, float *data)
{

    SNDFILE *sfr;

    //open file
    if(!(sfr = sf_open(filename, SFM_READ, sfinfo))){
        printf("Error : Could not open output file .\n");
        return 0;
    }


    // format should be Wav or Aiff
    if(sfinfo->format == 0x010002 || sfinfo->format == 0x020002){// 16bit

    }else{
        printf("Error : Could not open! file format is not WAV or AIFF\n");
        return 0;
    }

    if(sfinfo->samplerate!=44100){
        printf("Error : Samplingrate is not 44.1kHz : %d\n",sfinfo->samplerate);
        return 0;
    }

    if(sfinfo->channels != 1){
        printf("Error : Audio file is not MONO\n");
        return 0;
    }

    //allocate
    data = calloc(sfinfo->frames, sizeof(float));

    //read data
    sf_readf_float(sfr,data,sfinfo->frames);
    sf_close(sfr);
    return data; // success


}

int AudioFileWriter(char *filename, SF_INFO *sfinfo, float *data)
{

    SNDFILE *sfr;
    SF_INFO sfinfo_out = *sfinfo;


    //open file
    if(!(sfr = sf_open(filename, SFM_WRITE, &sfinfo_out))){
        printf("Error : Could not open output file .\n");
        return 0;
    }


    if(!(sf_writef_float(sfr, data, sfinfo->frames))){
        printf("Error : 0 data has been written.\n");
        return 0;
    }

    sf_close(sfr);
    return 1;

}

void AudioFree(float *data)
{
    free(data);
}

main関数

オーディオデータの読み込みは以下の通りに行う事ができるるようになりました。

main.C
    //Libsndfile
    SF_INFO sfinfo;
    //for audio data
    float *data = NULL;
    if(!(data = AudioFileLoader(AUDIO_IN,&sfinfo,data))){
        printf("Could not open Audio file\n");
        return 0;
    }

AudioFileLoader() の中で data　はallocateされます。
よって、プログラムの最後でメモリを解放するのを忘れないでください。

読み込んだファイルのサンプル数から、FFTフレームがいくつ取れるかを計算します。

main.C
    //discard surplus.
    int frame_num = (int)(sfinfo.frames/FFTSIZE);

端数は捨てます。

書き出し用のオーディオデータを格納する為のメモリ領域を確保します。

main.C
    //for output
    float *outdata = calloc(frame_num*FFTSIZE, sizeof(float));

FFTのセットアップは以下のようにします。
fft_strという構造体を作りました。

main.C
    //set up for FFT
    fft_str fft;
    FFT_init(&fft,FFTSIZE);

fft_strは以下の通り。

main.C
typedef struct {
    FFTSetup fftsetup;
    vDSP_Length log2n;
    DSPSplitComplex splitComplex;

    int fftsize;
    int ffthalfsize;
}fft_str;

FFT分析用のメモリ領域を確保します。サイズはFFTSIZEです。

main.C
    //for FFT analysis
    float *buffer;
    buffer = malloc(FFTSIZE*sizeof(float));

そして、FFTの処理を、先ほど計算した frame_unm　回数、繰り返します。

main.C
for(i=0;i<frame_num*OVERLAP -1 ;i++){

        //copy data
        for(j=0;j<FFTSIZE;j++)
            buffer[j] = data[i*FFTSIZE/OVERLAP + j];


        FFT(buffer,&fft);
         /*
          ここに処理を書く
        */
        iFFT(buffer,&fft);

        //add data
        for(j=0;j<FFTSIZE;j++)
            outdata[i*FFTSIZE/OVERLAP + j] += buffer[j];

    }

FFT() でFFTを行い、 iFFT() にてiFFTを行います。
ここで、

main.C
//copy data
for(j=0;j<FFTSIZE;j++)
      buffer[j] = data[i*FFTSIZE/OVERLAP + j];

の部分では、オーディオデータの読み込み時にFFTSIZEの半分の長さだけシフトさせ、オーバラップを行っています。 OVERLAP = 2　です。オーバーラップについては、以下のサイトに詳しく説明がされています。

http://www.toyo.co.jp/page.jsp?id=2399

このプログラムでは、最初のフレームは0~2048サンプル、次が1024~3072サンプル、その次が2048~4096サンプルといった具合に、半フェーズずれてFFT、iFFTがなされます。

最後にiFFT処理が施されたオーディオデータは outdata に加算されます。(加算せずに単純に代入するだけだと、窓関数を使用している関係で入力時と同じ音データが形成されません。波形が窓関数形に波波になってしまいます。)

そして最後に outdata をwavデータに書き出します。

main.C
if(!AudioFileWriter(AUDIO_OUT,&sfinfo,outdata)){
        printf("Could not Write Audio file\n");
        return 0;
    }

メモリ解放
vDSP関連のメモリ領域は FFT_free() にて一括解放できるようにしました。

main.C
    FFT_free(&fft);
    free(data);
    free(buffer);

出力されたオーディオデータを聴いてみると、入力時のデータが再現されているのが分かります。

iFFT()

iFFT()　関数は以下の通りです。

main.C
void iFFT(float *wavedata, fft_str *fft)
{
    vDSP_fft_zrip(fft->fftsetup, &fft->splitComplex, 1, fft->log2n, FFT_INVERSE);

    //transform to real
    vDSP_ztoc( &fft->splitComplex, 1, ( COMPLEX * )wavedata, 2, fft->ffthalfsize );
}

単純に、 FFT_FORWARD 部分を FFT_INVERSE に変え、vDSP_ctoz の代わりに vDSP_ztoc を使用します。これで、 wavedata メモリ領域内にiFFTの処理結果が格納される事になります。

Filter

入力信号をフィルタリングして、特定の周波数成分を減衰させる為には、FFTによって得られた実数部と虚数部を強度(magnitude)と位相(Phase)部に変換し、強度の値を変化させる必要があります。Filteringのプログラムは、main関数の FFT() と iFFT() の間に書きます。

実数部と虚数部から強度と位相に変換する為には以下のようにします。

main.C
magnitude = sqrt(real*real + img*img);
phase = atan2(img, real);

また、強度と位相から実数部と虚数部に戻すやり方は以下の通りです。

main.C
real = magnitude * cos(phase);
img = magnitude * sin(phase);

よって、例えば始めの100個のFFTビンのmagnitudeを0にする場合は

main.C
float magnitude_array, phase_array;
for(j=0;j<100;j++){

        magnitude_array = sqrt(fft.splitComplex.realp[j]*fft.splitComplex.realp[j] +
                               fft.splitComplex.imagp[j]*fft.splitComplex.imagp[j]);
        phase_array = atan2(fft.splitComplex.imagp[j],fft.splitComplex.realp[j]);

        magnitude_array *= 0;

        fft.splitComplex.realp[j] = magnitude_array*cos(phase_array);
        fft.splitComplex.imagp[j] = magnitude_array*sin(phase_array);

}

となり、 magnitude_array に0 ~ 1.0までの値を掛け合わせる事によってフィルターの度合いを変化させる事ができます。

main.C
magnitude_array *= 0;

時間が絶つに従って、序々にスイープさせる効果をつけるなら、

main.C
magnitude_array *= (float)i/frame_num;

で、できます。

Scalingについて

前回も触れましたが、vDSPではFFT処理後に実部と虚部をスケーリングする必要がありました。FFT.cの以下の部分です。

main.C
    //scaling
    float scale = 1.0/(fft->fftsize);
    vDSP_vsmul(fft->splitComplex.realp, 1, &scale, fft->splitComplex.realp, 1, fft->ffthalfsize);
    vDSP_vsmul(fft->splitComplex.imagp, 1, &scale, fft->splitComplex.imagp, 1, fft->ffthalfsize);

試しにこれをコメントアウトしてiFFTを行うと、wavデータのレンジをオーバーし、データが適切に書き出されません。
オーバーラップを2以上にしている場合、OVERLAP分更に割る必要があります。

次回

次回はFFTによって得られたスペクトラルから得られる様々な音信号の情報を分析したいと思います。

34
28
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
34
28