应用FFT实现信号频谱分析

1、电信类课程试验报告学院：理学院系别：电子信息工程课程名称：数字信号处理姓名：05656学号：05698日期：11月28日第：四次实验实验名称：应用FFT实现彳百号频谱分析一、实验目的(1)能够熟练掌握快速离散傅里叶变换(FastFourierTransform,FFT)的原理及应用FFT进行频谱分析的基础方法。(2)对离散傅里叶变换的主要性质及FFT在数字信号处理中的重要作用有进一步的了解。二、主要函数简介1.离散傅里叶变换(DFT)及其主要性质DFT表示离散信号的离散频谱，DFT的主要性质中有奇偶对称特性、虚实特性等。通过实验可以加深理解。实序列的DFT具有偶对称的实部和奇对称的虚部，这可以

2、证明如下：由定义，可得N1Xm=ZxkwNmk田N1N2n-2几=Zxkcos(km)-jZxksin(km)k=eNk=eNN-1XN-m=ZMkwNN5)kkfN-4=ZxkwNkWiNknk=QN-4=ZxkwN'mk=eN"2nM2n=Lxkcos(km)-jLxksin(km)k=eNk=eN.一*所以：xk=xN-k实序列DFT的这个特性，在本实验中可以通过实指数序列及三角序列看出来。对于单一频率的三角序列来说，它的DFT谱线也是单一的，这个物理意义可以从实验中得到验证，在理论上可以推导如下：设：.2二._Xk=sin(k)Rn"N其DFT为：N2二皿X

3、m-xkeNkz0N2二2N：km八sin(k)eNk30N1N1j221k%(eN-e2jkJ22k422JkmN)eN1Nj2JI(m)k(eN2jk日2i(m1)k七N)从而-j2二kN)=01N$nNx(i)=.(1-e)=-j2jkJ2j2X(N-2)=0X(N-1)=!其"*力=_'=jN2Jkf2j2.、.*以上这串式中X0反应了xk的支流分量，X1是xk的一次谐波，又根据虚实特性XN-1=X1而其他分量均为零。2二、当周期减小时显然sin(一rk)RNk的谱只应该在k=3及k=N-3才有分量，实验者可以通过可上述相同的步N3骤加以理论证明。2二由于cos(k)

4、*Rnk及sin(N性质可在本实验中得到验证。nk),RNk相位差一，所以它的DFT只包括实部而没有虚部，以上这些22.利用DFT对信号进行频谱分析DFT的重要应用之一是对时域连续信号的频谱进行分析，称为傅里叶分析，时域连续信号离散傅里叶分析的基本步骤如图51所示。凡R)a|上尸尸</可AID|忑)DFT毋乃工用於)下w(理)图5-1时域连续信号离散傅里叶分析的处理步骤其中消混叠低通滤波器LPF(预滤波器)的引入，是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠的影响。实际工作中，时域离散信号xk的时宽是很长的甚至是无限长的(例如语言或音乐信号)。由于DFT的需要，必须把xk限

5、制在一定的时间间隔之内，即进行数据截断。数据的截断相当于加窗处理。因此，在计算xk的DFT之前，用一个时域有限的窗函数wnk加到xk上是非常必要的。Xc(t)通过A/D变换器转成取样序列xk进行加窗处理，即vk=xkwNko加窗对频域的影响，用周期卷积表示。卜11J1Jjn2二X(ej)X(j(n飞)X(j(-)Tn二二Tn-.二二TT,ji1其中X(ej§或X(eT)为xjt)的频谱。在实际应用中，消混叠低通滤波器的阻带不可能式无限衰减的，故由Xc(j8)周期延拓得到的X(ejC)由非零混叠，即出现混叠现象。由于进行dft的需要，必须对序列xk进行加窗处理，即vk=xk'W

6、Nk,加窗对频域的影响，用周期卷积表示。V(ej')=1X(e»w(ej(W<)du2二,一奠最后是进行DFT应算。加窗后的DFT为N二_j纽mkVm=二vkeN0-m-N-1k=0其中假设窗函数长L小于或等于DFT长度N。有限长序列vk=xkwNk的dft相当于vk傅里叶变换的等间隔取样。Vm=V(ejVm便是Sc(t)的离散频率函数。因为DFT频率间隔为2L,且模拟频率q和数字频率建间的关系为T=oT,所以离散频率点对应的模拟N频率为2"co=NT显然频率分辨率4f为NT2二2二利用DFT计算频谱，只名出频谱Cm=m或0m=m的频率分量，即频率的取样值，而

7、不可能mNmNT得到连续的频谱函数。如果在两个离散的谱线之间有一个特别大的频谱分量，就无法检测出来了。为了在保持原来频谱形状不变的情况下，使谱线加密，即使频域取样点数增加，从而使原来看不到的频谱分量变得可以看到，可以通过在信号数据的末端补加一些零值点，使DFT计算周期内点数增加，但又不改变原有的记录数据的方法来实现。3.快速离散傅里叶变换(FFT)快速离散傅里叶变换是计算离散傅立叶变换的一种快速算法，为了提高应算速度，FFT将DFT的计算逐次分解成较小点数的DFT。按时间抽取(Decimation-In-Time,DIT)FFT算法把输入序列xk按起k值为偶数或是奇数分解成越来越短的序列。按频

8、域抽取(Decimation-In-Time,DIT)FFT算法是把输出序列xm按其m值是偶数或是奇数来分解成越来越短的序列。具体推导过程及原理可参见数字信号处理教科书。三、实验程序指数序列程序2.1clearallN=100;n=0:N-1;xn=0.9.An;XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1)plot(n,xn)xlabel('n');ylabel('x(n)');title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-

9、1;k=k*(2/100)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=100');2.2clearallN=100;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=100')subplot(1

10、,2,2)k=0:length(magXK)-1;title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/100)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=100');2.3 复合函数clearallN=100;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK=fft(xn,N);magXK=abs(XK);phaXK=a

11、ngle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/100)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=100

12、');3.1clearallN=32;n=0:N-1;xn=0.9.An;XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n,);ylabel(,x(n)');title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/32)stem(k,magXK,'.');x

13、label('k');ylabel('|X(k)|');title('X(k)N=32');3.2N=16iBiarlTQInklpTibjtlii>141-3.3N=644.1 N=N+2=34clearallN=34;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=34')

14、subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=34')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/34)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=34');4.2 N=32clearallN=32;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2

15、,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/32)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=32');4.3 N=64clearallN

16、=64;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=64')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=64')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/64)stem(k,magXK,'.');xlabe

17、l('k');ylabel('|X(k)|');title('X(k)N=64');5.1 N=32clearallN=32;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=32')subplot(1,2,2)k=0:length(m

18、agXK)-1;title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/32)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=32');5.2 N=64clearallN=64;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=64')subplot(1,2,2)k=0:lengt