2013/02/07

FFTによるconvolution(畳み込み)実験

「DFTによる畳み込み」記事と同じ内容をFFTでもやってみた。 内容は下図の通りで、1,1,0,0 の4つの入力信号と 0,1,0,0 のフィルタを、それぞれFFT(高速フーリエ変換)して周波数領域にしてしまう。そして周波数領域で複素数の乗算をして、iFFT(逆高速フーリエ変換)で戻すというもの。フィルタの内容は1サンプル遅れるディレイなので、結果は 0,1,1,0 となり、時間軸での畳み込みと同じ結果が出る。
FFTで重要なのはサンプル数が2^nであるということ。今回は4という少ないサンプル数で試している。 この記事の例はサンプル数が少なすぎて、現実感がない。実際の入力信号は、かなり長いのが普通だし、リアルタイムであれば、長さが不明ということになる。フィルタは普通数千サンプルぐらいで、長くても数万サンプル。この長さの違うサンプルを畳み込む必要がある。その方法は overlap-add method という畳み込み手法で、次回の記事で書いてみたいと思う。この方法を使えば、高速なFIRデジタルフィルタを実現できる。実際に多くのエフェクトで利用されている。Audacity の Equlaization や Convolution Reverb などがよい例。

FFT convolution 実験用ソースコード

/* fft_test.c */
#include <stdio.h>
#include <math.h>

/* FFT Cooley-Tukey型FFTアルゴリズム */
void fft(int *n, double *t, double *ar, double *ai) {
    int m, mh, j, k, irev, i=0;
    double wr, wi, xr, xi;
    for (j=1; j<*n-1; j++){
      for(k=*n >>1; k>(i^=k); k>>=1){}
        if(j < i){
            xr = *(ar+j);
            xi = *(ai+j);
            *(ar+j) = *(ar+i);
            *(ai+j) = *(ai+i);
            *(ar+i) = xr;
            *(ai+i) = xi;
        }
    }
    for(mh=1; (m=mh<<1)<=*n; mh=m){
        irev = 0;
        for(i=0; i<*n; i+=m){
            wr = cos(*t * irev);
            wi = sin(*t * irev);
            for(k=*n>>2; k>(irev^=k); k>>=1){}
            for(j=i; j<mh+i; j++){
                k = j + mh;
                xr = *(ar+j) - *(ar+k);
                xi = *(ai+j) - *(ai+k);
                *(ar+j) += *(ar+k);
                *(ai+j) += *(ai+k);
                *(ar+k) = wr * xr - wi * xi;
                *(ai+k) = wr * xi + wi * xr;
            }
        }
    }
}

int main(void) {
  int i,j;
  /* 入力信号 実部 */
  double sr[] ={1,1,0,0};
  int N = sizeof(sr)/8;
  /* 入力信号 虚部 */
  double si[N];
  for(i=0; i<sizeof(sr)/8; i++){
     si[i]=0;
  }

  /* 畳み込む関数 実部 */
  double sr2[] ={0,1,0,0};
  double si2[N];
  /* 畳み込む関数 虚部 */
  for(i=0; i<sizeof(sr)/8; i++){
     si2[i]=0;
  }

  /* フーリエ変換後の数値を入れる配列を用意 */
  double p[N];/* 振幅スペクトル */
  double p2[N];/* 振幅スペクトル */
  double hr3[N];/* 実部 */
  double hi3[N];/* 虚部 */
  double t = 2*M_PI/N;/* シータ */

 /* FFT 入力信号 */
 fft(&N,&t,&sr[0],&si[0]);
 printf("\nN:\tFFT1 Real\tFFT1 Imaginary\tFFT1 Power\n");
 for(i=0; i<N; i++){
     si[i] *= -1;
     /* 振幅スペクトル */
     p[i] = sqrt(sr[i]*sr[i] + si[i]*si[i]);
     /* 実部 虚部 振幅スペクトル 表示 */
     printf("%d:\t%12.11f\t%12.11f\t%12.11f\n"
            ,i,sr[i],si[i],p[i]);
 }

 /* FFT フィルタ */
 fft(&N,&t,&sr2[0],&si2[0]);
 printf("\nN:\tFFT2 Real\tFFT2 Imaginary\tFFT2 Power\n");
 for(i=0; i<N; i++){
     si2[i] *= -1;
     /* 振幅スペクトル */
     p2[i] = sqrt(sr2[i]*sr2[i] + si2[i]*si2[i]);
     /* 実部 虚部 振幅スペクトル 表示 */
     printf("%d:\t%12.11f\t%12.11f\t%12.11f\n"
           ,i,sr2[i],si2[i],p2[i]);
 }

 /* 畳み込み */
 printf("\n");
 for(i=0; i<N; i++){
   if(si[i]==0){
     if(sr[i]==0){
       hr3[i] = 0;
       hi3[i] = 0;
     }else{
       hr3[i] = sr[i]*sr2[i];
       hi3[i] = si2[i];
     }
   }else{
     hr3[i] = sr[i]*sr2[i]+(si[i]*si2[i])*-1;
     hi3[i] = sr[i]*si2[i]+si[i]*sr2[i];
   }
 } 
 printf("N:\tconvolution R\tconvolution i\n");
 for(i=0; i<N; i++){
     printf("%d:\t%12.11f\t%12.11f\n",i,hr3[i],hi3[i]);
 }

 /* iFFT  FFTと同じものを使う */
 fft(&N,&t,&hr3[0],&hi3[0]);
 printf("\nN:\tiFFT Real\tiFFT Imaginary\n");
 for(i=0; i<N; i++){
     hi3[i] *= -1;
     hi3[i] /= N;
     hr3[i] /= N;
     printf("%d:\t%12.11f\t%12.11f\n",i,hr3[i],hi3[i]);
 }

 return 0;
}
コンパイルして実行すると以下のような結果が表示される。1サンプル遅れているのが確認できる。
N: FFT1 Real       FFT1 Imaginary  FFT1 Power
0: 2.00000000000   0.00000000000   2.00000000000
1: 1.00000000000   -1.00000000000  1.41421356237
2: 0.00000000000   0.00000000000   0.00000000000
3: 1.00000000000   1.00000000000   1.41421356237

N: FFT2 Real       FFT2 Imaginary  FFT2 Power
0: 1.00000000000   0.00000000000   1.00000000000
1: 0.00000000000   -1.00000000000  1.00000000000
2: -1.00000000000  0.00000000000   1.00000000000
3: -0.00000000000  1.00000000000   1.00000000000

N: convolution R   convolution i
0: 2.00000000000   0.00000000000
1: -1.00000000000  -1.00000000000
2: 0.00000000000   0.00000000000
3: -1.00000000000  1.00000000000

N: iFFT Real       iFFT Imaginary
0: 0.00000000000   0.00000000000
1: 1.00000000000   0.00000000000
2: 1.00000000000   -0.00000000000
3: 0.00000000000   -0.00000000000


sound programming 目次はこちら