快速傅里叶变换FFT原理及源程序.docx
- 文档编号:9984106
- 上传时间:2023-02-07
- 格式:DOCX
- 页数:9
- 大小:16.14KB
快速傅里叶变换FFT原理及源程序.docx
《快速傅里叶变换FFT原理及源程序.docx》由会员分享,可在线阅读,更多相关《快速傅里叶变换FFT原理及源程序.docx(9页珍藏版)》请在冰豆网上搜索。
快速傅里叶变换FFT原理及源程序
资料范本
本资料为word版本,可以直接编辑和打印,感谢您的下载
快速傅里叶变换(FFT)原理及源程序
地点:
__________________
时间:
__________________
说明:
本资料适用于约定双方经过谈判,协商而共同承认,共同遵守的责任与义务,仅供参考,文档可直接下载或修改,不需要的部分可直接删除,使用时请详细阅读内容
《测试信号分析及处理》课程作业
快速傅里叶变换
程序设计思路
快速傅里叶变换的目的是减少运算量,其用到的方法是分级进行运算。
全部计算分解为级,其中;在输入序列中是按码位倒序排列的,输出序列是按顺序排列;每级包含个蝶形单元,第级有个群,每个群有个蝶形单元;每个蝶形单元都包含乘和系数的运算,每个蝶形单元数据的间隔为,i为第i级;同一级中各个群的系数分布规律完全相同。
将输入序列按码位倒序排列时,用到的是倒序算法——雷德算法。
自然序排列的二进制数,其下面一个数总比上面的数大1,而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。
若已知某数的倒序数是,求下一个倒序数,应先判断的最高位是否为0,与进行比较即可得到结果。
如果,说明最高位为0,应把其变成1,即,这样就得到倒序数了。
如果,即的最高位为1,将最高位化为0,即,再判断次高位;与进行比较,若为0,将其变位1,即,即得到倒序数,如果次高位为1,将其化为0,再判断下一位……即从高位到低位依次判断其是否为1,为1将其变位0,若这一位为0,将其变位1,即可得到倒序数。
若倒序数小于顺序数,进行换位,否则不变,防治重复交换,变回原数。
注:
因为0的倒序数为0,所以可从1开始进行求解。
程序设计框图
(1)倒序算法——雷德算法流程图
(2)FFT算法流程
FFT源程序
voidfft(x,n)
intn;
doublex[];
{inti,j,k,l,m,n1,n2;
doublec,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 {tr=x[j]; 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.28318530718/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; } } } } 计算实例及运行结果 设输入序列为 其离散傅里叶变换为 这里。 选n=512,计算离散傅里叶变换。 所用软件为Turboc2.0,操作界面如图1所示 图1Turboc2.0操作界面 程序运行结束后的界面如图2所示 图2程序运行后的界面 例子的具体程序如下: #include #include #include #definepi3.14159265359 voidfft(x,n) intn; doublex[]; {inti,j,k,l,i1,i2,i3,i4,n4,m,n1,n2; doublea,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; } 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.28318530718/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; } } } } main() {FILE*p; inti,j,n; doubledt=0.001; doublex[512]; p=fopen("d: \123.c","w"); n=512; for(i=0;i {x[i]=sin(200*pi*i*dt); } for(i=0;i {fprintf(p,"%10.7f",x[i]); fprintf(p,"\n"); printf("%10.7f",x[i]); printf("\n"); } fft(x,n); fprintf(p,"\nDISCRETEFOURIERTRANSFORM\n"); printf("\nDISCRETEFOURIERTRANSFORM\n"); fprintf(p,"%10.7f",x[0]); printf("%10.7f",x[0]); fprintf(p,"%10.7f+J%10.7f\n",x[1],x[n-1]); printf("%10.7f+J%10.7f\n",x[1],x[n-1]); for(i=2;i { fprintf(p,"%10.7f+J%10.7f",x[i],x[n-i]); fprintf(p,"%10.7f+J%10.7f",x[i+1],x[n-i-1]); fprintf(p,"\n"); printf("%10.7f+J%10.7f",x[i],x[n-i]); printf("%10.7f+J%10.7f",x[i+1],x[n-i-1]); printf("\n"); } fprintf(p,"%10.7f",x[n/2]); printf("%10.7f",x[n/2]); fprintf(p,"%10.7f+J%10.7f\n",x[n/2-1],-x[n/2+1]); for(i=2;i { fprintf(p,"%10.7f+J%10.7f",x[n/2-i],-x[n/2+i]); fprintf(p,"%10.7f+J%10.7f",x[n/2-i-1],-x[n/2+i+1]); fprintf(p,"\n"); printf("%10.7f+J%10.7f",x[n/2-i],-x[n/2+i]); printf("%10.7f+J%10.7f",x[n/2-i-1],-x[n/2+i+1]); printf("\n"); } } 将程序运行后所得数据绘制成曲线图(其中FFT变换的数据要先取绝对值后再画图)如下 由上图可知,变换后的图开在频率100Hz处出现一个峰值,这与理论上的结果一致。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 快速 傅里叶变换 FFT 原理 源程序