快速傅里叶变换(FFT)基础学习知识原理及其源程序

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

快速傅里叶变换

一、程序设计思路

快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M级，其中Mlog2N；在输入序列xi中是按码位倒序排列的，输出序列Xk是按顺序排列；每级包含

NN个蝶形单元，第i级有i个群，每个22rr群有2i1个蝶形单元； 每个蝶形单元都包含乘WN和WN系数的运算，每个蝶形

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

将输入序列xi按码位倒序排列时，用到的是倒序算法——雷德算法。 自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。 若已知某数的倒序数是J，求下一个倒序数，应先判断J的最高位是否为0，

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

2N即J，这样就得到倒序数了。如果kJ，即J的最高位为1，将最高位化为

2NNN0，即J，再判断次高位；与k进行比较，若为0，将其变位1，即J，

244即得到倒序数，如果次高位为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} //得到流程图的共几级 n1=n-1;

for(j=0,i=0;i{if(i-!

} 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;ifor(l=1;l<=m;l++) // 控制蝶形结级数 {n4=n2; n2=2*n4; n1=2*n2; e=6.28318530718/n1; for(i=0;i-!

}

四、计算实例及运行结果

设输入序列x(i)为

x(i)sin200it,(i0,1,2,...,n1)

其离散傅里叶变换为

ikX(k)x(i)WN,(k0,1,2,...,n1)

i0N1这里Wej2n。选n=512,计算离散傅里叶变换X(k)。

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

图1 Turbo c 2.0操作界面

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

-!

图2 程序运行后的界面

例子的具体程序如下： #include #include #include

#define pi 3.14159265359 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;in1=n-1;

for(j=0,i=0;i-!

for(i=0;ifor(l=1;l<=m;l++) {n4=n2; n2=2*n4; n1=2*n2; e=6.28318530718/n1; for(i=0;imain() {FILE *p; int i,j,n;

double dt=0.001; double x[512];

p=fopen(\"d:\\123.c\ n=512;

for(i=0;i{x[i]=sin(200*pi*i*dt);

-!

}

for(i=0;i{ fprintf(p,\"%10.7f\ fprintf(p,\"\\n\"); printf(\"%10.7f\ printf(\"\\n\"); }

fft(x,n);

fprintf(p,\"\\n DISCRETE FOURIER TRANSFORM\\n\"); printf(\"\\n DISCRETE FOURIER TRANSFORM\\n\"); fprintf(p,\"%10.7f\ printf(\"%10.7f\ fprintf(p,\"%10.7f+J%10.7f\\n\ printf(\"%10.7f+J%10.7f\\n\ for(i=2;i将程序运行后所得数据绘制成曲线图（其中FFT变换的数据要先取绝对值后再画图）如下

-!

由上图可知，变换后的图开在频率100Hz处出现一个峰值，这与理论上的结果一致。