还剩7页未读,继续阅读
文本内容:
《测试__分析及处理》课程作业快速傅里叶变换
1、程序设计思路快速傅里叶变换的目的是减少运算量,其用到的方法是分级进行运算全部计算分解为级,其中;在输入序列中是按码位倒序排列的,输出序列是按顺序排列;每级包含个蝶形单元,第级有个群,每个群有个蝶形单元;每个蝶形单元都包含乘和系数的运算,每个蝶形单元数据的间隔为,i为第i级;同一级中各个群的系数分布规律完全相同将输入序列按码位倒序排列时,用到的是倒序算法——雷德算法自然序排列的二进制数,其下面一个数总比上面的数大1,而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的若已知某数的倒序数是,求下一个倒序数,应先判断的最高位是否为0,与进行比较即可得到结果如果,说明最高位为0,应把其变成1,即,这样就得到倒序数了如果,即的最高位为1,将最高位化为0,即,再判断次高位;与进行比较,若为0,将其变位1,即,即得到倒序数,如果次高位为1,将其化为0,再判断下一位……即从高位到低位依次判断其是否为1,为1将其变位0,若这一位为0,将其变位1,即可得到倒序数若倒序数小于顺序数,进行换位,否则不变,防治重复交换,变回原数注因为0的倒序数为0,所以可从1开始进行求解
2、程序设计框图
(1)倒序算法——雷德算法流程图2FFT算法流程
3、FFT源程序voidfftxnintn;doublex[];{intijklmn1n2;doublecc1ess1ttr;forj=1i=1;in/2;i++{m=i;j=2*j;ifj==nbreak;}//得到流程图的共几级n1=n-1;forj=0i=0;in1;i++{ifij//如果ij即进行变址{tr=x[j];x[j]=x[i];x[i]=tr;}k=n/2;//求j的下一个倒位序whilekj+1//如果kj+1表示j的最高位为1{j=j-k;//把最高位变成0k=k/2;//k/2,比较次高位,依次类推,逐个比较,直到某个位为0}j=j+k;//把0改为1}fori=0;in;i+=2{tr=x[i];x[i]=tr+x[i+1];x[i+1]=tr-x[i+1];}n2=1;forl=1;l=m;l++//控制蝶形结级数{n4=n2;n2=2*n4;n1=2*n2;e=
6.28318530718/n1;fori=0;in;i+=n1//控制同一蝶形结运算,即计算系数相同蝶形结{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;forj=2;j=n4-1;j++//控制计算不同种蝶形结,即计算系数不同的蝶形结{i1=i+j;i2=i-j+n2;i3=i+j+n2;i4=i-j+n1;cc=cosa;ss=sina;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;}}}}
4、计算实例及运行结果设输入序列为其离散傅里叶变换为这里选n=512计算离散傅里叶变换所用软件为Turboc
2.0,操作界面如图1所示图1Turboc
2.0操作界面程序运行结束后的界面如图2所示图2程序运行后的界面例子的具体程序如下#include__th.h#includestdio.h#includestdlib.hvoidfftxnintn;doublex[];{intijkli1i2i3i4n4mn1n2;doubleaeccsstrt1t2;forj=1i=1;in/2;i++{m=i;j=2*j;ifj==nbreak;}n1=n-1;forj=0i=0;in1;i++{ifij{tr=x[j];x[j]=x[i];x[i]=tr;}k=n/2;whilekj+1{j=j-k;k=k/2;}j=j+k;}fori=0;in;i+=2{tr=x[i];x[i]=tr+x[i+1];x[i+1]=tr-x[i+1];}n2=1;forl=1;l=m;l++{n4=n2;n2=2*n4;n1=2*n2;e=
6.28318530718/n1;fori=0;in;i+=n1{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;forj=2;j=n4-1;j++{i1=i+j;i2=i-j+n2;i3=i+j+n2;i4=i-j+n1;cc=cosa;ss=sina;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;}}}}__in{FILE*p;intijn;doubledt=
0.001;doublex
[512];p=fopend:\
123.cw;n=512;fori=0;in;i++{x[i]=sin200*pi*i*dt;}fori=0;in;i++{fprintfp%
10.7fx[i];fprintfp\n;printf%
10.7fx[i];printf\n;}fftxn;fprintfp\nDISCRETEFOURIERTRANSFORM\n;printf\nDISCRETEFOURIERTRANSFORM\n;fprintfp%
10.7fx
[0];printf%
10.7fx
[0];fprintfp%
10.7f+J%
10.7f\nx
[1]x[n-1];printf%
10.7f+J%
10.7f\nx
[1]x[n-1];fori=2;in/2;i+=2{fprintfp%
10.7f+J%
10.7fx[i]x[n-i];fprintfp%
10.7f+J%
10.7fx[i+1]x[n-i-1];fprintfp\n;printf%
10.7f+J%
10.7fx[i]x[n-i];printf%
10.7f+J%
10.7fx[i+1]x[n-i-1];printf\n;}fprintfp%
10.7fx[n/2];printf%
10.7fx[n/2];fprintfp%
10.7f+J%
10.7f\nx[n/2-1]-x[n/2+1];fori=2;in/2;i+=2{fprintfp%
10.7f+J%
10.7fx[n/2-i]-x[n/2+i];fprintfp%
10.7f+J%
10.7fx[n/2-i-1]-x[n/2+i+1];fprintfp\n;printf%
10.7f+J%
10.7fx[n/2-i]-x[n/2+i];printf%
10.7f+J%
10.7fx[n/2-i-1]-x[n/2+i+1];printf\n;}}将程序运行后所得数据绘制成曲线图(其中FFT变换的数据要先取绝对值后再画图)如下由上图可知,变换后的图开在频率100Hz处出现一个峰值,这与理论上的结果一致。