快速傅里叶变换(FFT)原理及源程序

《测试信号分析及处理》课程作业

快速傅里叶变换

一、程序设计思路

快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M 级，其中N M 2log =；在输入序列()i x 中是按码位倒序排列的，输出序列()k X 是按顺序排列；每级包含2N 个蝶形单元，第i 级有i N 2

个群，每个群有12-i 个蝶形单元； 每个蝶形单元都包含乘r N W 和r N W -系数的运算，每个蝶形

单元数据的间隔为12-i ，i 为第i 级； 同一级中各个群的系数W 分布规律完全相同。

将输入序列()i x 按码位倒序排列时，用到的是倒序算法——雷德算法。 自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。

若已知某数的倒序数是J ，求下一个倒序数，应先判断J 的最高位是否为0，与2

N k =进行比较即可得到结果。如果J k >，说明最高位为0，应把其变成1，即2

N J +，这样就得到倒序数了。如果J k ≤，即J 的最高位为1，将最高位化为0，即2N J -，再判断次高位；与4N k =进行比较，若为0，将其变位1，即4N J +，即得到倒序数，如果次高位为1，将其化为0，再判断下一位……即从高位到低位依次判断其是否为1，为1将其变位0，若这一位为0，将其变位1，即可得到倒序数。若倒序数小于顺序数，进行换位，否则不变，防治重复交换，变回原数。 注：因为0的倒序数为0，所以可从1开始进行求解。

二、程序设计框图

（1）倒序算法——雷德算法流程图

(2)FFT算法流程

三、FFT源程序

void fft(x,n)

int n;

double x[];

{int i,j,k,l,m,n1,n2;

double c,c1,e,s,s1,t,tr;

for(j=1,i=1;i

{ m=i;

j=2*j;

if(j==n)break;

} //得到流程图的共几级n1=n-1;

for(j=0,i=0;i

{if(i

x[j]=x[i];

x[i]=tr;

}

k=n/2; //求j的下一个倒位序

while(k<(j+1)) //如果k<(j+1),表示j的最高位为1

{j=j-k; //把最高位变成0

k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 }

j=j+k; //把0改为1 }

for(i=0;i

{tr=x[i];

x[i]=tr+x[i+1];

x[i+1]=tr-x[i+1];

}

n2=1;

for(l=1;l<=m;l++) // 控制蝶形结级数

{n4=n2;

n2=2*n4;

n1=2*n2;

e=6./n1;

for(i=0;i

{tr=x[i];

x[i]=tr+x[i+n2];

x[i+n2]=tr-x[i+n2];

x[i+n2+n4]=-x[i+n2+n4];

a=e;

for(j=2;j<=(n4-1);j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结

{i1=i+j;

i2=i-j+n2;

i3=i+j+n2;

i4=i-j+n1;

cc=cos(a);

ss=sin(a);

a=a+e;

t1=cc*x[i3]+ss*x[i4];

t2=ss*x[i3]-cc*x[i4];

x[i4]=x[i2]-t2;

x[i3]=-x[i2]-t2;

x[i2]=x[i1]-t1;

x[i1]=x[i1]+t1;

}

}

}

}

四、计算实例及运行结果

设输入序列)(i x 为

)1,...,2,1,0(,200sin )(-=∆=n i t i i x π

其离散傅里叶变换为

)1,...,2,1,0(,)()(1

0-==∑-=n k W i x k X N i ik N 这里n j e W π

2-=。选n=512,计算离散傅里叶变换)(k X 。

所用软件为Turbo c 2.0，操作界面如图1所示

图1 Turbo c 2.0操作界面

程序运行结束后的界面如图2所示

图2 程序运行后的界面例子的具体程序如下：

#include

#include

#include

#define pi 3.

void fft(x,n)

int n;

double x[];

{int i,j,k,l,i1,i2,i3,i4,n4,m,n1,n2;

double a,e,cc,ss,tr,t1,t2;

for(j=1,i=1;i

{ m=i;

j=2*j;

if(j==n)break;

}

n1=n-1;

for(j=0,i=0;i

{if(i

{tr=x[j];

x[j]=x[i];

x[i]=tr;

}

k=n/2;

while(k<(j+1))

{j=j-k;

k=k/2;

}

j=j+k;

}