对于FFT和IFFT的算法和频谱分析的研究.docx

1、对于FFT和IFFT的算法和频谱分析的研究实验总成绩： 装 订 线报告份数： 西安邮电大学 通信与信息工程学院 科研训练论文专业班级： 学生姓名： 学号(班内序号): 2013 年 9 月 22日装 订 线 对于FFT和IFFT的算法和频谱分析的研究（The algorithms and spectrum analysis of FFT and IFFT）摘要：快速傅里叶变换(FastFourier Tranformation，FFT)是将一个大点数N的DFT分解为若干小点的DFT的组合。将运算工作量明显降低，从而大大提高离散傅里叶变换(DFT)的计算速度,从而更加适合进行实时运算。因各个科学

2、技术领域广泛的使用了FFT技术它大大推动了信号处理技术的进步，现已成为数字信号处理强有力的工具，本论文将比较全面的叙述各种快速傅里叶变换算法原理、特点，并完成了基于MATLAB的实现。最后通过FFT和IFFT的两个应用升余弦滚降和确定函数的频谱分析来分别验证FFT和IFFT的正确性和优越性。关键词：MATLAB FFT IFFT 升余弦滚降函数 三角函数 三角波函数Abstract：Fast Fourier Transformation is mainly achieved by deviding theDFT which is consisted of N points into a lot

3、 of little size of DFT,thus the workload could be reduced to some extent.It is obvious that the reduce of workload may improve the computing speed,so it is fit to get the results immediately.Because of its wide application in various fields of science , the computing speed is considerably improved.F

4、FT and IFFT promote the advancement of signal processing especially Digital Signal Processing,so it has become a powerful tool in DSP.This article will introduce the theory,the feature of FFT and IFFT,then we will get it achieved based on MATLAB.Finally we will get two applications of FFT and IFFT i

5、ncluding raised cosine roll-off and the spectrum analyzing of two determined function to verify the FFT and IFFT correctness and superiority respectively. Key words: MATLAB ,FFT ,IFFT,Raised cosine function,Trigonometric, Triangular wave function引言1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在计算数学杂志上发表了“机器计算傅立

6、叶级数的一种算法”的文章，这是一篇关于计算DFT的一种快速有效的计算方法的文章。它的思路建立在对DFT运算内在规律的认识之上。这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换（Fast Fourier Transform），简称FFT，1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法26使运算效率进一步提高。FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者

7、说数字系统中应用离散傅立叶变换，可以说是进了一大步。随着科学的进步，FFT算法的重要意义已经远远超过傅里叶分析本身的应用。FFT算法之所以快速，其根本原因在于原始变化矩阵的多余行，此特性也适用于傅里叶变换外的其他一些正交变换，例如，快速沃尔什变换、数论变换等等。在FFT的影响下，人们对于广义的快速正交变换进行了深入研究，使各种快速变换在数字信号处理中占据了重要地位。因此说FFT对数字信号处理技术的发展起了重大推动作用。1.FFT的算法 1.1 FFT算法的基本思想1设离散的有限长时间序列x(n)，0nN-1，则其离散傅立叶变换为：这样，矩阵W 中有许多相同的元素，从而可以简化DFT的运算过程F

8、FT 算法有许多形式，笔者只讨论最基本的时间抽取基-2FFT 算法1.2 算法分析一个N 点长序列，直接用DFT 方法需要复数乘法N次；复数加法N(N-1)次。而由图2 可知，采用FFT 则只需要复数乘法次；复数加法 次。当 时， 这样，运算速度提高了1-2 个数量级图1 为FFT 算法和直接DFT 算法所需运算量与计算点数N 的关系曲线显然，N 越大时，优越性越明显但当N 相当大时，利用单机串行进行FFT 运算同样满足不了实时系统的需要 1 1.3 算法的程序实现思想及分析 首先检验待变换的序列的元素个数是否为2的幂次方个，如果不是的话则将其补零使之成为2的整数幂次方个。然后利用已经编好的位

