信号与系统课程设计——FFT的计算机实现快速傅立叶变换(FFT)的计算机实现.doc

信号与系统课程设计fft的计算机实现快速傅里叶变换（fft）的计算机实现摘要：本文是信号与系统课程的课程设计，旨在熟悉fft的计算过程，结合dft物理意义和实验结果加深对傅立叶变换的理解。文章首先用matlab对一个简单信号进行fft仿真，得出频谱图；其次完成了fft的c语言实现，结合matlab作图及数据处理功能得出了c实现下的fft结果；最后，讨论分析实验结果。关键词：dft、基-2按时间抽取fft算法、matlab、c、频谱、物理意义1. 算法描述1) dft的运算量2) 减少运算的方法：v 化长序列为短序列。如将长度为n的序列分解为两个长度为n/2的序列v 利用的性质（注：本文中的c程序未用到此性质） 3) c程序采用基-2按时间抽取的fft算法设输入序列长度为 (m为正整数)，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的fft算法,也称为coolkey-tukey算法。若n不满足条件，则人为地加上若干零值，使。2. 实验设计与步骤为简单起见,同时不失一般性,本实验采用三个余弦成分和一个支流偏置成分叠加所得的信号作为信号源。3. matlab的fft仿真和c的fft实现1) 关于频谱的理论分析：2) matlab的fft仿真v matlab程序代码function output_args = fft2( x0,x1,f1,w1,x2,f2,w2,x3,f3,w3,fs )%这是一个自定义函数，输入被采样的信号的参数和采样频率，然后输出原信号波形、采样信号序列、采样序列幅度谱和相位谱。函数规定信号由一个直流成分（大小为x0），三个余弦成分（各自频率分别为f1,f2,f3,幅值为x1,x2,x3,初相位为w1,w2,w3）,采样频率为fs，采样时间为1秒。x0,x1,f1,w1,x2,f2,w2,x3,f3,w3,fs作为参数输入t=0:1/fs:1-1/fs;%定义采样时刻n=length(t);%采样序列长度s=x0+x1*cos(2*pi*f1*t+pi/180*w1)+x2*cos(2*pi*f2*t+pi/180*w2)+x3*cos(2*pi*f3*t+pi/180*w3);%采样信号y=fft(s);%快速傅立叶变换tt=0:1/10000:1;%原信号描点ss=x0+x1*cos(2*pi*f1*tt+pi/180*w1)+x2*cos(2*pi*f2*tt+pi/180*w2)+x3*cos(2*pi*f3*tt+pi/180*w3);%原信号figure;plot(tt,ss);grid;title(原始信号);%原信号波形figure;subplot(3,1,1);plot(0:n-1,s(1:n),-o);xlim(0 n-1);grid;title(采样信号);subplot(3,1,2);plot(0:n-1,abs(y(1:n),-o);xlim(0 n-1);grid;xlabel(k);ylabel(幅度);title(理想采样信号的幅度谱);subplot(3,1,3);plot(0:n-1,angle(y(1:n),-o);grid;xlabel(k);ylabel(相位);axis(0 n-1 -pi pi);title(理想采样信号的相位谱);end;v 改变采样频率 依次键入：%其中fft2()是自定义函数,对信号在1秒内以频率fs进行采样fft2(1,1,1,90,2,2,180,3,3,180,32)fft2(1,1,1,90,2,2,180,3,3,180,16)fft2(1,1,1,90,2,2,180,3,3,180,8)fft2(1,1,1,90,2,2,180,3,3,180,4)原始信号如图：图1 图2、fs =32hz图3、fs =16hz图4、fs =8hz图5、fs =4hz注意到：1、幅度谱，除却k=0外，图形呈轴对称分布；2、相位谱，n为偶数时除却k=0和n/2外（n为奇数时，除却k=0外），图形呈中心对称分布。理论上,由于复指数的周期性，长度为n（假设n为偶数，n为奇数时类似）的的dft频谱分析，可把k=n/2,n/2+1,n-1看作是负频率-n/2,-n/2+1,-n/2+2,-2,-1,特别的，对于实序列而言，正负频率成分的幅值必须相等，而初相位必须相反。下面取fs=8hz和4hz部分计算结果分析，其他情况可类似处理。可看到：fs=8hz时，由频谱结合物理意义计算的得到的余弦成分的幅值和初相位与原信号一致；fs=4hz时，由频谱结合物理意义计算的得到的余弦成分的幅值和初相位与原信号不同。对此，可根据抽样定理解释。因为f1=1hz，f2=2hz， f3=3hz，所以连续时间信号的截至频率fm=f3=3hz，所以抽样频率fs=4hz时，fs %改变信号的直流成分（由1变为2），原信号和频谱如下fft2(2,1,1,90,2,2,180,3,3,180,32)图9 图10 fs=32hz变化趋势：幅度谱中k=0对应的幅度变为2n=64(原为32)，其余处无明显变化，相位谱中，k=0,1,2,3对应的相位值都不变。 键入： %改变第三个余弦成分的幅值（由3改为4），原信号图和频谱图如下fft2(2,1,1,90,2,2,180,4,3,180,32)图11 图12 fs=32hz变化趋势：幅度谱中k=3对应的幅度变为64（4*32/2=64,原来这个值为3*32/2=48），其余处无明显变化，相位谱中k=0,1,2,3对应的相位值都未发生变化。 键入： %改变第三个余弦成分的频率（由3hz改为4hz），原信号图和频谱如下fft2(2,1,1,90,2,2,180,4,4,180,32)图13 图14 fs=32hz变化趋势：幅度谱中，k=4处幅度为n/2*4=64,k=3处幅度变为0，k=0,1,2处幅度不变；相位谱中，k=4处相位变为180，k=0,1,2,3处相位值不变。 键入：%改变第三个余弦成分的初相位（由180改为270），原信号和频谱如下图fft2(2,1,1,90,2,2,180,4,4,270,32)图15 图16 fs=32hz变化趋势：幅度谱波形无明显变化；相位谱中，k=4处相位值变为-90度，即270度。由此可见，改变信号的波幅、频率和相位，幅度谱和相位谱将发生相应的变化，即恰合频谱的物理意义。3) c语言实现fftv c代码/*此代码用于fft计算，采用基-2按时间抽取的fft算法decimation-in-time(dit)(coolkey-tukey)。为方便起见，同时不失一般性，把含有一个直流成分和三个正弦成分的信号作为被采样信号，信号的输入由函数input()完成，若想改变信号波形只需改变input()函数代码。代码分别定义了复数结构体、复数运算和码位倒读函数。计算结果分别在命令窗口和文件e:keshe.txt中输出，keshe.txt文件数据可用于matlab作图分析*/#include#include/*圆周率*/double pi=3.141592653589793;int m,n;/*定义复数*/struct complex_double real;double img;*a,*b;/*定义2的幂计算*/int x_2(int a)int i,r;r=1;for(i=0;i=0;j-)ii+=(int)(i/x_2(j)*x_2(m-1-j);/*先把i用二进制表示，然后码位倒读*/i-=(int)(i/x_2(j)*x_2(j);return(ii);/*定义复数乘法*/struct complex_ multi(struct complex_ a,struct complex_ b)struct complex_ c;c.real=a.real*b.real-a.img*b.img;c.img=a.real*b.img+a.img*b.real;return(c);/*定义复数幂运算*/struct complex_ w_k(struct complex_ w,int k)int i;struct complex_ x,y;x=w;if(k=0)x.real=1;x.img=0;return(x);elsefor(i=1;ik;i+)y=multi(x,w);x=y;return(x);/*定义复数加法*/struct complex_ add(struct complex_ a,struct complex_ b)struct complex_ c;c.real=a.real+b.real;c.img=a.img+b.img;return(c);/*离散序列的输入*/void input()int nn,i;float f1,f2,f3,a0,a1,a2,a3,w1,w2,w3;m=0;/*在一秒钟内，采点数*/printf(采点数:);scanf(%d,&nn);printf(n);/*信号的直流成分*/printf(信号的直流成分:n);scanf(%f,&a0);/*信号有三个余弦成分*/*余弦成分1*/printf(第1个正弦成分的频率:n);scanf(%f,&f1);printf(幅值:n);scanf(%f,&a1);printf(初相位：n);scanf(%f,&w1);/*余弦成分2*/printf(第2个正弦成分的频率:n);scanf(%f,&f2);printf(幅值:n);scanf(%f,&a2);printf(初相位：n);scanf(%f,&w2);/*余弦成分3*/printf(第3个正弦成分的频率:n);scanf(%f,&f3);printf(幅值:n);scanf(%f,&a3);printf(初相位：n);scanf(%f,&w3);for(i=0;x_2(i)nn;i+)m=i+1;n=x_2(m);/*n为不小于nn的最小的2的幂*/a=(struct complex_*)calloc(n,sizeof(struct complex_);/*动态开辟n个单位*/b=(struct complex_*)calloc(n,sizeof(struct complex_);/*动态开辟n个单位*/*把采样信号用码位倒读的方法存入计算机中*/for(i=0;inn;i+)areverse(i).real=a0+a1*cos(2*pi*f1*i/n+w1/180*pi)+a2*cos(2*pi*f2*i/n+w2/180*pi)+a3*cos(2*pi*f3*i/n+w3/180*pi);/*其中reverse()为码位倒读函数*/areverse(i).img=0;/*动态空间空闲处，补零*/for(i=nn;in;i+)areverse(i).real=0;areverse(i).img=0;/*主函数*/void main()int i,j,k,m,n,q,r;file *fp;/*文件指针，指向存储fft结果的文件*/struct complex_ d,w,wk,x,y;j=0;input();m=1;n=n;for(i=0;im;i+)/*m级运算*/n=n/2;m=m*2;q=x_2(i+1);/*q=2(i+1)*/w.real=cos(-2*pi/q);/*旋转因子实部*/w.img=sin(-2*pi/q);/*旋转因子虚部*/ for(j=0;jn;j+)/*n群运算*/ for(k=0;km/2;k+)wk=w_k(w,k); bj*m+k=add(aj*m+k,multi(wk,aj*m+k+m/2);/*add()为复数相加函数，multi()为复数相乘函数*/ for(k=m/2;km;k+) wk=w_k(w,k); bj*m+k=add(aj*m+k-m/2,multi(wk,aj*m+k); /*add()为复数相加函数，multi()为复数相乘函数*/ for(r=0;rn;r+) ar=br; /*将结果写入文件e:keshe.txt*/if(fp=fopen(e:keshe.txt,w+)=null)printf(cannot open file strike any key exit!);getch();exit(1);for(n=0;nn;n+)fprintf(fp,%16.14f %16.14fn,an.real,an.img);for(n=0;nn;n+)printf(%16.14f+%16.14fin,an.real,an.img);free(a);/*释放内存空间*/free(b);/*释放内存空间*/v 改变采样频率实验数据将由keshe.txt调入到matlab中，做出频谱图。 参数设置1图17 fs=32hz(c实现)图18 fs=32hz（matlab仿真）注意到两图的幅度谱图，很相似，而相位谱却相差很大。用c实现的fft得到的相位谱也似乎没有什么对称性，但可以看到两个相位谱在k=0,1,2,3和31,30,29处的相位值是基本一致的，而其他点处，相位值差异则很大。（表1，列出了具体的数值计算结果）表1kfft（c）fft（matlab）fft结果相位值fft结果相位值032.0000000000000 + 0.00000000000000i032.0000000000000 + 0.00000000000000i010.00000000000000 + 16.0000000000000i1.5707962.99760216648792e-15 - 16.0000000000000i1.5707962-32.0000000000000 + 1.00000000000000e-14i3.141593-32.0000000000000 - 9.76996261670138e-15i3.1415933-48.0000000000000 + 3.00000000000000e-14i3.141593-48.0000000000000 - 4.04121180963557e-14i3.1415934-1.00000000000000e-14 + 0.00000000000000i3.141593-1.16001965116118e-14 - 9.42055475210265e-16i3.0605650.00000000000000 + 0.00000000000000i0-3.99680288865056e-15 - 1.77635683940025e-15i2.72336861.00000000000000e-14 + 1.00000000000000e-14i0.7853987.54951656745106e-15 - 4.88498130835069e-15i0.57430571.00000000000000e-14 + 1.00000000000000e-14i0.7853981.11022302462516e-15 - 5.21804821573824e-15i1.3611568-1.00000000000000e-14 + 1.00000000000000e-14i2.356194-1.06581410364015e-14 - 7.10542735760100e-15i2.5535990.00000000000000 - 1.00000000000000e-14i-1.57081.11022302462516e-15 + 7.21644966006352e-15i-1.41815100.00000000000000 + 0.00000000000000i0-2.22044604925031e-15 - 2.22044604925031e-15i2.356194110.00000000000000 + 0.00000000000000i01.11022302462516e-14 - 5.32907051820075e-15i0.4475212-1.00000000000000e-14 + 0.00000000000000i3.141593-9.71608556119124e-15 - 9.42055475210265e-16i3.044936130.00000000000000 + 0.00000000000000i03.55271367880050e-15 - 1.33226762955019e-15i0.35877114-1.00000000000000e-14 + 0.00000000000000i3.1415933.55271367880050e-15 + 0.00000000000000i015-2.00000000000000e-14 + 1.00000000000000e-14i2.677945-1.58761892521397e-14 - 1.15463194561016e-14i2.512796161.00000000000000e-14 + 0.00000000000000i0-3.55271367880050e-15 + 0.00000000000000i3.14159317-2.00000000000000e-14 - 1.00000000000000e-14i-2.67795-1.58761892521397e-14 + 1.15463194561016e-14i-2.512818-1.00000000000000e-14 + 0.00000000000000i3.1415933.55271367880050e-15 + 0.00000000000000i019-1.00000000000000e-14 + 0.00000000000000i3.1415933.55271367880050e-15 + 1.33226762955019e-15i-0.3587720-1.00000000000000e-14 + 0.00000000000000i3.141593-9.71608556119124e-15 + 9.42055475210265e-16i-3.04494211.00000000000000e-14 + 0.00000000000000i01.11022302462516e-14 + 5.32907051820075e-15i-0.44752220.00000000000000 + 0.00000000000000i0-2.22044604925031e-15 + 2.22044604925031e-15i-2.35619230.00000000000000 + 1.00000000000000e-14i1.5707961.11022302462516e-15 - 7.21644966006352e-15i1.41814724-1.00000000000000e-14 - 1.00000000000000e-14i-2.35619-1.06581410364015e-14 + 7.10542735760100e-15i-2.55359250.00000000000000 - 1.00000000000000e-14i-1.57081.11022302462516e-15 + 5.21804821573824e-15i-1.36116261.00000000000000e-14 - 1.00000000000000e-14i-0.78547.54951656745106e-15 + 4.88498130835069e-15i-0.574327-1.00000000000000e-14 + 0.00000000000000i3.141593-3.99680288865056e-15 + 1.77635683940025e-15i-2.7233728-1.00000000000000e-14 + 0.00000000000000i3.141593-1.16001965116118e-14 + 9.42055475210265e-16i-3.0605629-48.0000000000000 - 5.00000000000000e-14i-3.14159-48.0000000000000 + 4.04121180963557e-14i-3.1415930-32.0000000000000 - 2.00000000000000e-14i-3.14159-32.0000000000000 + 9.76996261670138e-15i-3.14159311.00000000000000e-14 - 16.0000000000000i-1.57082.99760216648792e-15 + 16.0000000000000i-1.5708结合图17、18和表1，可以发现两个相位谱在k=0,1,2,3和31,30,29处一致，而在其他点相差很大，注意到这些点的幅值很小（理论上应该为0，实际上却由于计算机的误差，而使之非0），可认为为0，可忽略这些点，而主要考虑主频率成分。事实上，结合图表可发现，c程序所得的结果也并非完全无对称性的，并且可注意到所有不对称点的相位值接近圆周率或0。这种差异是由于计算机的误差导致的，在一般情况下，这不会造成太大的影响，而在相位处于边界值（和0）时，这种误差将计算结果发生很大变化。注：在先前编写的c程序中，把复数实部和虚部的数据类型设置为float时，相位谱的不对称更加严重，甚至c和matlab仿真的两个相位谱在k= 31,30,29都不一致，此处实验结果是把数据类型改为double类型后（即提高精度后）所得的，仍存在这种不对称，理论上可通过不断提高精度进行改进，但只要是某频率成分初相位为或0，就存在这种严重不对称的风险。 参数设置2图19 fs=16hz(c实现) 图20 fs=16hz（matlab仿真） 可发现c实现下的fft幅度谱图随采样频率变化的趋势与matlab仿真实验类似，而由于参数1实验相同的原因，c实现的fft相位谱并不对称。 参数设置3图21 fs=8hz(c实现) 图22 fs=8hz（matlab仿真）类似的，可发现幅度谱图随采样频率变化的趋势与matlab仿真实验类似。而c和matlab实现下，k=6时的相位值一个为3.14，一个为-3.14，具体原因同上。而和相差，理论上初相位为是无差异的，但是体现在相位谱上就是造成不对称。 参数设置4图23 fs=4hz(c实现)图24 fs=4hz（matlab仿真）同matlab仿真结果一样出现混叠现象。 参数设置5图25 fs=14hz(c实现) 图26 fs=14hz（matlab仿真）注意到幅度谱和相位谱都有很大区别，其中相位谱的区别是由两种不同的计算机语言规则造成的（正如之前提到的），而幅度谱的区别是由算法不同造成的， c程序是采用补零的方法把序列长度凑成，然后把长度为n的序列输出为结果，这样做的后果是改变了进行dft的序列， 得到的自然不是所求的结果，所以这个c程序，只能针对序列长度满足的情况。v 改变被采样信号波形 参数设置6图27 fs=32hz(c实现)图28 fs=32hz(matlab仿真)二者幅度谱相同，相位谱只在中间部分有差异（此差异对实际分析影响很小，且也无实际意义，可忽略），而在k=0,1,2,3和31,30,29处两幅图的情况一致，故验证了之前关于出现图17和图18二者相位谱如此大差异原因的猜想。 参数设置7图29 fs=32hz(c实现)图30 fs=32hz(matlab仿真)二者幅度谱图相同，且较上一参数时发生的变化趋势与matlab仿真试验中一样；相位谱图在k=0,1,2,3和31,30,29处相位值相同。 参数设置8图31 fs=32hz(c实现)图32 fs=32hz(matlab仿真)二者幅度谱图相同，且较上一参数时发生的变化趋势与matlab仿真试验中一样；相位谱图在k=0,1,2,4和31,30,28处相位值相同。至此，考察了序列长度和信号谐波成分对matlab的fft仿真和c实现的fft的影响，并结合了dft的物理意义，对实验现象做出了解释。对比matlab和c实现下的fft结果的差异，做出了自己的解释。4. 课程设计心得与自我评价这个课程设计的主要部分是在暑假完成的，课程设计是关于fft的计算机实现，要求掌握fft算法，并分别用matlab进行fft仿真和编写实现fft的c代码。记得之前学习电路理论、复变函数、数理方程和信号系统时都遇到过把时域上的问题转化为频域上的问题解决的情况，电路中主要是laplace变换的实际应用，而复变函数里面主要讲的是各种变换的技术新性问题，数理方程里面研究的是一些数学方程的物理模型，而信号与系统课程则引入了离散的概念。首先是时域上的离散，由对连续时间进行信号采样，即x(t)xn，把连续时间信号的各种性质在离散时间信号做推广，如ctftdtft，laplace变换z-变换；其次是频域上的离散，对连续频率进行信号采样，即dtftdft，而fft则是处理dft的一种快速算法。课程设计完成的过程中遇到许多问题：首先是matlab