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

説明用としてDFTを使った畳み込みプログラムを作ってみた。本来はDFTを高速にしたアルゴリズムであるFFTで作るものだが、流れを確認するために簡単なDFTにしてみた。

式とかはこちらを参考にした。
http://ja.wikipedia.org/wiki/畳み込み

下のプログラムでは、入力サンプルは4つで{1,1,0,0}というもの。フィルタも4つで{0,1,0,0}というもの。これは1サンプル遅れて出力されることを意味する。つまりディレイ。出力結果は{0,1,1,0}となる。この例は単純であるが、実際には周波数成分のコントロールや残響のコントロールなど幅広い応用が考えられる。

これを実現する方法で一番シンプルな方法は時間領域で畳み込むこと。しかし、データ量が増えると処理速度の問題が出てしまう。そこで、下図のように一度DFTで周波数領域に置き換えて、そこで演算を行い、その結果をiDFTで戻して出力をするという一見面倒な処理をするのだが、FFTのおかげで時間領域の畳み込みよりも高速に処理することができる。ここでは学習向けにDFTを使っているが、本来はFFT。FFTによる畳み込みはFIRの各種デジタルフィルター、EQ、HPF、Convolution Reverbなどで広く利用されている。

同じ内容をFFTでもやってみた。

DFT convolution 実験用ソースコード

このプログラムは数式をそのままにして自己流で作ったもの。参考にしているプログラムはないです。ということで間違っている可能性あり。また上から順に説明するために処理を関数にしてないです。複素数計算のところは怪しい。本来はもっとスマートに書くのだろう。またcomplex.hの利用も考えたけど、直接いじっている感じがしなかったのでやめた。またやたらと配列を使っているが、これも実用的なプログラムの場合はポインタを使う。まぁいろいろ突っ込みどころの多そうなサンプルです。そのかわり、とても素直で理解しやすいと思う。
/* dftcon.c */
#include <stdio.h>
#include <math.h>
/*************************************************************/
int main(void){
  int i,j;
/*************************************************************/
/* 入力信号 実部 */
  double sr[] ={
    1,
    1,
    0,
    0
    };
  int N = sizeof(sr)/8;
  /* 入力信号 虚部 すべて0を入力 */
  double si[N];
  for(i=0; i<sizeof(sr)/8; i++){
     si[i]=0;
  };
/*************************************************************/
/* フィルタ 実部 */
  double sr2[] ={
     0,
     1,
     0,
     0
     };
  double si2[N];
  /* フィルタ 虚部 すべて0を入力 */
  for(i=0; i<sizeof(sr)/8; i++){
     si2[i]=0;
  };
/*************************************************************/
/* 変換後の数値を入れる配列を用意 */
  double hr[N];/* 実部 */
  double hi[N];/* 虚部 */
  double p[N];/* 振幅スペクトル */

  double hr2[N];/* 実部 */
  double hi2[N];/* 虚部 */
  double p2[N];/* 振幅スペクトル */

  double hr3[N];/* 実部 */
  double hi3[N];/* 虚部 */
/*************************************************************/
/* DFT処理 入力信号 フーリエ変換 */
 double t = 2*M_PI/N;
 for(i=0; i<N; i++){ 
   hr[i] = 0;
   hi[i] = 0;
   for(j=0; j<N; j++){ 
     hr[i] += (sr[j]*cos(t*i*j))+(si[j]*sin(t*i*j));/* 実部 */
     hi[i] += (si[j]*cos(t*i*j))-(sr[j]*sin(t*i*j));/* 虚部 */
   }
 }
 /* 振幅を求めつつ結果を表示 */
 printf("\nN:\tDFT1 Real\tDFT1 Imaginary\tDFT1 Power\n");
 for(i=0; i<N; i++){
     /* 振幅スペクトル */
     p[i] = sqrt(hr[i]*hr[i] + hi[i]*hi[i]);
     /* 実部 虚部 振幅スペクトル 表示 */
     printf("%d:\t%12.11f\t%12.11f\t%12.11f\n",i,hr[i],hi[i],p[i]);
 }
/*************************************************************/
/* DFT処理 フィルタ フーリエ変換 */
 for(i=0; i<N; i++){ 
   hr2[i] = 0;
   hi2[i] = 0;
   for(j=0; j<N; j++){ 
     hr2[i] += (sr2[j]*cos(t*i*j))+(si2[j]*sin(t*i*j));/* 実部 */
     hi2[i] += (si2[j]*cos(t*i*j))-(sr2[j]*sin(t*i*j));/* 虚部 */
   }
 }
 /* 振幅を求めつつ結果を表示 */
 printf("\nN:\tDFT2 Real\tDFT2 Imaginary\tDFT2 Power\n");
 for(i=0; i<N; i++){
     /* 振幅スペクトル */
     p2[i] = sqrt(hr2[i]*hr2[i] + hi2[i]*hi2[i]);
     /* 実部 虚部 振幅スペクトル 表示 */
     printf("%d:\t%12.11f\t%12.11f\t%12.11f\n",i,hr2[i],hi2[i],p2[i]);
 }
/*************************************************************/
/* 周波数領域で演算(畳み込み) 複素数計算*/
 printf("\n");
 for(i=0; i<N; i++){
   if(hi[i]==0){
      if(hr[i]==0){
         hr3[i] = 0;
         hi3[i] = 0;
      }else{
         hr3[i] = hr[i]*hr2[i];
         hi3[i] = hi2[i];          
      }
   }else{
      hr3[i] = hr[i]*hr2[i]+(hi[i]*hi2[i])*-1;
      hi3[i] = hr[i]*hi2[i]+hi[i]*hr2[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]);
 }
/*************************************************************/
/* iDFT処理 逆フーリエ変換 */
 for(i=0; i<N; i++){
   sr[i] = 0;/*初期化*/
   si[i] = 0;/*初期化*/
   for(j=0; j<N; j++){ 
     sr[i] += (hr3[j]*cos(t*i*j))-(hi3[j]*sin(t*i*j));/* 実部 */
     si[i] += (hi3[j]*cos(t*i*j))+(hr3[j]*sin(t*i*j));/* 虚部 */
   }
   sr[i] /= N;
   si[i] /= N;
 }
 printf("\nN:\tiDFT Real\tiDFT Imaginary\n");
 /* 畳み込んだ結果を表示 */
 for(i=0; i<N; i++){
     printf("%d:\t%12.11f\t%12.11f\n",i,sr[i],si[i]);
 }
/*************************************************************/
 return 0;
}
コンパイルして実行すると以下のような結果が表示される。1サンプル遅れているのが確認できる。
N:      DFT1 Real       DFT1 Imaginary  DFT1 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:      DFT2 Real       DFT2 Imaginary  DFT2 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:      iDFT Real       iDFT Imaginary
0:      -0.00000000000  -0.00000000000
1:      1.00000000000   -0.00000000000
2:      1.00000000000   0.00000000000
3:      0.00000000000   0.00000000000


sound programming 目次はこちら