基4快速傅里叶变换
一、功能
计算复序列的基4快速傅里叶变换。
二、方法简介
序列\(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^{-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%，加法量也略有减少。
三、使用说明
是用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;
}
}