9、倒序子程序输出位倒序序号，将输入序列不断分组，进行处理后从新分组，直至完成最后的处理即可输出变换后的结果。需要注意的是，这里fft的变换后结果的元素个数可能与原输入序列的个数不一样，因为如果不是2的幂次方个的话输入序列后面是要补零的。（主程序为fft_dit和fft_dif）1.4流程示意图整个FFT 频谱分析与显示过程可用图2 所示的流程图示意. fft按频域抽取算法 fft按时域抽取算法 图22.IFFT的频谱变换基本原理在实验中为了简化算法我们直接利用前面已经编好的fft_dit或fft_dif这两个现成函数来实现，基本原理如下：将要变换的频谱序列先取共轭然后将其送入前面的函数，将变化后

10、的结果再取共轭即实现IFFT的功能。主程序为ifft_my 图33 升余弦滚降3.1升余弦滚降原理要实现无码间干扰基带传输时，系统必须满足奈奎斯特准则即：对于上述公式，我们分3种情况来说明其含义：（1）Ts1/2W情况，Z（f）由频率间隔为1/Ts的X（f）曲线无频率重叠地周期性复制并相加构成的，它还是周期性频谱。在这种情况下，有一特定频谱可满足无码间干扰传输的条件，它就是已获广泛应用的升余弦谱。升余弦滤波器的传递函数表示式为：称为滚降因子，取值为01。在=0时，滤波器的带宽W为1/(2Ts),称为奈奎斯特带宽；=0.5时，滤波器的截止频率W=（1+）/（2Ts）=0.75Rs;=1时，滤波器

11、的截止频率W=Rs 23453.2 升余弦滚降的实现及其与快速傅里叶逆变换的关系 升余弦滚降函数是在基带无码间干扰传输中经常用到的频域函数。其主要特性是升余弦滚降函数经过频域平移叠加后能够成一个在各个频域幅度恒定的频域函数。我们在实现升余弦滚降时可以首先在频域实现=0，=0.5，=0.75和=1四种升余弦滚降函数，然后通过自己编写的ifft_my进行傅里叶逆变换，并将变换后的时域结果反映在图上，分别对比不同是对应的时域上的不同的时域特性。3.3 变换前的时域特性和变换后的频域特性 图44 利用fft和ifft进行具体函数频谱分析的实例4.1 三角函数的频谱分析及其信号的恢复 在实际信号的分析中

12、，三角函数是非常常见和基本的信号，在这里我们就对三角函数进行分析。在分析中我们会碰到要分析的函数的采样点数不是2的整数次幂个，我们要对它进行补零处理。4.2 三角波函数的分析 在分析中因为可能会用到自己编写的T2F函数和F2T函数，但是这两个函数仍然是建立在自己编写的fft_dit和ifft_my的基础上的。只是对输入变量进行了更加完善的处理，一个是补零，另一个是对频域进行了搬移，将原来pi2pi的部分搬到-pi0，因为大家在看频谱时比较习惯频谱是关于零对称的。4.3 分析结果5 结论5.1 关于fft和ifft函数的分析 在fft和ifft的实现过程中，的确能降低运算的次数，但是也正好印证了

13、一个很著名的理论，“时间换空间，空间换时间”。在fft和ifft的算法中，我们降低运算的次数是以占用更多的空间换来的。每次将要变换的的序列进行分组，然后对每个组进行处理，虽然降低了运算次数，但是也增加了运算空间的占用。5.2 关于升余弦滚降函数的结论 通过图形我们可以观察出，对于不同的函数，时域主要部分占用的宽度是一样的。但是随着的增大，时域的起伏是越来越小的。因此我们可以认为，越大，时域起伏越小，因此在非线性系统中传输时的失真越小。但是与此同时带来的缺点是占用的越大，频域占用的带宽是越大的，越利于在限带系统中传输，频率的利用率也就越低。5.3 关于三角函数和三角波函数的频谱分析结论 （1）根

