FFT简介(来源：百度百科)：
快速傅立叶变换(FFT)是离散傅立叶(DFT)的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m),即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
Iterative_FFT
参数 n ——a的元素个数
a[0,1,…,(n-1)]——系数
局部变量
wd——单位元的n次原根
w——当前点的估值多项式
y——变换结果
y = BitReverse(a)
for j = 1 to logn
d =
wd =
w = 1
for k = 0 to d/2-1
for m = k to n-1 step d
t = w×y[m + d/2]
x = y[m]
y[m] = x + t
y[m + d/2] = x – t
endfor
w = w×wd
endfor
endfor
return y
C语言实现：
void Iterative_FFT(double *ir,double *ii,double *or,double
*oi,int len)
{
int i,j,n,d,k,m;
double
t;
double wdr,wdi,wr,wi;
double wtr,wti;
double tr,ti;
double xr,xi;
n=pow2(len);
for(i=0;i
{
j = BitReverse(i,len);
or[j] = ir[i];
oi[j] = ii[i];
}
for(j=1;j<=len;j++)
{
d = pow2(j);
t = 2*pi/d;
wdr = cos(t); wdi = sin(t);
wr = 1; wi = 0;
for (k=0;k<=d/2-1;k++)
{
for (m=k;m<=n-1;m+=d)
{
tr = wr*or[m + d/2] - wi*oi[m + d/2];
ti = wr*oi[m + d/2] + wi*or[m + d/2];
xr = or[m]; xi = oi[m];
or[m] = xr + tr; oi[m] = xi + ti;
or[m+d/2] = xr - tr; oi[m+d/2] = xi -
ti;
}
wtr = wr*wdr - wi*wdi;
wti = wr*wdi + wi*wdr;
wr = wtr;
wi = wti;
}
}
}
FFT的串行递归实现伪代码：
Recursive_FFT
参数 n ——a的元素个数
a[0,1,…,(n-1)]——系数
局部变量
wn——单位元的n次原根
w——当前点的估值多项式
a0——偶数项系数
a1——奇数项系数
y——变换结果
y0——a0的FFT结果
y1——a1的FFT结果
if n = =1 then return a
else
wn =
w = 1
a0 = (a[0],a[2],…a[n-2])
a1 = (a[2],a[3],…a[n-1])
y0 = Recursive_FFT(a0,n/2)
y1 = Recursive_FFT(a1,n/2)
for k = 0 to n/2-1
y[k] = y0[k] + w×y1[k]
y[k +
n/2] = y0[k]－w×y1[k]
w =
w×wn
endfor
return y
endif
C语言实现：
void Recursive_FFT(double *ir,double *ii,double *or,double
*oi,int len)
{
int i,k;
int n;
double wnr,wni;
double wr,wi;
double wtr,wti;
double temp0_ir[MAX_N], temp0_ii[MAX_N], temp1_ir[MAX_N],
temp1_ii[MAX_N];
double temp0_or[MAX_N], temp0_oi[MAX_N], temp1_or[MAX_N],
temp1_oi[MAX_N];
memset(temp0_ir, 0 ,sizeof(temp0_ir));
memset(temp0_ii, 0 ,sizeof(temp0_ii));
memset(temp1_ir, 0 ,sizeof(temp1_ir));
memset(temp1_ii, 0 ,sizeof(temp1_ii));
memset(temp0_or, 0 ,sizeof(temp0_oi));
memset(temp0_or, 0 ,sizeof(temp0_oi));
n=pow2(len);
if(len == 0)
{
or = ir;
oi = ii;
}
else
{
wnr =
cos(2*pi/n);
wni = sin(2*pi/n);
//printf("wnr:%.2f wni:%.2f\n",wnr,wni);
wr = 1;
wi = 0;
for(i=0; i
{
temp0_ir[i/2] = ir[i];
temp0_ii[i/2] = ii[i];
temp1_ir[i/2] = ir[i+1];
temp1_ii[i/2] = ii[i+1];
}
//for(i=0; i<2; i++)
//{
// printf("%.2f %.2f %.2f
%.2f\n",temp0_ir[i],temp0_ii[i],temp1_ir[i],temp1_ii[i]);
//}
Recursive_FFT(temp0_ir, temp0_ii, temp0_or, temp0_oi, len-1);
Recursive_FFT(temp1_ir, temp1_ii, temp1_or, temp1_oi, len-1);
for(k=0; k<=n/2-1; k++)
{
or[k] = temp0_or[k] + (wr*temp1_or[k] - wi*temp1_oi[k]);
oi[k] = temp0_oi[k] + (wr*temp1_oi[k] + wi*temp1_or[k]);
or[k+n/2] = temp0_or[k] - (wr*temp1_or[k] - wi*temp1_oi[k]);
oi[k+n/2] = temp0_oi[k] - (wr*temp1_oi[k] + wi*temp1_or[k]);
wtr = wr*wnr - wi*wni;
wti = wr*wni + wi*wnr;
wr = wtr;
wi = wti;
}
}
}
京公网安备 11010802041100号