1、快速傅里叶变换FFT原理及源程序之欧阳育创编测试信号分析及处理课程作业时间:2021.02.04创作:欧阳育快速傅里叶变换一、程序设计思路快速傅里叶变换的目的是减少运算量,其用到的方法是分级进行运算。全部计算分解为级,其中;在输入序列中是按码位倒序排列的,输出序列是按顺序排列;每级包含个蝶形单元,第级有个群,每个群有个蝶形单元; 每个蝶形单元都包含乘和系数的运算,每个蝶形单元数据的间隔为,i为第i级; 同一级中各个群的系数分布规律完全相同。将输入序列按码位倒序排列时,用到的是倒序算法雷德算法。 自然序排列的二进制数,其下面一个数总比上面的数大1,而倒序二进制数的下面一个数是上面一个数在最高位加
2、1并由高位向低位仅为而得到的。 若已知某数的倒序数是,求下一个倒序数,应先判断的最高位是否为0,与进行比较即可得到结果。如果,说明最高位为0,应把其变成1,即,这样就得到倒序数了。如果,即的最高位为1,将最高位化为0,即,再判断次高位;与进行比较,若为0,将其变位1,即,即得到倒序数,如果次高位为1,将其化为0,再判断下一位即从高位到低位依次判断其是否为1,为1将其变位0,若这一位为0,将其变位1,即可得到倒序数。若倒序数小于顺序数,进行换位,否则不变,防治重复交换,变回原数。注:因为0的倒序数为0,所以可从1开始进行求解。二、程序设计框图(1)倒序算法雷德算法流程图(2)FFT算法流程三、F
3、FT源程序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;in/2;i+) m=i; j=2*j; if(j=n)break; /得到流程图的共几级 n1=n-1; for(j=0,i=0;in1;i+) if(ij) /如果ij,即进行变址 tr=xj; xj=xi; xi=tr; k=n/2; /求j的下一个倒位序 while(k(j+1) /如果k(j+1),表示j的最高位为1 j=j-k; /把最高位变成0 k=k/2; /k/2,比较次高位,依次类推,逐个比较,
4、直到某个位为0 j=j+k; /把0改为1 for(i=0;in;i+=2) tr=xi; xi=tr+xi+1; xi+1=tr-xi+1; n2=1; for(l=1;l=m;l+) / 控制蝶形结级数 n4=n2; n2=2*n4; n1=2*n2; e=6.28318530718/n1; for(i=0;in;i+=n1) /控制同一蝶形结运算,即计算系数相同蝶形结 tr=xi; xi=tr+xi+n2; xi+n2=tr-xi+n2; xi+n2+n4=-xi+n2+n4; a=e; for(j=2;j=(n4-1);j+) /控制计算不同种蝶形结,即计算系数不同的蝶形结 i1=i+
5、j; i2=i-j+n2; i3=i+j+n2; i4=i-j+n1; cc=cos(a); ss=sin(a); a=a+e; t1=cc*xi3+ss*xi4; t2=ss*xi3-cc*xi4; xi4=xi2-t2; xi3=-xi2-t2; xi2=xi1-t1; xi1=xi1+t1; 四、计算实例及运行结果设输入序列为其离散傅里叶变换为这里。选n=512,计算离散傅里叶变换。所用软件为Turbo c 2.0,操作界面如图1所示图1 Turbo c 2.0操作界面程序运行结束后的界面如图2所示图2 程序运行后的界面例子的具体程序如下:#include#include#include
6、#define pi 3.14159265359void 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;in/2;i+) m=i; j=2*j; if(j=n)break; n1=n-1; for(j=0,i=0;in1;i+) if(ij) tr=xj; xj=xi; xi=tr; k=n/2; while(k(j+1) j=j-k; k=k/2; j=j+k; for(i=0;in;i+=2) tr=xi; xi=tr+xi+1; xi+1=
7、tr-xi+1; n2=1; for(l=1;l=m;l+) n4=n2; n2=2*n4; n1=2*n2; e=6.28318530718/n1; for(i=0;in;i+=n1) tr=xi; xi=tr+xi+n2; xi+n2=tr-xi+n2; xi+n2+n4=-xi+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*xi3+ss*xi4; t2=ss*xi3-cc*xi4; xi4=xi2-t2; xi3=-x
8、i2-t2; xi2=xi1-t1; xi1=xi1+t1; main()FILE *p; int i,j,n; double dt=0.001; double x512; p=fopen(d:123.c,w); n=512; for(i=0;in;i+) xi=sin(200*pi*i*dt); for(i=0;in;i+) fprintf(p,%10.7f,xi); fprintf(p,n); printf(%10.7f,xi); printf(n); fft(x,n); fprintf(p,n DISCRETE FOURIER TRANSFORMn);printf(n DISCRETE
9、FOURIER TRANSFORMn); fprintf(p,%10.7f,x0); printf(%10.7f,x0); fprintf(p,%10.7f+J%10.7fn,x1,xn-1); printf(%10.7f+J%10.7fn,x1,xn-1); for(i=2;in/2;i+=2) fprintf(p,%10.7f+J%10.7f,xi,xn-i); fprintf(p,%10.7f+J%10.7f,xi+1,xn-i-1); fprintf(p,n); printf(%10.7f+J%10.7f,xi,xn-i); printf(%10.7f+J%10.7f,xi+1,xn-
10、i-1); printf(n); fprintf(p,%10.7f,xn/2); printf(%10.7f,xn/2); fprintf(p,%10.7f+J%10.7fn,xn/2-1,-xn/2+1); for(i=2;in/2;i+=2) fprintf(p,%10.7f+J%10.7f,xn/2-i,-xn/2+i); fprintf(p,%10.7f+J%10.7f,xn/2-i-1,-xn/2+i+1); fprintf(p,n); printf(%10.7f+J%10.7f,xn/2-i,-xn/2+i); printf(%10.7f+J%10.7f,xn/2-i-1,-xn/2+i+1); printf(n); 将程序运行后所得数据绘制成曲线图(其中FFT变换的数据要先取绝对值后再画图)如下由上图可知,变换后的图开在频率100Hz处出现一个峰值,这与理论上的结果一致。时间:2021.02.04创作:欧阳育
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1