快速傅立叶变换FFT的计算机实现.doc

《快速傅立叶变换FFT的计算机实现.doc》由会员分享，可在线阅读，更多相关《快速傅立叶变换FFT的计算机实现.doc（12页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流快速傅立叶变换FFT的计算机实现.精品文档.信号与系统课程设计FFT的计算机实现快速傅里叶变换（FFT）的计算机实现赖智鹏华中科技大学电气与电子工程学院0809班U200811806Email: 592425891摘要：本文是信号与系统课程的课程设计，旨在熟悉FFT的计算过程，结合DFT物理意义和实验结果加深对傅立叶变换的理解。文章首先用MATLAB对一个简单信号进行FFT仿真，得出频谱图；其次完成了FFT的C语言实现，结合MATLAB作图及数据处理功能得出了C实现下的FFT结果；最后，讨论分析实验结果。关键词：DFT、基-2按时间抽取FFT

2、算法、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的F

3、FT实现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

4、/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(原始信号);%原信号波

5、形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 改变采样频率 依次键入：%其中FF

6、T2()是自定义函数,对信号在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（假

7、设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时，f

8、s %改变信号的直流成分（由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对应的相位值

9、都未发生变化。 键入： %改变第三个余弦成分的频率（由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度，即2

10、70度。由此可见，改变信号的波幅、频率和相位，幅度谱和相位谱将发生相应的变化，即恰合频谱的物理意义。3) C语言实现FFTv C代码/*此代码用于fft计算，采用基-2按时间抽取的FFT算法Decimation-in-Time(DIT)(Coolkey-Tukey)。为方便起见，同时不失一般性，把含有一个直流成分和三个正弦成分的信号作为被采样信号，信号的输入由函数input()完成，若想改变信号波形只需改变input()函数代码。代码分别定义了复数结构体、复数运算和码位倒读函数。计算结果分别在命令窗口和文件e:keshe.txt中输出，keshe.txt文件数据可用于matlab作图分析*/#

11、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,st

12、ruct 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_ A

13、dd(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(第

14、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+)

15、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+w

16、3/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*

17、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

18、: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中，做出频谱图。 参数设

19、置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.0000

20、000000000i1.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.000000000

21、00000i3.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.000000000000

22、00e-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.0000

23、0000000000i0-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.

24、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.5527

25、1367880050e-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.552713

26、67880050e-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

27、-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.1102

28、2302462516e-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.1

29、41593-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

30、-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）时，这种误差将计算结果发生很大变化。注：在先前编写

31、的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实现

32、) 图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仿真）注意到幅度谱和相位谱都有很大区别，其中相位谱的区别是由两种不同的计算机语言规则造成的（正如之前提到的），而幅度谱的区别

33、是由算法不同造成的， 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(M

34、ATLAB仿真)二者幅度谱图相同，且较上一参数时发生的变化趋势与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. 课程设计心得与自我评价这

35、个课程设计的主要部分是在暑假完成的，课程设计是关于FFT的计算机实现，要求掌握FFT算法，并分别用MATLAB进行FFT仿真和编写实现FFT的C代码。记得之前学习电路理论、复变函数、数理方程和信号系统时都遇到过把时域上的问题转化为频域上的问题解决的情况，电路中主要是Laplace变换的实际应用，而复变函数里面主要讲的是各种变换的技术新性问题，数理方程里面研究的是一些数学方程的物理模型，而信号与系统课程则引入了离散的概念。首先是时域上的离散，由对连续时间进行信号采样，即x(t)xn，把连续时间信号的各种性质在离散时间信号做推广，如CTFTDTFT，Laplace变换Z-变换；其次是频域上的离散，

36、对连续频率进行信号采样，即DTFTDFT，而FFT则是处理DFT的一种快速算法。课程设计完成的过程中遇到许多问题：首先是MATLAB的仿真。这时候问题不大，直接利用现成的内置函数。其次是FFT的C实现。此处在代码的编写和程序的调试耗费了较多的时间，主要问题在于对算法的理解不够以及对C语言的一些问题未完全搞懂。后来在C实现的FFT结果与MATLAB的仿真结果对比时，遇到较多问题：1、初相位为180度（或-180度）或0的频率成分在C实现的频谱中，与MATLAB仿真结果可能有很大的出入，体现为二者相位相差一个正负号；2、C实现的结果与MATLAB仿真结果在非主频率成分下的初相位相差很大。而经过仔细

37、分析得出结论：1、问题1中初相位为180度（或-180度）时出现的异常与精度选取有关，实际应用中要么改进计算精度以消除这一异常，或者避免出现初相位为180（或-180）度的情况；2、问题2中出现的问题也是精度选取和计算机具体计算过程造成的，实际运用中由于这些频率成分的幅值可忽略不计，所以讨论它们相位是没必要的。这次课程设计运用了信号与系统、C语言、MATLAB语言的知识，使我加深了对DFT的物理意义的理解，熟练了C和MATLAB的操作。最重要的是锻炼了发现问题和解决问题能力：课程设计要求学生就某些给定的问题进行研究，而这一给定的问题需先分解为几个更小的问题，而在解决每一个小问题的过程中，又会发

38、现一些新的问题，（这些新发现的细节问题有时需要特别注意，因为可以积累一些处理问题的经验），联系前后解决（或至少给出自己的解释）新发现的问题。做课程设计和平时的课本学习的一个很大的差别可能在于：平时的学习是提出一个问题，然后又反反复复（可能是多种不同的解法）的告诉你怎么解决这个问题，主要强调的是被动接受；课程设计则是要求自己独立去解决问题，而且一定程度上说很多问题还是自己发现并提出来的，强调的更多的是主动汲取，是一种拓展学习。可能平时的学习也可借鉴课程设计的经验罢，打破被动接受的常规，主动去思考研究问题。自我评价：比较满意。虽然觉得课程设计完成的没有太大的亮点，私下也以为这个题目较其他几个要简单

39、些（这也是我选择这个题目的原因，另一个原因则是对算法研究和编程比较感兴趣），但至少是自己花了不少时间独立完成的，过程中也学到了不少，所以很高兴看到自己的“作品”。参考文献：1、【美】卡门（Kamen,E.W.）等.信号与系统基础.北京：科学出版社，20022、郑平安.程序设计基础（C语言版）第二版.北京：清华大学出版社，20063、【美】A.V.奥本海姆，R.W.谢弗.离散时间信号处理（第二版）.西安：西安交通大学出版社，20014、【加】Simon Haykin，【美】Barry Van Veen.信号与系统（第二版）.北京：电子工业出版社，20045、丛玉良，王宏志.数字信号处理原理及其MATLAB实现.北京：电子工业出版社，20056、百度文库