サルでも分かるFFT（3）

ようやくFFTが分かった。結構長い道のりだった。

まずは、普通のO(n^2)のDFTを行列表現で書いてみます。このとき、eのべき乗のところを式で書くとイメージが大変なので、極座標系の矢印で書いてあげると分かりやすいです。

この式をじっくり見ると、「サルでも分かるFFT(2)」で紹介したk成分と(N-k)成分の対称性が見えます。Fourier visualizerがはき出した極座標系の周波数スペクトラムとの対比がよく分かると思います。

さて、FFTですが、まず偶数番と奇数番の周波数成分に分けます。偶数番成分に注目すると、

と書けますが、矢印の部分の4*8の行列は真ん中を中心に対象であることが分かります。よってこれは、

と書けます。この時点で演算数が半分になりました。これを繰り返せば分割統治の容量でO(N^2)からO(N log N)まで落とせそうです。

　同様に、奇数番のものは以下のように表すことができます。

これをゴリゴリやり続けると、最終的に以下の要素に辿り着きます。

ここまでくれば、あとは近い塊から順に例の「バタフライダイアグラム」(参照元：wikipedia、作成者: たまなるたみ氏)に当てはめればいいです。

左側のx[i]の並びですが、ビット逆転（ビット列に直して、reverseをする）が使われています。これは何故？と思いましたが、偶数と奇数の成分に逐次分けていったことを思い出せば理解できます。

まず、上半分は、偶数です。さらにその中の上半分の2個は、4で割れる数（下二桁のビットが00のもの）です。逆に下半分の2個は、ビットが10のものです。

これを再帰的に考えていけば、LSBとMSBの関係が逆転していることが分かります。それでビット逆転を使うと都合がいいわけです。

ほんで、C++のソース。最初pythonで書こうとしていましたが、C++のcomplexライブラリの方が使いやすかったので、C++で書きました。FFTとInverse FFTは、非常に似た処理を行うので、共通部分はgenFFT()にまとめています。ソースのコンパクトさよりも分かりやすさを追求したので、比較的読みやすい（はず？）です。

view plainprint?

1. typedef double Amp;               // amplitude at each sampling time  
2. typedef complex<double> Freq;     // constituent frequencies  
3.   
4. int bit_reverse(int x, int n) {  
5.     int rev = 0;  
6.   
7.     for (--n; n > 0; n >>= 1, x >>= 1)  
8.         rev = (rev << 1) | (x & 1);  
9.   
10.     return rev;  
11. }  
12.   
13. vector<Freq> genFFT(const vector<Freq> &x, int sgn) {  
14.     int N = x.size();  
15.     vector<Freq> X(N);  
16.   
17.     // set init values with bit reversed  
18.     for (int i = 0; i < N; i++) {  
19.         int rev = bit_reverse(i, N);  
20.         X[i] = x[rev];  
21.     }  
22.   
23.     // butterfly diagrams  
24.     for (int b = 2; b <= N; b <<= 1) {  
25.         for (int i = 0; i < N; i++) {  
26.             if (i % b >= b/2)  
27.                 continue;  
28.   
29.             int k = i + b/2;  
30.             Freq Wi = exp(Freq(0, sgn*2*PI*(i%b)/b));  
31.             Freq Wk = exp(Freq(0, sgn*2*PI*(k%b)/b));  
32.   
33.             X[i] += Wi*X[k];  
34.             X[k] = Wk*X[k] + (X[i] - Wi*X[k]);  
35.         }  
36.     }  
37.   
38.     return X;  
39. }  
40.   
41. vector<Freq> FFT(const vector<Amp> &x) {  
42.     int N = x.size();  
43.   
44.     // convert real numbers to complex numbers  
45.     vector<Freq> x_t(N);  
46.     for (int i = 0; i < N; i++)  
47.         x_t[i] = Freq(x[i], 0);  
48.   
49.     return genFFT(x_t, -1);  
50. }  
51.   
52. vector<Amp> IFFT(const vector<Freq> &X) {  
53.     int N = X.size();  
54.   
55.     vector<Freq> x_t = genFFT(X, 1);  
56.   
57.     // convert complex numbers to real numbers  
58.     // also normalize the amplituds (devide by N)  
59.     vector<Amp>x(N);  
60.     for (int i = 0; i < N; i++)  
61.         x[i] = x_t[i].real()/N;  
62.   
63.     return x;  
64. }  

typedef double Amp;               // amplitude at each sampling time
typedef complex<double> Freq;     // constituent frequencies

int bit_reverse(int x, int n) {
    int rev = 0;

    for (--n; n > 0; n >>= 1, x >>= 1)
        rev = (rev << 1) | (x & 1);

    return rev;
}

vector<Freq> genFFT(const vector<Freq> &x, int sgn) {
    int N = x.size();
    vector<Freq> X(N);

    // set init values with bit reversed
    for (int i = 0; i < N; i++) {
        int rev = bit_reverse(i, N);
        X[i] = x[rev];
    }

    // butterfly diagrams
    for (int b = 2; b <= N; b <<= 1) {
        for (int i = 0; i < N; i++) {
            if (i % b >= b/2)
                continue;

            int k = i + b/2;
            Freq Wi = exp(Freq(0, sgn*2*PI*(i%b)/b));
            Freq Wk = exp(Freq(0, sgn*2*PI*(k%b)/b));

            X[i] += Wi*X[k];
            X[k] = Wk*X[k] + (X[i] - Wi*X[k]);
        }
    }

    return X;
}

vector<Freq> FFT(const vector<Amp> &x) {
    int N = x.size();

    // convert real numbers to complex numbers
    vector<Freq> x_t(N);
    for (int i = 0; i < N; i++)
        x_t[i] = Freq(x[i], 0);

    return genFFT(x_t, -1);
}

vector<Amp> IFFT(const vector<Freq> &X) {
    int N = X.size();

    vector<Freq> x_t = genFFT(X, 1);

    // convert complex numbers to real numbers
    // also normalize the amplituds (devide by N)
    vector<Amp>x(N);
    for (int i = 0; i < N; i++)
        x[i] = x_t[i].real()/N;

    return x;
}

FFTが出来たので、次回「多倍長整数の除算をFFTで高速化する」でFFTシリーズをシメたいと思います。