快速傅里葉變換

1、功能

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

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}}\)，將序列\(x(n)\)按序號\(n\)的奇偶分紅兩組，即
\[ \left.\begin{matrix}\begin{align*}x_{1}(n)=&x(2n)\\ x_{2}(n)=&x(2n+1)\end{align*}\end{matrix}\right\} \qquad n=0,1,...,\frac{N}{2}-1 \]
所以，\(x(n)\)的傅里葉變換可寫成
\[ \begin{align*}X(k) &= \sum_{n=0}^{N/2-1}x(2n)W^{2nk}_{N} + \sum_{n=0}^{N/2-1}x(2n+1)W^{(2n+1)k}_{N}\\&= \sum_{n=0}^{N/2-1}x_{1}(n)W^{nk}_{N/2} + W_{N}^{k}\sum_{n=0}^{N/2-1}x_{2}(n)W^{nk}_{N/2}\end{align*} \]
由此可得\(X(k)=X_{1}(k)+W_{N}^{k}X_{2}(k), \qquad k = 0,1,...,\frac{N}{2}\)，式中
\[ \begin{align*}X_{1}(k)&=\sum_{n=0}^{N/2-1}x(2n)W^{2nk}_{N}\\X_{2}(k)&=\sum_{n=0}^{N/2-1}x(2n+1)W^{(2n+1)k}_{N}\end{align*} \]
他們分別是\(x_1(n)\)和\(x_2(n)\)的\(N/2\)點DFT。上面的推導代表：一個\(N\)點DFT被分解爲兩個\(N/2\)點DFT，這兩個\(N/2\)點DFT又可合成一個\(N\)點DFT。但上面給出的公式僅能獲得\(X(k)\)的前\(N/2\)點的值，要用\(X_{1}(k)\)和\(X_{2}(k)\)來表達\(X(k)\)的後半部分的值，還必須運用權係數\(W_N\)的週期性與對稱性，即
\[ W_{N/2}^{n(k+N/2)}=W_{N/2}^{nk}, \quad W_{N}^{(k+N/2)}=-W_{N}^{k} \]
所以，\(X(k)\)的後\(N/2\)點的值可表示爲
\[ \begin{align*}X(k+\frac{N}{2})&=X_{1}(k+\frac{N}{2})+W_{N}^{k+N/2}X_{2}(k+\frac{N}{2})\\&=X_{1}(k)-W_{N}^{k}X_{2}(k), \ k=0,1,...,\frac{N}{2}-1\end{align*} \]
經過上面的推導能夠看出，\(N\)點的DFT能夠分解爲兩個\(N/2\)點DFT，每一個\(N/2\)點DFT又能夠分解爲兩個\(N/4\)點DFT。依此類推，當\(N\)爲2的整數次冪時（\(N=2^M\)），因爲每分解一次下降一階冪次，因此經過\(M\)次分解，最後所有成爲一系列2點DFT運算。以上就是按時間抽取的快速傅里葉變換（FFT）算法。

序列\(X(k)\)的離散傅里葉反變換定義爲
\[ x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_{N}^{-nk}, \qquad n=0,1,...,N-1 \]
它與離散傅里葉正變換的區別在於將\(W_N\)改變爲\(W_N^{-1}\)，並多了一個除以\(N\)的運算。由於\(W_N\)和\(W_N^{-1}\)對於推導按時間抽取的快速傅里葉變換算法並沒有實質性區別，所以可將FFT和快速傅里葉反變換（IFFT）算法合併在同一程序中。

3、使用說明

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

/************************************
    x       ---一維數組，長度爲n，開始時存放要變換數據的實部，最後存放變換結果的實部。
    y       ---一維數組，長度爲n，開始時存放要變換數據的虛部，最後存放變換結果的虛部。
    n       ---數據長度，必須是2的整數次冪。
    sign    ---當sign=1時，子函數計算離散傅里葉正變換；當sign=-1時，子函數計算離散傅里葉反變換
************************************/
#include "math.h"

void fft(double *x, double *y, int n, int sign)
{
    int i, j, k, l, m, n1, n2;
    double c, c1, e, s, s1, t, tr, ti;
    for(j = 1, i=1; i < 16; i++) {
        m = i;
        j = 2 * j;
        if(j == n)
            break;
    }
    n1 = n - 1;
    for(j = 0, i = 0; i < n1; i++) {
        if(i < j) {
            tr = x[j];
            ti = j[j];
            x[j] = x[i];
            y[j] = j[i];
            x[i] = tr;
            y[i] = ti;
        }
        k = n / 2;
        while(k < (j + 1)) {
            j = j - k;
            k = k / 2;
        }
        j = j + k;
    }
    n1 = 1;
    for(l = 1; l <= m; l++) {
        n1 = 2 * n1;
        n2 = n1 / 2;
        e = 3.14159265359 / n2;
        c = 1.0;
        s = 0.0;
        c1 = cos(e);
        s1 = -sign * sin(e);
        for(j = 0; j < n2; j++) {
            for(i = j; i < n; i += n1) {
                k = i + n2;
                tr = c * x[k] - s * y(k);
                ti = c * y[k] + s * x[k];
                x[k] = x[i] - tr;
                y[k] = y[i] - ti;
                x[i] = x[i] + tr;
                y[i] = y[i] + ti;
            }
            t = c;
            c = c * c1 - s * s1;
            s = t * s1 + s * c1;
        }
    }
    if(sign == -1) {
        for(i = 0; i < n; i++) {
            x[i] /= n;
            y[i] /= n;
        }
    }

}