基4快速傅里葉變換

1、功能

計算復序列的基4快速傅里葉變換。

2、方法簡介

序列\(x(n)(n=0,1,...,N-1)\)的離散傅里葉變換定義爲
\[ X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}, \qquad k=0,1,...,N-11 \]
其中\(W_{N}^{nk}=e^{-j\frac{2\pi nk}{N}}\)，若\(N=4^M\)，則將序列\(x(n)\)分紅四個\(N/4\)點的序列\(x_1(n) 、x_2(n) 、x_3(n) 、x_4(i)(n=0,1, …,N/4—1)\), 即
\[ x(n) =x_1(n) +x_2(n) +x_3(n) +x_4(n) \]
式中
\[ \left\{\begin{matrix}\begin{align*}x_{1}(n)&=x(n), & (n=0,1,...,N /4 - 1)\\ x_{2}(n)&=x(n+\frac{N}{4}), & (n=0,1,...,N /4 - 1)\\ x_{3}(n)&=x(n+\frac{N}{2}), & (n=0,1,...,N /4 - 1)\\ x_{4}(n)&=x(n+\frac{3N}{4}), & (n=0,1,...,N /4 - 1)\end{align*}\end{matrix}\right. \]
把\(x(n)\)代入DFT表達式中，則有
\[ \begin{align*}X(k)&=\sum_{n=0}^{N/4-1}[x_{1}(n)W_{N}^{nk}+x_{2}(n)W_{N}^{(n + N/4)k}+x_{3}(n)W_{N}^{(n + N/2)k}+x_{4}(n)W_{N}^{(n + 3N/4)k}]\\&=\sum_{n=0}^{N/4-1}[x_{1}(n)+x_{2}(n)W_{N}^{N/4k}+x_{3}(n)W_{N}^{N/2k}+x_{4}(n)W_{N}^{3N/4k}]W_{N}^{nk}\\\end{align*}\\(k=0,1,...,N-1) \]
把\(X(k)\)按頻率抽取，得
\[ \left\{\begin{matrix}\begin{align*}X(4k)&=\sum_{n=0}^{N/4-1}[x_{1}(n)+x_{2}(n)+x_{3}(n)+x_{4}(n)]W_{N/4}^{nk}\\X(4k+1)&=\sum_{n=0}^{N/4-1}[x_{1}(n)-jx_{2}(n)-x_{3}(n)+jx_{4}(n)]W_{N}^{n}W_{N/4}^{nk}\\X(4k+2)&=\sum_{n=0}^{N/4-1}[x_{1}(n)-x_{2}(n)+x_{3}(n)-x_{4}(n)]W_{N}^{2n}W_{N/4}^{nk}\\X(4k+3)&=\sum_{n=0}^{N/4-1}[x_{1}(n)+jx_{2}(n)-x_{3}(n)-jx_{4}(n)]W_{N}^{3n}W_{N/4}^{nk}\end{align*}\end{matrix}\right. \]
經過上面得分解，可求得全部\(X(k)\)值，其基本運算式爲
\[ \left\{\begin{matrix}\begin{align*}f_{1}(n)&=x_{1}(n)+x_{2}(n)+x_{3}(n)+x_{4}(n)\\f_{2}(n)&=[x_{1}(n)-jx_{2}(n)-x_{3}(n)+jx_{4}(n)]W_{N}^{n}\\f_{3}(n)&=[x_{1}(n)-x_{2}(n)+x_{3}(n)-x_{4}(n)]W_{N}^{2n}\\f_{4}(n)&=[x_{1}(n)+jx_{2}(n)-x_{3}(n)-jx_{4}(n)]W_{N}^{3n}\end{align*}\end{matrix}\right.，(k=0,1,...,\frac{N}{4}-1) \]
這樣，就將一個\(N\)點得DFT轉化爲四個\(N/4\)點DFT來計算。依此類推，直至分解到最後一級。以上就是按頻率抽取的基4快速傅里葉變換算法。與基2FFT相比，基4FFT的乘法量約減小25%，加法量也略有減小。

3、使用說明

是用C語言實現基4快速傅里葉變換（FFT）的方法以下：

/************************************
    x       ---一維數組，長度爲n，開始時存放要變換數據的實部，最後存放變換結果的實部。
    y       ---一維數組，長度爲n，開始時存放要變換數據的虛部，最後存放變換結果的虛部。
    n       ---數據長度，必須是4的整數次冪。
************************************/
#include "math.h"

void fft(double *x, double *y, int n)
{
    int i, j, k, m, il, i2, i3, nl, n2;
    double a, b ,c ,e ,rl ,r2 ,r3 ,r4 ,s1 ,s2 ,s3 ,s4;
    double col, co2, co3, sil, si2, si3;

    for(j = 1; i = 1; i < 10; i++) {
        m = i;
        j = 4 * j;
        if(j == n) break;
    }
    n2 = n;
    for(k = 1; k <= m; k++) {
        n1 = n2;
        n2 = n2 / 4;
        e = 6.28318530718 / nl;
        a = 0;
        for(j = 0; j < n2; j++) {
            b = a + a;
            c = a + b;
            co1 = cos(a);
            co2 = cos(b);
            co3 = cos(c);
            sil = sin(a);
            si2 = sin(b);
            si3 = sin(c);
            a = (j + l) * e;
            for(i = j; i < n; i = i +n1) {
                il = i + n2;
                i2=il+n2;
                i3 = i2 + n2;
                nl = n-1;
                r1 = x[i] + x[i2];
                r3 = x[i] - x[i2];
                s1 = y[i] + y[i2];
                s3 = y[i] - y[i2];
                r2 = x[il] + x[i3];
                r4 = x[il] — x[i3];
                s2 = y[il] + y[i3];
                s4 = y[il] - y[i3];
                x[i] = rl — r2;
                r2 = r1 — r2;
                rl = r3 — s4;
                r3 = r3 + s4;
                y[i] = s1 + s2;
                s2 = s1 - s2;
                s1 = s3 + r4;
                s3 = s3 — r4;
                x[il] = col * r3 + sil * s3;
                y[il] = col* s3 — sil * r3;
                x[i2] = co2 * r2 + si2 * s2;
                y[i2] = co2 * s2 — si2 * r2;
                x[i3] = co3 * rl + si3 * s1;
                y[i3] = co3 * s1 — si3 * r1;
            }
        }
    }
    n1 = n - 1;
    for(j = 0, i = 0; i < n1; i++) {
        if(i < j) {
            rl = x[j］；
            s1 = y[j];
            x[j] = x[i];
            y[j] = y[i];
            x[i] = rl;
            y[i] = sl;
        }
        k = n / 4;
        while( 3 * k < (j + 1)) {
            j = j - 3 * k;
            k = k / 4;
        }
        j = j + k;
    }
}