離散フーリエ変換(1) - 素直な実装

故あってDFTについて書きます。

数学的な話は雑にやっていきます。　わからないので。

数学においてフーリエ変換（フーリエへんかん、英語: Fourier transform; FT）は実変数の複素または実数値函数を別の同種の函数に写す変換である。変換後の函数はもとの函数に含まれる周波数を記述し、しばしばもとの函数の周波数領域表現 (frequency domain representation) と呼ばれる。(Wikipedia)

$f(x)$ は $-\infty \le x +\infty$ で定義されている連続な関数とすると、その周波数 $\xi$ における成分 $F \left( \xi \right)$ は、

$F \left( \xi \right) = \int^{\infty}_{-\infty} f \left( x \right) e^{-2 \pi i x \xi} dx$

で得られます。

ただし、実用上は離散化/量子化された有限長のデータがよく使われます。

ということで離散フーリエ変換です。

離散フーリエ変換（りさんフーリエへんかん、英語: discrete Fourier transform、DFT）とは
離散群上のフーリエ変換であり、信号処理などで離散化されたデジタル信号の周波数解析などに
よく使われる。

長さ $n$ の時間領域のデータ列 $x_0, \dots , x_{n-1}$ を以下の処理で長さ $n$ の周波数領域のデータ $F_0, \dots, F_{n-1}$ に変換します。

$F_j = \sum^{n-1}_{k=0}x_k e^{-\frac{2 \pi i}{n}jk}$

ざっくりいうと $F \left( \frac{2 \pi}{n} j \right)$ が $F_j$ に対応していると考えて良いと思います。

これを素直に実装してみましょう。

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

using namespace std;

static const double PI = 3.14159265358979;

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> > x;
  cout << "Original Data" << endl;
  for(size_t i = 0; i < 10; ++i) {
    x.push_back(complex<double>(static_cast<double>(i), 0));
    cout << x.at(i) << endl;
  }
  cout << endl;
  
  // DFT
  vector< complex<double> > result = dft(x);

  cout << "DFT" << endl;
  for(size_t i = 0; i < result.size(); ++i) {
    cout << result.at(i) << endl;
  }
      
  // DFTの結果の複素共役をさらにDFTしてデータ数で割ると元データを取り出せる(逆変換)
  for(size_t i = 0; i <result.size(); ++i) {
    result.at(i) = conj(result.at(i));
  }

  cout << endl << "IDFT" << endl;
  vector< complex<double> > idft = dft(result);
  for(size_t i = 0; i < idft.size(); ++i) {
    cout << idft.at(i) / static_cast<double>(idft.size()) << endl;
  }
  return 0;
}

結果

Original Data
(0,0)
(1,0)
(2,0)
(3,0)
(4,0)
(5,0)
(6,0)
(7,0)
(8,0)
(9,0)

DFT
(45,0)
(-5,15.3884)
(-5,6.88191)
(-5,3.63271)
(-5,1.6246)
(-5,-1.52504e-13)
(-5,-1.6246)
(-5,-3.63271)
(-5,-6.88191)
(-5,-15.3884)

IDFT
(1.20259e-13,9.87654e-14)
(1,8.3844e-14)
(2,5.88862e-14)
(3,3.92575e-14)
(4,3.01981e-14)
(5,1.93623e-14)
(6,9.23706e-15)
(7,-1.75859e-14)
(8,-1.7053e-14)
(9,-3.71259e-14)

0,1, ..., 9の列がDFT、逆DFTを経て（大体）元に戻ったのが確認できると思います。

今回はここまでにします。