實序列快速傅里葉變換（二）

1、功能

用\(N\)點復序列快速傅立葉變換來計算\(2N\)點實序列的離散傅立葉變換。

2、方法簡介

假設\(x(n)\)是長度爲\(2N\)的實序列，其離散傅立葉變換爲
\[ X(k)=\sum_{n=0}^{2N-1}x(n)W_{2N}^{nk} \ , \ k=0,1,...,2N-1 \]
爲有效地計算傅立葉變換\(X(k)\), 咱們將\(x(n)\)分爲偶數組和奇數組，造成兩個新序列\(x(n)\)和\(g(n)\)，即
\[ \left\{\begin{matrix}\begin{align*}f(n)&=x(2n)\\ g(n)&=x(2n+1)\end{align*}\end{matrix}\right. , n=0,1,...,N-1 \]
而後將\(f(n)\)和\(g(n)\)組成一個復序列\(h(n)\)
\[ h(n)=f(n)+jg(n), \ n = 0,1,...,N-1 \]
用FFT計算\(h(n)\)的\(N\)點傅立葉變換\(H(k)\), 而且\(H(k)\)可表示爲
\[ H(k)=F(k)+jG(k), \ n = 0,1,...,N-1 \]
由上容易推出
\[ \left\{\begin{matrix}\begin{align*}F(k)&=\frac{1}{2}[H(k)+H^{*}(N-k)]\\ G(k)&=-\frac{j}{2}[H(k)-H^{*}(N-k)]\end{align*}\end{matrix}\right. , n=0,1,...,N-1 \]
求得\(F(k)\)和\(G(k)\)後，利用下面的蝶形運算計算\(x(n)\)的離散傅立葉變換\(X(k)\)
\[ \left\{\begin{matrix}\begin{align*}X(k)&=F(k)+G(k)W_{2N}^{k}\\ X(k+n)&=F(k)-G(k)W_{2N}^{k}\end{align*}\end{matrix}\right. , n=0,1,...,N-1 \]
這種實序列FFT算法比相同長度的復序列FFT算法大約可減小一半的運算量。

3、使用方法

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

/************************************
    x       ----長度爲n。開始時存放要變換的實數據，最後存放變換結果的前n/2+1個值，
                其存儲順序爲[Re(0),Re(1),...,Re(n/2),Im(n/2-1),...,Im(1)]。
                其中Re(0)=X(0),Re(n/2)=X(n/2)。根據X(k)的共軛對稱性，很容易寫
                出後半部分的值。
    n       ----數據長度，必須是2的整數次冪，即n=2^m。
************************************/
#include "math.h"
#include "stdlib.h"
#include "fft.c"

void rlfft(double *x, int n)
{
    int i, nl;
    double a, c, e, s, fr, fi, gr, gi, *£, *g;
    f = malloc(n / 2 * sizeof(double));
    g = malloc(n / 2 * sizeof(double));
    n1 = n / 2;
    e = 3.141592653589793 / nl;
    for(i = 0; i < n1; i++) {
        f[i] = x[2 * i];
        g[i] = x[2 * i + 1];
    }
    fft(f, g, n1, 1);
    x[0] = f[0] + g[0];
    x[n1] = f[0] - g[0];
    for(i = 1; i < n1; i++) {
        fr = (f[i] + f[n1 - i] / 2);
        fi = (g[i] - g[n1 - i] / 2);
        gr = (g[n1 - i] + g[i] / 2);
        gi = (f[n1 - i] - f[i] / 2);
        a = i * e;
        c = cos(a);
        s = sin(a);
        x[i] = fr + c * gr + s * gi;
        x[n - i] = fi + c * gi - s * gr;
    }
    free(f);
    free(g);
}

fft.c文件參見快速傅里葉變換