14、据理论分析，三角函数的频谱因该是几个冲激函数的组合，三角波函数的频谱应该是Sa函数平方的形式。而在实际分析中我们通过观察频谱图可以看出分析的结果基本符合理论分析结果。这也侧面印证了我们自己编写的fft_dit和ifft_my的正确性和有效性。 （2）在通过频谱恢复时域时，我们观察到恢复后的三角函数的时域函数有所失真，三角波函数基本无失真。原因可能是在T2F和F2T函数中可能因为补零改变了元素个数，因此引起了恢复后的函数时域出现失真。附录1：1蒋冬初,何 飞.FFT算法的并行处理研究.湖南城市学院学报(自然科学版) .2005-06-30 2Oppenheim A V, Schafer R W.

15、Digit al Signal Proces sing M . N Y: Prent ice-Hall , 1975.3Proak is J G, Manolakis D G. Digital Signal Proces sing-PrinciplesM . Algorithms and applicat ions ( 2nd Edit ion) , N Y:Print ice-Hall, 1995.4崔灵智. M atlab 在数字信号处理课程设计中的应用 J. 山东水利职业学院刊, 2008. 3: 11- 12.5樊昌信 曹丽娜 通信原理（第六版） 国防工业出版社 2006.86卢光跃

16、黄庆东 包志强 数字信号处理处理及应用 人民邮电出版社 2012.67王立宁 乐光新 詹菲 MATLAB与通信仿真 人民邮电出版社 2000.01.01 小组成员：组长：08谷群（03101045）组员：31王敏（03101069） 03陈若寰（03101040） 04边钰涵（03101041）任务分配：资料查阅：王敏程序设计：谷群文本编辑：边钰涵技术校对：陈若寰附录：程序代码1 .自己实现的fft(1)fft_dit.m 按时域抽取快速傅里叶变换function num3=fft_dit(in1,k)if nargin=1 N=length(in1);else N=k;end in=zero

