プログラムの国のフーリエ変換

やあ子どもたち。花粉症は大丈夫か。今年は花粉症どころかPM2.5の方が空気中の微粒子として何やら気になるから無用の外出は禁物だぞ。
今日はプログラムの世界でのフーリエ変換の話をしよう。コンピュータを使ったプログラムで出来るのは離散フーリエ変換(DFT: DiscreteFourierTransform)だけだよ。だからプログラムの世界に住んでいる俺達にとって、フーリエ変換といったら離散フーリエ変換のことなんだから気をつけてくれよ。
さてフーリエ変換ってなんだ。まずは入力が実関数のみと限定した場合（たいていの場合はそうだ）、プログラムにすると、こんなだ。

void dft(
         const vector< double >& a,// [入力] 実数数値列としての信号
         vector< double >* Re,// [出力] 実数部
         vector< double >* Im// [出力] 虚数部
         )
{
    // dft
    const auto& N = a.size();
    for( int j=0; j<N; ++j )
    {
        double Re_sum = 0;
        double Im_sum = 0;
        for( int i=0; i<N; ++i )
        {
            double tht = 2*M_PI/N * j * i;
            Re_sum += a[i] *cos( tht );
            Im_sum += a[i] *sin( tht );
        }// i
        Re->push_back( Re_sum );
        Im->push_back( Im_sum );
    }// j
    return;
}

まずわかるのは、「実数値の離散的な列、vector< double > をフーリエ変換したら、同じ配列長さの、vector< double >が2本帰ってくる」ということ。2本というのは、sin成分と、cos成分の2本。まずは簡単のため、sin成分についてだけ考えていこう。
入力信号を処理していくイメージを図にしてみたよ。冒頭のソースコードを見てまず思い浮かぶのはこんな図だろう。

このノイズの入ったノコギリ波の繰り返し信号が、今回のフーリエ変換の入力信号 vector< double > a だよ。これを、配列の添字(i)を横軸に、（この例では300要素）、そしてその横軸の添字値に対応するaの要素値（a[i]）を縦軸に見立てて描画したものと思ってほしい。
赤い信号は、この同じ大きさの配列空間の中に定義される、1番波長の長いsinカーブ、そして2番めに波長の長いsinカーブの様子だ。それぞれ、波長がちょうど配列全体と同じ大きさ、波長がちょうど配列全体の1/2の大きさという具合だ。
次に、こうして考えたsinカーブと入力信号aとの相関値を計算していくのだ。これは上記コードの中で、Re_sum,Im_sumに値を足しこんで行ってる部分に相当する処理だよ。相関値ってなんだ。相関値計算のイメージはこんな感じ。

相関というのはこのグラフで言うと、入力信号とsinカーブどうしで、同じX軸値に対応するY軸値どうしを掛けあわせて作った新しいグラフを作り、全てのY軸値を足しあわせたものだよ。図の緑色の部分の面積（負の方向に突き出た面積は負の面積値として考える）が、相関値そのものとなる。上の図は1番目と2番目のsinカーブと、今回のノコギリ波入力信号との相関値計算のイメージを表したものだよ。それぞれのサインカーブの場合について、１つずつ、相関値が計算されていくわけです。
ちょうど物理の時間に習った、弦の振動もしくは閉管の中の音波の状況に似ているけど、まさに同じ感じでこのように波形を考えるわけなんだね。
さてこの調子で、N番目のsinカーブ（=波長が配列全体の1/Nの大きさをもつ波）に至るまでsinカーブを作成しては、入力信号との相関値を計算するということをしていく。参考までに、以下にそれぞれ、23番目と230番目のサインカーブの場合の図も載せておくぞ。

でまあ一応、cosカーブの方のイメージ図も乗せておくぞ。cosについても上記のsinの時と全く同じように、相関値を計算していくのだ。

こうして順次計算した相関値を配列にまとめたものが、そのまま、入力信号aに対するフーリエ変換の出力結果になる。上述のプログラムで言うとcosカーブの分はReという変数に、sinカーブの分はImという出力変数に入る。
Reと、Imの横軸はなんだろうと考える時、それはまさに、今見てきた、sinやcosカーブの波長の長さそのものに対応していることがわかるはずだ。そのものずばりで言えば周波数だ。単位は、整数K。その意味は、「入力信号の全横軸定義域をK分割した波長に対応する周波数」だ。どうだい実にすっきりした話じゃないかこれ。
フーリエ変換とはどんな入力情報に対し何を計算しどんな出力をしているのかという説明は以上の限りだ。
あと、配列長さ20の場合に同じ事をやった図も載せておこう。こっちの方が、あくまで離散的な数値計算をやってるだけというイメージが湧くかも知れない。とにかく、配列長さ300だろうが、20だろうが、話としては全く変わらないのだということを理解してほしい。

f:id:nurs:20130311003123p:image:w180:h120

f:id:nurs:20130311003227p:image:w180:h120

ところでみんなは、std::complexを知っているかなそれは過去日記でも触れているのだけど、こいつを使うと上述のプログラムがさらにこんなに短く書けてしまうのだ。

void dft( const vector< complex<double> >& a,
                   vector< complex< double > >* pout )
{
    const auto& N = a.size();
    pout->reserve( N );
    for( int j=0; j<N; ++j )
    {
        complex< double > sum(0,0);
        for( int i=0; i<N; ++i )
        {
            double tht = 2*M_PI/N * j * i;
            sum += a[i] * exp( complex< double >( 0, -tht ) );
        }// i
        pout->push_back( sum );
    }// j
    return;
}

またこう書くと、入力を実関数に限定しないより一般的なツールにもなってくるという利点があるよ。やー、コードが簡潔であること、コードが小さいって実にいいことだね。

こうなってくると、「複素数値の離散的な列、vector< complex > をフーリエ変換したら、同じ配列長さの、vector< complex >が1本、帰ってくる」ということが言えそうだね。

（おぼえておこう！）
複素数値の離散的な列、vector< complex<double> > をフーリエ変換したら、同じ配列長さの、vector< complex<double> >が1本、帰ってくる

さて以上で第一部は終了だよ。続いて第二部。第二部はフーリエ変換、DFTの計算を超高速に計算するFFTのというアルゴリズムの紹介だよ。