一、功能
用(N)点复序列快速傅立叶变换来计算(2N)点实序列的离散傅立叶变换。
二、方法简介
假设(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{egin{matrix}egin{align*}f(n)&=x(2n)\ g(n)&=x(2n+1)end{align*}end{matrix}
ight. , 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{egin{matrix}egin{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}
ight. , n=0,1,...,N-1
]
求得(F(k))和(G(k))后，利用下面的蝶形运算计算(x(n))的离散傅立叶变换(X(k))
[left{egin{matrix}egin{align*}X(k)&=F(k)+G(k)W_{2N}^{k}\ X(k+n)&=F(k)-G(k)W_{2N}^{k}end{align*}end{matrix}
ight. , n=0,1,...,N-1
]
这种实序列FFT算法比相同长度的复序列FFT算法大约可减少一半的运算量。
三、使用方法
是用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文件参见快速傅里叶变换