離散フーリエ変換(3) - 畳み込み定理

Wikipedia 畳み込み
http://ja.wikipedia.org/wiki/%E7%95%B3%E3%81%BF%E8%BE%BC%E3%81%BF

にあるように $F$ を $f_i$ をフーリエ変換したもの、 $G$ を $g_i$ をフーリエ変換したものとすると $f_i$ と $g_i$ の畳み込みをフーリエ変換した結果は、 $F_i \times G_i$ となります。

前回の式(*)で、各 $1 \le \tau < n$ に対して $C_\tau$ を求めるのに $n$ 回の乗算が必要なのですべての $C_\tau$ を求めるためには全部で $n^2$ 回の乗算が必要です。

これと同じ結果が、 $f_i$ のフーリエ変換、 $g_i$ のフーリエ変換、 $F_i$ と $G_i$ の乗算、最後にフーリエ逆変換で求めることができるよということです。

ここで、フーリエ変換と逆変換は高速フーリエ変換で高速に計測することが出来るので $n^2$ 回の乗算を素直にやるより高速に計算することができますよということです。

今回はとりあえず普通のDFTでやります。

#include <iostream>
#include <vector>
#include <complex>

using namespace std;

static const double PI = 3.14159265358979;

// (1)でやったやつ
vector< complex<double> > dft(const vector< complex<double> > &x) {
  vector< complex<double> > result;
  for(size_t j = 0; j < x.size(); ++j) {
    complex<double> F_j(0.0, 0.0);
    for(size_t k = 0; k < x.size(); ++k) {
      double angle = -2.0 * PI * j * k / x.size();
      F_j += x.at(k) * complex<double>(cos(angle), sin(angle));
    }
    result.push_back(F_j);
  }

  return result;
}

int main() {
  vector< complex<double> > orig;
  for(int i = 0; i < 30; ++i) {
    if(i >= 10 && i < 20)
      orig.push_back(complex<double>(i % 10, 0));
    else
      orig.push_back(complex<double>(2.0, 0));
  }

  vector< complex<double> > tmpl;
  for(int i = 0; i < 30; ++i) {
    if(i < 10)
      tmpl.push_back(complex<double>(i % 10, 0));
    else
      tmpl.push_back(complex<double>(0, 0));
  }

  vector< complex<double> > orig_f = dft(orig);
  vector< complex<double> > tmpl_f = dft(tmpl);

  vector< complex<double> > orig_tmpl_f;
  for(int i = 0; i < 30; ++i) {
    orig_tmpl_f.push_back(conj(orig_f.at(i)) * tmpl_f.at(i));
  }

  vector< complex<double> > cross_corr = dft(orig_tmpl_f);
  for(int i = 0; i < 30; ++i) {
    cout << cross_corr.at(i) / 30.0 << endl;
  }

  return 0;
}

この結果が以下です

(90,3.56598e-12)
(72,3.16807e-12)
(65,2.82322e-12)
(68,2.47079e-12)
(80,2.01605e-12)
(100,1.41351e-12)
(127,1.03834e-12)
(160,9.4739e-13)
(198,7.42754e-13)
(240,1.45898e-13)
(285,-5.36223e-13)
(258,-1.06629e-12)
(230,-1.26382e-12)
(202,-1.25056e-12)
(175,-1.6712e-12)
(150,-2.20552e-12)
(128,-2.69817e-12)
(110,-3.03923e-12)
(97,-3.47124e-12)
(90,-4.01693e-12)
(90,-4.54368e-12)
(90,-3.0847e-12)
(90,-1.84552e-12)
(90,-2.42532e-13)
(90,6.82121e-14)
(90,-2.12215e-13)
(90,7.48912e-13)
(90,1.19182e-12)
(90,3.8464e-13)
(90,6.36646e-13)

前回の結果が

なのでちゃんと一致していることが分かります。

次回は2次元のDFTをやっていきたいとおもいます。