快速傅里叶变换FFT的matlab实现和FFT的简单应用

信号与系统课程设计1快速傅里叶变换FFT的matlab实现和FFT的简单应用（1周江280130502908级英才实验三班）（2石磊280130102108级英才实验三班）【摘要】在信号处理中，DFT（离散傅里叶变换）的计算具有举足轻重的地位。但是基于其复杂的计算，直接应用起来十分麻烦，基于此，本文利用Matlab软件对有限长度信号的DFT进行改进，提出FFT（快速傅里叶变换），并利用FFT对所给连续时间和离散时间信号做了频谱分析。关键词：DFT，FFT，有限长度信号，频谱分析。一、前言：傅里叶变换在信号处理中具有十分重要的作用，但是基于离散时间的傅里叶变换具有很大的时间复杂度，根据傅里叶变换理论，对一个有限长度且长度为N的离散信号，做傅里叶变换的时间复杂度为2()ON，当N很大时，其实现的时间是相当惊人的（比如当N为410时，其完成时间为810（为计算机的时钟周期）），故其实现难度是相当大的，同时也严重制约了DFT在信号分析中的应用，故需要提出一种快速的且有效的算法来实现。正是鉴于DFT极其复杂的时间复杂度，1965年..JWCooley和..JWTukey巧妙地利用NW因子的周期性和对称性，提出了一个DFT的快速算法，即快速傅里叶变换（FFT），从而使得DFT在信号处理中才得到真正的广泛应用。本文基于时间抽选奇偶分解，利用Matlab软件实现快速傅里叶变换。基于所编的FFT源程序应用的一个实例，本文对有限长度离散时间和连续时间信号进行频谱分析。二、FFT的具体实现、2.1DFT的算法和时间复杂度对于一个长度为N的离散信号序列[]xn，其DFT变换为10()[]NnkNnXkxnW（1）其中2jnknkNNWe。对任意01mN，101(1)0()[][0][1]...[1]NnmmmNmNNNNnXmxnWxWxWxNW（2）信号与系统课程设计2完成（2）式的计算要做N次复数乘法和1N次复数加法，那么对于（1）就要做2N次复数乘法和(1)NN次复数加法，其时间复杂度是2()ON，计算工作量和运算时间是非常巨大的。故要对其进行性质分析，进而做出快速的算法实现DFT，以便计算机能够快速实现，从而在实际信号处理工作中得到广泛的应用。2.2FFT的提出和Matlab实现2.2.1基于时间抽选的FFT变换原理设序列[]xn的长度为2MN（如果不是，则只需要添加多个0值项进行补充），这种基于序列长度为2M的FFT变换称为‘基2FFT’。本文实现这种基于时间抽选原理的基2FFT。首先将()xn按奇偶分成两个子序列，即：[][2][][21]anxnbnxn，0,1,2...,12Nn（3）令(),()AkBk分别表示[],[]anbn的/2N各点的傅里叶变换，即：/21/212200()()[][]NNnknkNNnnAkBkanWbnW=/21/21/2/200[][]NNknknNNnnanWbnW（4）所以，经过分析可得到：()()()kNXkAkWBk，0,1,...,/21kN(/2)()()kNXNkAkWBk，0,1,...,/21kN（5）利用（5）式计算()Xk，只需要2/2/2NN次乘法。然后由于122MN任然是偶数，故又可以将[]an和[]bn分解成奇偶部分，采用上面同样的算法，重复进行下去，直到信号只剩下一个值时，便得到了信号的傅里叶变换。2.2.2以上算法便是信号长度为2M，基于时间抽选的快速傅里叶变换FFT。基于此算法，下面利用Matlab进行源程序的编译，下面便是FFT的源代码：程序开始functiony=myfft(xr,n)p=0:n-1;%开始倒位序信号与系统课程设计3nu=log2(n);p1=p;b=zeros(1,n);fort=1:nu;p2=floor(p1/2);b=b*2+(p1-2*p2);p1=p2;end;yr(p+1)=xr(b+1);xr=yr;%倒位序结束t=0:n/2-1;%计算因子w开始(只计算w0到wn/2-1)forv=0:n/2-1;w=exp(-2*i*pi*t/n);end;%计算因子w结束form=1:nu;%计算x(k)开始h=2^(m-1);k=1;while(kn+1)fort=1:h;y=bitshift(k-1,nu-m,nu)+1;%求w的幂次数xch(k)=xr(k)+w(y)*xr(k+h);k=k+1;end;fort=1:h;y=bitshift(k-1-h,nu-m,nu)+1;%求w的幂次数xch(k)=xr(k-h)-xr(k)*w(y);k=k+1;end;信号与系统课程设计4end;xr=xch;end;%计算x(k)结束y=xr%输出变换后的结果程序结束附件清单1和2验证了其正确性。三、利用快速傅里叶变换myfft函数实现频域分析实例：3.1设[]xn是由两个正弦信号即白噪声的叠加，请用傅里叶变换对其作频域分析。[]xn的产生：N=2^8;f1=.1;f2=.2;fs=1;a1=5;a2=3;w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);其傅里叶变换和频域分析，以及仿真图详见附件3。3.2利用快速傅里叶变换FFT对连续信号作谱分析。设()cos(200)sin(100)cos(50)xtttt，利用傅里叶变换分析其频谱结构，选择不同的截取长度T，观察存在的截取效应，并试用加窗的方法减少谱间干扰。选取的参数为：（1）频率400,1/fHzTf；（2）采样信号序列[]()()xnxnTwn，()wn是窗函数，选取两种窗函数：矩形窗函数()()NwnRn和Hamming窗，后者在程序中调用函数w(n)=hamming(N)产生成都为N的Hamming窗函数列向量wn；（3）对()xn作2048点FFT作为x(t)的近似连续频谱()Xjf.其中N为采样点数，NfT，T为截取时间长度，取三种长度：0.04s，0,16s,0.32s。信号与系统课程设计5程序源代码及仿真图见附件4。四、总结4.1对myfft实现快速傅里叶变换的评价本文基于时间抽选奇偶分解算法对DFT进行了改进，提出了快速傅里叶变换FFT，并用Matlab实现了FFT，并用其对所给信号进行了频谱分析。在第一个实例里，利用myfft函数对离散信号[]xn进行了傅里叶变换得到它的频域函数和时域函数，从仿真图可以看出myfft实现的傅里叶变换时成功的。在第二个实例里面，也是利用myfft对信号序列进行傅里叶变换，并对信号利用窗函数减少频谱间的干扰，较准确地分析了所给信号的频谱。4.2其它FFT算法简介本文其实还只是快速傅里叶变换的一个简单实现而已。快速傅里叶变换，顾名思义，它应该还有其它实现的算法。本文用的是基于时间抽取奇偶分解的算法，在常规的快速傅里叶变换实现中，还有其它两种非常重要的算法：基于频率抽选奇偶分解算法和混合基算法。其中混合基算法是运用最广泛的一种FFT算法，并在很多领域取得了广泛的应用。由于笔者水平有限，便不再对上述两种算法进行深入的讨论，留待以后慢慢尝试！4.3FFT的进一步深入应用由于傅里叶变换在信号处理中具有举足轻重的作用，故FFT有着非常广泛的应用。下一步的工作便是利用myfft算法设计一个fir滤波器，多所给参杂了噪声的语音信号进行滤波，由于这个设计需要用到其它数字信号处理的知识，可能要等到下学期学习了数字信号课以后才能实现了。4.4课程设计之后的感悟只能用一个字来形容：无奈，累！本来期末时间都要复习考试，可是还要做这个课程设计，花掉我们很多的时间和精力，最后还搞出来一个不是十分满意的东西。哎，一个字：苦不堪言！不过，也并非只是感觉苦不堪言，也有感觉高兴的地方：学了Matlab了，小有收获，还熬了夜，提前体验了一把工作加班到深夜的滋味，很不爽，但又很爽！五、参考文献：[1]余成波，陶红艳。数字信号处理及MATLAB实现，北京：清华大学出版社，2008[2]（美）EdwardW.Kamen,BonnieS.Heck著，高强译。信号与系统基础教程，北京：电子工业出版社，2007[3]曹弋，赵阳。MATLAB实用教程，北京：电子工业出版社，2007信号与系统课程设计6六、附件清单附件1：functiony=myfft(xr,n)xr=[j*pi*1/8j*pi*2/8j*pi*3/8j*pi*4/8j*pi*5/8j*pi*6/8j*pi*7/8j*pi*8/8];n=8;p=0:1:n-1;nu=log2(n);p1=p;b=zeros(1,n);fort=1:nu;p2=floor(p1/2);b=b*2+(p1-2*p2);p1=p2;end;yr(p+1)=xr(b+1);xr=yr;%µ¹Î»Ðò½áÊøt=0:n/2-1;%¼ÆËãÒò×Ów¿ªÊ¼(Ö»¼ÆËãw0µ½wn/2-1)forv=0:n/2-1;w=exp(-2*i*pi*t/n);end;%¼ÆËãÒò×Ów½áÊøform=1:nu;%¼ÆËãx(k)¿ªÊ¼h=2^(m-1);k=1;while(kn+1)fort=1:h;y=bitshift(k-1,nu-m,nu)+1;%ÇówµÄÃÝ´ÎÊýxch(k)=xr(k)+w(y)*xr(k+h);k=k+1;end;fort=1:h;y=bitshift(k-1-h,nu-m,nu)+1;%ÇówµÄÃÝ´ÎÊýxch(k)=xr(k-h)-xr(k)*w(y);k=k+1;end;end;xr=xch;%¼ÆËãx(k)½áÊøend;y=xr得到结果：信号与系统课程设计7附件2一个验证myfft正确与否的小程序：代码：N=8;n=0:N-1;xn=cos(pi*n/8);Xk=myfft(xn,N);plot(Xk);stem(n,abs(Xk),'.');axis([0,20,0,20]);ylabel('|Xk|');title('8点FFT变换');仿真图：信号与系统课程设计8附件3:源代码：%产生两个正弦加白噪声N=2^8;f1=.1;f2=.2;fs=1;a1=5;a2=3;w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);%应用FFT求频谱subplot(2,2,1);plot(x(1:N/4));title('原始信号');f=-0.5:1/N:0.5-1/N;x=myfft(x);y=ifft(x);subplot(2,2,2);plot(f,fftshift(abs(x)));title('频域信号');subplot(2,2,3);plot(real(x(1:N/4)));信号与系统课程设计9title('时域信号');仿真图：附件4源代码：fs=400;T=1/fs;%采样频率和采样间隔Tp=0.04;N=Tp*fs;%采样点数N1=[N,4*N,8*N];%设定三种截取长度form=1:3n=1:N1(m);xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T);Xk=myfft(xn,4096);fk=[0:4095]/4096/T;subplot(3,2,2*m-1);plot(fk,abs(Xk)/max(abs(Xk)));ifm==1title('矩形窗截取');endend信号与系统课程设计10%hamming窗截取form=1:3n=1:N1(m);wn=hamming(N1(m));xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T).*wn';Xk=myfft(xn,4096);fk=[0:4095]/4096/T;subplot(3,2,2*m);plot(fk,abs