一、功能
计算复序列的快速傅里叶变换。
二、方法简介
序列(x(n)(n=0,1,...,N-1))的离散傅里叶变换定义为
[X(k)=sum_{n=0}^{N-1}x(n)W_{N}^{nk}, qquad k=0,1,...,N-1
]
其中(W_{N}^{nk}=e^{-jfrac{2pi nk}{N}})，将序列(x(n))按序号(n)的奇偶分成两组，即
[left.egin{matrix}egin{align*}x_{1}(n)=&x(2n)\ x_{2}(n)=&x(2n+1)end{align*}end{matrix}
ight} qquad n=0,1,...,frac{N}{2}-1
]
因此，(x(n))的傅里叶变换可写成
[egin{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})，式中
[egin{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)点的值可表示为
[egin{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}-1end{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）算法合并在同一程序中。
三、使用说明
是用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;
}
}
}