17、s(1,N); in(1:length(in1)=in1; num=zeros(1,N); num3=zeros(1,N); k=log2(N); for n=1:N num(n)=in(weidaoxu(n,N)+1); end for n=1:k for c=1:2(k-n) num1=zeros(1,2n); num2=zeros(1,2n); num1=num(c-1)*2n+1:c*2n); for m=1:2(n-1) num2(m)=num1(m)+num1(m+2(n-1)*exp(-i*(m-1)*2*pi/2n); num2(m+2(n-1)=num1(m)-num1(m+

18、2(n-1)*exp(-i*(m-1)*2*pi/2n); end num(c-1)*2n+1:c*2n)=num2; end endnum3=num;(2)fft_dif.m 按频域抽取快速傅里叶变换function num3=fft_dif(in1,k)if nargin=1 N=length(in1);else N=k;end in=zeros(1,N); in(1:length(in1)=in1; num=in; num3=zeros(1,N); k=log2(N); for n=1:k for c=1:2(n-1) num1=zeros(1,2(k+1-n); num2=zeros(

19、1,2(k+1-n); num1=num(c-1)*2(k+1-n)+1:c*2(k+1-n); for m=1:2(k-n) num2(m)=num1(m)+num1(m+2(k-n) num2(m+2(k-n)=(num1(m)-num1(m+2(k-n)*exp(-i*(m-1)*2*pi/2(k+1-n); end num(c-1)*2(k+1-n)+1:c*2(k+1-n)=num2; end end for n=1:N num3(n)=num(weidaoxu(n,N)+1); end（3）weidaoxu.m 位倒序函数function out=niweidaoxu(in,lon

20、g) l=log2(long); num1=zeros(1,l); %out=0; for i=1:l num1(i)=mod(in,2); out=(in-mod(in,2)/2; end（4）ifft_my.m 快速傅里叶逆变换function out=ifft_my(in) num1=conj(in); num2=fft_dit(num1);out=conj(num2)/length(in);2.升余弦滚降函数的时域频域分析 (1)主函数%通过FFT和IFFT来分析升余弦滚降函数的时域和频域的特点clear all;a=0;N=128;l=2;figure(1)f,out_f=sheng

21、yuxiangunjiang(a,N,l);subplot(221);plot(f,out_f);xlabel(f);ylabel(a=0时的频谱幅度);axis(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(222);plot(f,out_t_);xlabel(t);ylabel(a=0时的时域幅度);axis(-0.25 0.25 -0.1 0.5);a=0.5;f,out_f=shengyuxiangunji

22、ang(a,N,l);subplot(223);plot(f,out_f);xlabel(f);ylabel(a=0.5时的频谱幅度);axis(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(224);plot(f,out_t_);xlabel(t);ylabel(a=0.5时的时域幅度);axis(-0.25 0.25 -0.1 0.5);a=0.75;figure(2)f,out_f=shengyuxiangu

23、njiang(a,N,l);subplot(221);plot(f,out_f);xlabel(f);ylabel(a=0.75时的频谱幅度);axis(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(222);plot(f,out_t_);xlabel(t);ylabel(a=0.75时的时域幅度);axis(-0.25 0.25 -0.1 0.5);a=1;f,out_f=shengyuxiangunjiang(

24、a,N,l);subplot(223);plot(f,out_f);xlabel(f);ylabel(a=1时的频谱幅度);axis(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(224);plot(f,out_t_);xlabel(t);ylabel(a=1时的时域幅度);axis(-0.25 0.25 -0.1 0.5);（2）shengyuxiangunjiang.m 产生升余弦滚降频谱的函数function

25、 f,out_f=shengyuxiangunjiang(a,N,l) f=-l:2*l/N:l-2*l/N; out_f=zeros(1,N); out_f(1:N/(2*l)*(1-a)=0; k=-pi/2:l*pi/(a*N):pi/2-l*pi/(a*N) out_f(N/(2*l)*(1-a)+1:N/(2*l)*(1+a)=(sin(k)+1)/2; out_f(N/(2*l)*(1+a)+1:N/2)=1; for m=1:N/2 out_f(N+1-m)=out_f(m); end (3)rand01.m 产生随机的01序列function s=rand01(p,m,n) %

26、 输入参数: % p:0-1分布中1的概率 % m,n:产生的随机变量样本个数mn % 输出:产生的随机变量样本矢量 x=rand(m,n); s=(sign(x-p+eps)+1)/2; % eps = 2(-52).3 确定函数的频谱分析 （1） pinpufenxi.m 主程序clear alldt=0.001;df=0.2;fs=1/dt;t=-10:dt:10-dt;y=sin(2*pi*100*t)+2*sin(2*pi*400*t);M1,m1,df1,f1=T2F(y,dt,df,fs);figure(1)subplot(211);plot(t,y);xlabel(t);yla

27、bel(三角函数频谱分析前的时域);axis(-0.05,0.05,-5,5);subplot(212);plot(f1,abs(M1);xlabel(f);ylabel(三角函数频谱分析后的频谱);figure(2)z=zeros(1,length(t);z=zeros(1,length(t)/4),0:5*4/length(t):5-5*4/length(t),zeros(1,length(t)/2);for n=1:length(t)/2 z(length(t)+1-n)=z(n);endM2,m2,df2,f2=T2F(z,dt,df,fs);subplot(211);plot(t,z

28、);xlabel(t);ylabel(三角波频谱分析前的时域);axis(-10,10,-1,5);subplot(212);plot(f2,abs(M2);xlabel(f);ylabel(三角波频谱分析后的频谱);axis(-1,1,0,30);figure(3)y1=F2T(M1,fs);yy1=y1(1:length(t);subplot(211);plot(t,-yy1);xlabel(t);ylabel(三角函数恢复后的时域函数);axis(-0.05,0.05,-5,5);y2=F2T(M2,fs);yy2=y2(1:length(t);subplot(212);plot(t,abs(yy2);xlabel(t);ylabel(三角波恢复后的时域函数);axis(-10,10,-1,5);（2） T2F.m 改进后的傅里叶变换 6function M,m,df1,f=T2F(m,ts,