数字信号处理上机实验指导书

1、.数字信号处理上机实验指导书电子与信息工程学院目录实验一离散信号产生及频谱的绘制4实验二时域采样与频域采样9实验三离散傅立叶变换及谱分析16实验四用FFT对信号作频谱分析21实验五双线性变换法设计IIR数字滤波器32实验六 FIR数字滤波器设计与软件实现38附录：实验用MATLAB语言工具箱函数简介45引言“数字信号处理”是一门理论和实验密切结合的课程，为了深入地掌握课程内容，应当在学习理论的同时，做习题和上机实验。上机实验不仅可以帮助学生深入地理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。所以，根据本课程的重点要求编写了四个实验。第一章是全书的基础内容，抽样定理、时域离散系统的时

2、域和频域分析以及系统对输入信号的响应是重要的基本内容。由于第一章大部分内容已经在前期信号与系统课程中学习完，所以可通过实验一帮助学生温习以上重要内容，加深学生对“数字信号处理是通过对输入信号的一种运算达到处理目的” 这一重要概念的理解。这样便可以使学生从信号与系统课程顺利的过渡到本课程的学习上来。第二章、三章DFT、FFT是数字信号处理的重要数学工具，它有广泛的使用内容。采用实验二、实验三加深理解DFT的基本概念、基本性质。FFT是它的快速算法，必须学会使用。数字滤波器的基本理论和设计方法是数字信号处理技术的重要内容。学习这一部分时，应重点掌握IIR和FIR两种不同的数字滤波器的基本设计方法。

3、IIR滤波器的单位冲激响应是无限长的，设计方法是先设计模拟滤波器，然后再通过SZ平面转换，求出相应的数字滤波器的系统函数。这里的平面转换有两种方法，即冲激响应不变法和双线性变换法，后者没有频率混叠的缺点，且转换简单，是一种普遍应用的方法。FIR滤波器的单位冲激响应是有限长的，设计滤波器的目的即是求出符合要求的单位冲激响应。窗函数法是一种基本的，也是一种重要的设计方法。学习完第七章后可以进行实验四。由于数字信号处理实验的主要目的是验证数字信号处理的有关理论，进一步理解巩固所学理论知识。所以，实现实验用的算法语言可以有许多种，但为了提高实验效率，要求学生用编程效率比C语言高好几倍的MATLAB语言

4、。下面介绍MATLAB的主要特点。（有关MATLAB的启动、程序运行和有关信号处理工具箱函数等内容将放到最后附录中介绍。）MATLAB是一种交互式的以矩阵为基本数据结构的系统。在生成矩阵对象时，不要求明确的维数说明。所谓交互式，是指MATLAB的草稿纸编程环境。即用户每输入一条命令并按回车键，MATLAB系统便解释执行之，并显示执行结果。根据该结果，用户立即知道刚输入的命令的正确性，或利用中间结果进行其他处理等。与C语言或FORTRON语言做科学数值计算的程序设计相比较，利用MATLAB可节省大量的编程时间。将其用于数字信号处理实验，则可大大提高实验效率，在有限的上机时间内，实验内容可增加几倍

5、。例如，C语言FFT子程序有70多行，而用MATLAB只调用一个fft函数即可实现对序列进行FFT计算。另外，MATLAB的工具箱及图形显示（打印）功能，可满足各层次人员直观、方便的进行分析、计算和设计工作，从而可大大节省时间。例如，序列的卷积、滤波，系统函数H(z)的幅频特性和相频特性等计算，均有现成的工具箱函数。而用其它算法语言完成这些计算的编程比较麻烦，且程序较长。由于上述特点，在美国一些大学里，MATLAB已成为辅助教学的有益工具。MATLAB已成功地用于数字信号处理课程中的问题分析、实验、滤波器设计及计算机模拟。附录中所介绍的信号处理工具箱函数及绘图函数基本可满足本教材所要求的上机实

6、验需要。对序列进行谱分析的MATLAB程序及运行结果见附录。实验一离散信号产生及频谱的绘制一、实验目的: （1）熟悉Matlab环境。（2）掌握 Matlab 中一些基本函数的建立方法（2）通过编程绘制的幅度相位谱加深理解系统的特性二、实验内容1、编写程序，产生以下离散序列：（1）f(n)=(n) (-3<n<4)n1=-3;n2=4;n0=0;n=n1:n2;x=n=n0;stem(n,x,'filled');axis(n1,n2,0,1.1*max(x);xlabel('时间(n)');ylabel('幅度x(n)');title

7、('单位脉冲序列');（2）f(n)=u(n) (-5<n<5)n1=-5;n2=5;n0=0;n=n1:n2;x=n>=n0;stem(n,x,'filled');axis(n1,n2,0,1.1*max(x);xlabel('时间(n)');ylabel('幅度x(n)');title('单位阶跃序列');box（3） (0<n<16)n1=16;a=-0.1;w=1.6*pi;n=0:n1;x=exp(a+j*w)*n);subplot(2,2,1);plot(n,real(x)

8、;title('复指数信号的实部');subplot(2,2,3);stem(n,real(x),'filled');title('复指数序列的实部');subplot(2,2,2);plot(n,imag(x);title('复指数信号的虚部');subplot(2,2,4);stem(n,imag(x),'filled');title('复指数序列的虚部');box %box on 加边框/ box off不加边框/ box 加或不加切换，注意只对上面一个图有效（4）f(n)=3sin(n/4)

9、 (0<n<20)f= 1/8;Um=3;nt=2;%显示的周期数目N=32; T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=Um*sin(2*f*pi*tn);subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(t)');subplot(2,1,2);stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');box on %对当前图形右边及上边加边框（5）一个连续的周期性方波信号

10、频率为200Hz，信号幅度在-1+1V之间，要求在图形窗口上显示其两个周期的波形。以4kHz的频率对连续信号进行采样，编写程序生成连续信号和其采样获得的离散信号波形。f=200;nt=2; %显示周期数N=20;T=1/f;dt=T/N; %每个周期显示20个离散值，4kHz的频率n=0:nt*N-1;tn=n*dt;x=square(2*f*pi*tn,25); %其中25为占空比subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(t)');subplot(2,1,2);stem(tn,x

11、);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');box（6）绘制的幅度谱和相位谱。clc;%本语句的作用是清除命令执行界面中所有的输出信息clear ;%清除workspace中所有的变量clf;%清除所有的绘图内容(如果本次程序执行前已经有绘图窗口存在，则可能将本程序将要绘制的图形绘制到之前的窗口中，可能导致疑惑)w=0:0.1:2*pi;%将连续w分成极小的间隔，%(1)取消下面一行可以做第1个幅度相位谱hjw=1+exp(-j*w)+exp(-2*j*w)+exp(-3*j*w)+exp(-4*j*w)+exp(-5*

12、j*w)+exp(-6*j*w);mag=abs(hjw);%取模/幅度值ang=angle(hjw);%相位figure(1);%第一个窗口subplot(211)%将第一个窗口分成2行1列的子图，分的方法是从左到右，从上到下plot(w,mag);%以w为横坐标，幅度谱为纵坐标绘制幅度谱set(gca,'ytick',-0.5:0.7:3);%设定y轴显示范围，从-0.5到3，每一个刻度间隔是0.4set(gca,'Xtick',0:pi/8:2*pi);%同y轴含义，此处对X轴标注set(gca,'XtickLabel','0

13、9;,'pi/8','pi/4','3*pi/8','pi/2','5*pi/8','3*pi/4','7*pi/8','pi','9*pi/8','10*pi/8','11*pi/8','3*pi/2','13*pi/8','7*pi/4','15*pi/8','2*pi');%设定横轴上坐标如何显示，w1=0:2*pi/7:2*pi;hj

14、w1=1+exp(-j*w1)+exp(-2*j*w1)+exp(-3*j*w1)+exp(-4*j*w1)+exp(-5*j*w1)+exp(-6*j*w1);mag1=abs(hjw1);hold on;stem(w1,mag1,'r');title('幅度谱');%设定第1幅子图的标题subplot(212);%第2个子图plot(w,ang);%显示相位%stem(w,mag)set(gca,'ytick',-pi:pi/4:pi);set(gca,'Xtick',0:pi/8:2*pi);set(gca,'Xti

15、ckLabel','0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','3*pi/4','7*pi/8','pi','9*pi/8','10*pi/8','11*pi/8','3*pi/2','13*pi/8','7*pi/4','15*pi/8','2*pi');title('

16、;相位谱');%设定第2幅子图的标题grid on %可以在同一个子图里面叠绘多张图（7）已知系统函数，用MATLAB绘出8阶系统函数的零极点图、幅频响应和相频响应曲线。b=1 0 0 0 0 0 0 0 -1; %H(z)的分子多项式系数矢量a=1; %H(z)的分母多项式系数矢量subplot(1,3,1);, zplane(b,a); %绘制H(z)的零极点图H,w=freqz(b,a); %计算系统的频率响应subplot(1,3,2); plot(w/pi,abs(H); %绘制幅频响应曲线axis(0,1,0,2.5);xlabel('omega/pi');

17、ylabel('|H(ejomega)|'); subplot(1,3,3); plot(w/pi,angle(H); %绘制相频响应曲线xlabel('omega/pi');ylabel('phi(omega)');实验二时域采样与频域采样一、实验目的（1）掌握模拟信号采样前后频谱的变化（2）要求掌握频率域采样会引起时域周期化的概念（3）频率域采样定理及其对频域采样点数选择的指导作用。二、实验原理与方法时域采样定理的要点是：a) 对模拟信号以间隔T进行时域等间隔理想采样，形成的采样信号的频谱是原模拟信号频谱以采样角频率（）为周期进行周期延拓。公

18、式为：b) 采样频率必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的频谱不产生频谱混叠。利用计算机计算上式并不方便，下面我们导出另外一个公式，以便用计算机上进行实验。理想采样信号和模拟信号之间的关系为：对上式进行傅立叶变换，得到：在上式的积分号内只有当时，才有非零值，因此：上式中，在数值上，再将代入，得到：上式的右边就是序列的傅立叶变换，即上式说明理想采样信号的傅立叶变换可用相应的采样序列的傅立叶变换得到，只要将自变量用代替即可。频域采样定理的要点是：a) 对信号x(n)的频谱函数X(ej)在0，2上等间隔采样N点，得到则N点IDFT得到的序列就是原序列x(n)以N为周期进行周期延拓后

19、的主值区序列，公式为：b) 由上式可知，频域采样点数N必须大于等于时域离散信号的长度M(即NM)，才能使时域不产生混叠，则N点IDFT得到的序列就是原序列x(n),即=x(n)。如果N>M，比原序列尾部多N-M个零点；如果N<M，z则=IDFT发生了时域混叠失真，而且的长度N也比x(n)的长度M短，因此。与x(n)不相同。在数字信号处理的应用中，只要涉及时域或者频域采样，都必须服从这两个采样理论的要点。 对比上面叙述的时域采样原理和频域采样原理，得到一个有用的结论，这两个采样理论具有对偶性：“时域采样频谱周期延拓，频域采样时域信号周期延拓”。因此放在一起进行实验。二、实验内容及步骤

20、（1）时域采样理论的验证。给定模拟信号，式中A=444.128，=50，=50rad/s，它的幅频特性曲线如图3.1图2.1 的幅频特性曲线现用DFT(FFT)求该模拟信号的幅频特性，以验证时域采样理论。按照的幅频特性曲线，选取三种采样频率，即=1kHz，300Hz，200Hz。观测时间选。为使用DFT，首先用下面公式产生时域离散信号，对三种采样频率，采样序列按顺序用，表示。因为采样频率不同，得到的，的长度不同， 长度（点数）用公式计算。选FFT的变换点数为M=64，序列长度不够64的尾部加零。X(k)=FFTx(n) ， k=0,1,2,3,-,M-1式中k代表的频率为 。要求： 编写实验程

21、序，计算、和的幅度特性，并绘图显示。观察分析频谱混叠失真。（2）频域采样理论的验证。给定信号如下：编写程序分别对频谱函数在区间上等间隔采样32和16点，得到：再分别对进行32点和16点IFFT，得到：分别画出、的幅度谱，并绘图显示x(n)、的波形，进行对比和分析，验证总结频域采样理论。提示：频域采样用以下方法容易变程序实现。 直接调用MATLAB函数fft计算就得到在的32点频率域采样 抽取的偶数点即可得到在的16点频率域采样，即。 当然也可以按照频域采样理论，先将信号x(n)以16为周期进行周期延拓，取其主值区（16点），再对其进行16点DFT(FFT),得到的就是在的16点频率域采样。四、

22、思考题： 如果序列x(n)的长度为M，希望得到其频谱在上的N点等间隔采样，当N<M时， 如何用一次最少点数的DFT得到该频谱采样？五、实验报告及要求a）运行程序打印要求显示的图形，。b) 分析比较实验结果，简述由实验得到的主要结论c) 简要回答思考题d) 附上程序清单和有关曲线。六、实验程序清单1 时域采样理论的验证程序清单% 时域采样理论验证程序exp2a.mTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)% Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;alph=pi*50*20.5;omega

23、=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn='xa(nT)'subplot(3,2,1);tstem(xnt,yn);%调用自编绘图函数tstem绘制序列图box on;title('(a) Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk);title('(a) T*FTxa(nT),Fs=1000Hz');xlabel('f(Hz)');ylabel

24、('幅度');axis(0,Fs,0,1.2*max(abs(Xk)%=% Fs=300Hz和 Fs=200Hz的程序与上面Fs=1000Hz完全相同。2 频域采样理论的验证程序清单%频域采样理论验证程序exp2b.mM=27;N=32;n=0:M;%产生M长三角波序列x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb;Xk=fft(xn,1024);%1024点FFTx(n), 用于近似序列x(n)的TFX32k=fft(xn,32);%32点FFTx(n)x32n=ifft(X32k);%32点IFFTX32(k)得到x3

25、2(n)X16k=X32k(1:2:N);%隔点抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16点IFFTX16(k)得到x16(n)subplot(3,2,2);stem(n,xn,'.');box ontitle('(b) 三角波序列x(n)');xlabel('n');ylabel('x(n)');axis(0,32,0,20)k=0:1023;wk=2*k/1024;%subplot(3,2,1);plot(wk,abs(Xk);title('(a)FTx(n)');xlabel(

26、'omega/pi');ylabel('|X(ejomega)|');axis(0,1,0,200)k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box ontitle('(c) 16点频域采样');xlabel('k');ylabel('|X_1_6(k)|');axis(0,8,0,200)n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box ontitle('(d) 16点IDFT

27、X_1_6(k)');xlabel('n');ylabel('x_1_6(n)');axis(0,32,0,20)k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box ontitle('(e) 32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis(0,16,0,200)n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box ontitle('(f)

28、 32点IDFTX_3_2(k)');xlabel('n');ylabel('x_3_2(n)');axis(0,32,0,20)%=tstem 程序清单=function tstem(xn,yn)%时域序列绘图函数% xn:信号数据序列，yn:绘图信号的纵坐标名称（字符串）n=0:length(xn)-1;stem(n,xn,'.');box onxlabel('n');ylabel(yn);axis(0,n(end),min(xn),1.2*max(xn)七、实验程序运行结果1 时域采样理论的验证程序运行结果exp2a

29、.m如图3.2所示。由图可见，采样序列的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时，在折叠频率110Hz附近频谱混叠更很严重。图2.22 时域采样理论的验证程序exp2b.m运行结果如图3.3所示。图2.3该图验证了频域采样理论和频域采样定理。对信号x(n)的频谱函数X(ej)在0，2上等间隔采样N=16时， N点IDFT得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列：由于N<M，所以发生了时域混叠失真，因此。与x(n)不相同，如

30、图图3.3(c)和(d)所示。当N=32时，如图图3.3(c)和(d)所示，由于N>M，频域采样定理，所以不存在时域混叠失真，因此与x(n)相同。实验三离散傅立叶变换及谱分析一、 实验目的1掌握离散傅里叶变换的计算机实现方法。2检验实序列傅里叶变换的性质。3掌握计算序列的圆周卷积的方法。4熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。5学习用DFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT。二、 实验内容1实现离散傅里叶变换。2计算序列圆周卷积。3计算实序列傅里叶变换并检验DFT性质。4实现连续信号傅里叶变换以及由不

31、同采样频率采样得到的离散信号的傅里叶变换。5实现补零序列的傅里叶变换。6实现高密度谱和高分辨率谱，并比较二者的不同。三、 实验报告要求见各程序要求%以下为4个扩展函数% （1）离散傅立叶变换 采用矩阵相乘的方法function Xk=dft(xn,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.(nk);Xk=xn*WNnk;%（2）逆离散傅立叶变换function xn=idft(Xk,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.(-nk);xn

32、=(Xk*WNnk)/N;% (3) 实序列的分解 % 实序列可分解为共扼对称分量% 和共扼反对称分量function xec,xoc=circevod(x)N=length(x);n=0:(N-1);xec=0.5*(x+x(mod(-n,N)+1); %根据上面的公式计算，mod函数为取余xoc=0.5*(x-x(mod(-n,N)+1);% (4) 序列的循环移位 function y=cirshftt(x,m,N)if length(x)>N error('N mustbe >= the length of x') %要求移位周期大于信号长度endx=x z

33、eros(1,N-length(x);n=0:1:N-1;n=mod(n-m,N);y=x(n+1);%例1 本例检验实序列的性质DFTxec(n)=ReX(k)DFTxoc(n)=ImX(k)%设 x(n)=10*(0.8).n 0<=n<=10 将x(n)分解为共扼对称及共扼反对称部分%实验报告要求：（1）将实验结果画出 （2）实验结果说明什么n=0:10;x=10*(0.8).n;xec,xoc=circevod(x);subplot(2,1,1);stem(n,xec); %画出序列的共扼对称分量title('Circular -even component'

34、;)xlabel('n');ylabel('xec(n)');axis(-0.5,10.5,-1,11)subplot(2,1,2);stem(n,xoc); %画出序列的共扼反对称分量title('Circular -odd component') xlabel('n');ylabel('xoc(n)');axis(-0.5,10.5,-4,4)figure(2)X=dft(x,11); %求出序列的DFTXec=dft(xec,11); %求序列的共扼对称分量的DFTXoc=dft(xoc,11); %求序列的

35、共扼反对称分量的DFTsubplot(2,2,1);stem(n,real(X);axis(-0.5,10.5,-5,50)title('RealDFTx(n)');xlabel('k'); %画出序列DFT的实部subplot(2,2,2);stem(n,imag(X);axis(-0.5,10.5,-20,20)title('ImagDFTx(n)');xlabel('k'); %画出序列DFT的虚部subplot(2,2,3);stem(n,Xec);axis(-0.5,10.5,-5,50)title('DFTxe

36、c(n)');xlabel('k');subplot(2,2,4);stem(n,imag(Xoc);axis(-0.5,10.5,-20,20)title('DFTxoc(n)');xlabel('k');% 例2 本例为计算序列的圆周卷积程序% 运行之前应在命令窗口输入 x1,x2,N 的值%实验报告要求：自己选择2个序列进行计算，将实验结果写出。if length(x1)>Nerror('N mustbe >= the length of x1')endif length(x2)>Nerror(&#

37、39;N mustbe >= the length of x2')endx1=x1 zeros(1,N-length(x1); %将x1,x2补0成为N长序列x2=x2 zeros(1,N-length(x2);m=0:1:N-1;x2=x2(mod(-m,N)+1); %该语句的功能是将序列翻褶，延拓，取主值序列H=zeros(N,N);for n=1:1:N %该for循环的功能是得到x2序列的循环移位矩阵 H(n,:)=cirshftt(x2,n-1,N); %和我们手工计算圆周卷积得到的表是一致的endy=x1*H' %用矩阵相乘的方法得到结果% 例3 本例验证采

38、样定理%令，绘制其傅立叶变换。用不同频率对其进行采样，分别画出离散时间傅立叶变换。已给出采样频率为时的的程序 %实验报告要求：（1）请写出时的程序 (2)将实验结果画出 (3)实验结果说明什么（采样频率不同结果有何不同）。Dt=0.00005; %步长为0.00005st=-0.005:Dt:0.005;xa=exp(-1000*abs(t); %取时间从-0.005s到0.005s这段模拟信号Wmax=2*pi*2000; %信号最高频率为2*2000K=500; %频域正半轴取500个点进行计算k=0:1:K;W=k*Wmax/K; % 求模拟角频率Xa=xa*exp(-j*t'*

39、W)*Dt; %计算连续时间傅立叶变换（利用矩阵运算实现） Xa=real(Xa); %取实部W=-fliplr(W),W(2:501); %将角频率范围扩展为从-到+Xa=fliplr(Xa),Xa(2:501);%A = 1 3 5 7 9 则fliplr(A)= 9 7 5 3 1 subplot(2,2,1);plot(t*1000,xa); %画出模拟信号，横坐标为时间（毫秒），纵坐标为幅度xlabel('time(millisecond)');ylabel('xa(t)');title('anolog signal');su

40、bplot(2,2,2);plot(W/(2*pi*1000),Xa*1000);%画出连续时间傅立叶变换 xlabel('frequency(kHZ)'); %横坐标为频率（kHz）ylabel('xa(jw)'); %纵坐标为幅度title('FT');%下面为采样频率5kHz时的程序Ts=0.0002; %采样间隔为n=-25:1:25;x=exp(-1000*abs(n*Ts); %离散时间信号K=500;k=0:1:K;w=pi*k/K; %w为数字频率X=x*exp(-j*n'*w); %计算离散时间傅立叶变换（序列的傅立叶变

41、换）X=real(X);w=-fliplr(w),w(2:K+1);X=fliplr(X),X(2:K+1);subplot(2,2,3);stem(n*Ts*1000,x); %画出采样信号（离散时间信号）xlabel('time(millisecond)');gtext('Ts=0.2ms'); %该语句可以将引号中的内容放置在figure中的任何地方，只需 %将十字的中心放在想放置内容的地方，然后按鼠标即可。ylabel('x1(n)');title('discrete signal');subplot(2,2,4);plo

42、t(w/pi,X); %画出离散时间傅立叶变换xlabel('frequency(radian)'); %横坐标为弧度ylabel('x1(jw)');title('DTFT');%例4 本例说明补零序列的离散傅立叶变换%序列，已给出序列的傅立叶变换程序和将原序列补零到10长序列的DFT%实验报告要求： （1）编写将序列补零到20长后计算DFT的程序(2)给出实验结果(3)实验结果说明什么（即序列补零后进行DFT,频谱有何变化）n=0:4;x=ones(1,5); %产生矩形序列k=0:999;w=(pi/500)*k;X=x*(exp(-j*p

43、i/500).(n'*k); %计算离散时间傅立叶变换Xe=abs(X); %取模subplot(3,2,1);stem(n,x);ylabel('x(n)'); %画出矩形序列subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %画出离散时间傅立叶变换N=10;x=ones(1,5),zeros(1,N-5); %将原序列补零为10长序列n=0:1:N-1;X=dft(x,N); %进行DFTmagX=abs(X);k=(0:length(magX)'-1)*N/length(magX);subpl

44、ot(3,2,3);stem(n,x);ylabel('x(n)'); %画出补零序列subplot(3,2,4);stem(k,magX); %画出DFT结果axis(0,10,0,5);ylabel('|X(k)|');%例5 本题说明高密度谱和高分辨率谱之间的区别，高密度谱是信号补零后得到的，虽然谱线相当密，但是因为信号有效长度不变，所以其分辨率也不变，因此还是很难看出信号的频谱成分。高分辨率谱是将信号有效长度加长，因此分辨率提高，可以看出信号的成分。%有一个序列为 （该序列周期计算可得40）%（1）下面给出有10个有效采样点序列的DFT程序%（2）请写出

45、将第一问中的10长序列补零到40长，计算其DFT%（3）采样n=0：39，计算有40个有效采样点的序列的DFT%实验报告要求： （1）请编写将有10个有效采样点的序列补零到40长后计算DFT的程序 (2) 请编写计算有40个有效采样点的序列的DFT的程序 (3) 将实验结果画出并分析实验结果说明什么M=10;n=0:M-1;x=2*cos(0.35*pi*n)+cos(0.5*pi*n);subplot(2,1,1);stem(n,x);title('没有足够采样点的信号');Y=dft(x,M);k1=0:M-1;w1=2*pi/M*k1;subplot(2,1,2);ste

46、m(w1/pi,abs(Y);title('信号的频谱');M=10;n=0:M-1;b=40;s=0:b-1;x=2*cos(0.35*pi*n)+cos(0.5*pi*n), zeros(1,30);subplot(2,1,1);stem(s,x);title('有足够采样点的信号');Y=dft(x,b);k1=0:b-1;w1=2*pi/b*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title('信号的频谱');实验四用FFT对信号作频谱分析一、实验目的1、理解用FFT对周期序列进行频谱分析时所面临的问题并掌

47、握其解决方法。2、掌握用时域窗函数加权处理的技术。3、理解用FFT对非周期信号进行频谱分析所面临的问题并掌握其解决方法。二、实验原理与计算方法、对周期序列进行频谱分析应注意的问题k k(a)时域周期整数倍截断 (b)时域非周期整数倍截断图 4-1 周期函数的幅频曲线k图 4-2 矩形窗的频谱对时间序列作FFT时，实际上要作周期延拓（如果取长序列的一段进行计算还要先作截断）。周期序列是无限长时间序列，如果截断区间刚好就是该序列周期的整数倍，那么在进行周期延拓后，将还原出原来的周期序列，由此可以较精确地计算出的该周期序列的频谱。反之，如果截断区间并不是该序列周期的整数倍，那么在进行周期延拓后，就不

48、可能还原出原来的周期序列，由此计算出的频谱与该周期序列的频谱存在误差，而且误差的大小与截断区间的选取直接相关，如图4-1所示，其中幅度频谱的量值表示为，以dB(分贝)为单位。这种误差是由于周期序列与矩形截断序列相乘在频域产生二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，如图4-2所示。除了k = 0处主瓣内集中了大部分的能量外，两旁的较小峰值处的旁瓣也分散了一部分能量，它与周期序列频谱卷积的结果使原来集中的频谱展宽，称为频率泄漏。如果对已知周期的信号作频谱分析，在进行时域截断时，完全可以选取其周期的整数倍裁取，从而可以避免这种频率泄漏的发生。不过，通常需要进行频谱分析的信号是周期未知的信号

49、，或随机信号，无法判断它的周期值，为了尽量避免频率泄漏对结果的影响，在作时间截断时，就应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。2、时域窗函数的应用作为截断函数，矩形窗在作时间截断时，对所截取区间内的信号不加以任何影响，而其它的窗函数都将对所截取区间内的信号作加权处理。除三角窗、Hanning窗和Hamming窗外，常用的窗函数还有很多，例如Parzen窗、Kaiser窗、Chebyshev窗、Tukey窗、Poisson窗、Caushy窗、Gaussian窗和Blackman窗等等。本次实验利用窗函数作时域加权截断。 0 t 0 k(a) 正弦函数的加权的非周期时域截断(b)减小了泄

50、漏的频谱图 4-3 采用Hanning窗加权后的时域截断和频谱图 4-3 中给出了采用Hanning窗对正弦函数作非整周期的时域加权截断后的波形和频谱，与图4-1(b)比较，泄漏已明显减少。应该指出，前面所给出的窗函数都是定义为以0点为中心、宽度为N +1的加权函数，在这里应用时，需要将其右移，成为区间内的加权函数。3、对非周期序列进行频谱分析应注意的问题（）混叠一般非周期信号作FFT之前要进行时域采样和周期延拓(无限长时间信号还应先截断再延拓)。根据Fourier变换理论，经等周期的冲激采样后，离散序列的频谱是原信号频谱以为周期的周期延拓。按照Nyquist采样定理，由采样引起周期延拓后频谱

51、之间不发生混叠的条件是：(1)原信号应该是有限带宽信号，设其频带宽度为fm；(2)频谱的周期，即采样周期应满足Nyquist条件。0 n 0 N/2 N(a)时域按周期Ts采样(b)频域一个周期内在N/2附近出现混叠图 4-4非周期函数采样后的幅频曲线由于实际上有限长时间信号的FT是频域的无限函数，因此采样所得的离散序列的频谱必定产生混叠，减小采样周期只能减小而不能消除混叠。对于时间有限函数，当采样周期较大时，也会在FFT得到的频域出现混叠，形成频谱失真，造成频谱分析结果与原信号的实际频谱的差异，也无法恢复出原信号。当然，实际工作中只要采样和截断产生的误差在许可的范围内就行了，但应该认识到混叠

52、是引起频谱分析误差的一个主要原因。还应该注意的是，离散Fourier变换的频域也是周期化的，区间内的样点实际上是负频率区的量值，因此如果出现混叠，就将在一个周期内出现，并发生在附近的区域，如图4-4所示。要减少混叠，就要尽量减小采样周期。（）泄漏周期函数截断引起的频率泄漏问题，在非周期函数截断处理后同样存在，这种误差是由于采样后的离散序列与矩形截断序列相乘在频域造成二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，它与离散序列频谱卷积的结果使原来集中于每一个样点处的频谱展宽，其影响在高频区(接近N/2的样点)特别明显，如图4-5所示。同样，为了尽量避免频率泄漏对结果的影响，在对非周期函数作时间

53、截断时，除尽量增加截断序列的宽度外，也应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。0 n 0 N/2 N(a)时域截断(b)频域一个周期内在N/2附近出现泄漏图4-5函数采样后作截断的幅频曲线在选取了适当的窗函数后，应当使窗函数的宽度与被处理的序列长度相同，如果作变换前还需要补零（例如为了作卷积运算或避免栅栏效应），则应将原序列与窗函数相乘后再补零，即补零的样点不用窗函数加权处理。（）栅栏效应非周期信号应具有连续的频谱，在对作抽样后进行DFT，得到的是离散的频谱。如果排除混叠和泄漏等误差的影响，所得的结果也只是的连续频谱上的个样值。这就象通过栅栏上的等间距缝隙观看到的另一边的景象，故此称

54、栅栏效应。被栅栏遮住的景象中有可能存在与显现出的频谱差异较大的变化，即显示信号特征的频谱分量。为了使被栅栏遮住的部分能尽可能地显现出来，可以采用增加频域样点密度的方法，即在不增加信号采样点的情况下，用时域补零加宽变换尺度N来实现，称为补零重构。例如原来信号采样得到12个样点，在其后面再加上4个零，使序列的总长度为16个样点。这样处理的结果原来信号的采样间隔和频率都没有改变，设采样频率为，经补零重构之后，采样频率仍然为，但是原来频域样点间宽度为/12，经补零重构之后频域样点间宽度为/16。这就使补零重构之后频域样点密度增加，而且显示出原来没有显露的一些频率位置的频谱。三、实验内容（1）将余弦函数

55、cos(2pt)以Ts=1/53 s抽样，对余弦序列做样点数为N=128的FFT，画出频谱曲线，观察并记录频率泄漏现象，然后用Hamming窗和三角窗作加权截断，观察并记录泄漏的衰减。实验代码：n=0:127;N=128;Ts=1./53;t=n.*Ts;xn=cos(2.*pi.*t);w1=hamming(N);w2=bartlett(N);H=xn.*w1'B=xn.*w2'figure(1);subplot(2,1,1);stem(n,xn,'.');title('xn函数')subplot(2,1,2);stem(abs(fft(xn,

56、N),'.');title('xn频谱曲线')figure(2);subplot(2,1,1);stem(n,H,'.');title('Hamming窗加权')subplot(2,1,2);stem(abs(fft(H,N),'.');title('Hamming窗加权频谱曲线')figure(3);subplot(2,1,1);stem(n,B,'.');title('三角窗加权')subplot(2,1,2);stem(abs(fft(B,N),'.');title('三角窗加权频谱曲线')实验截图：（2）将幅度为1，周期为2的方波信号，按Ts=1/37 s的间距抽样，做样点数N=128的FFT，画出频谱曲线，然后用Hamming窗和三角窗作加权截断，观察并记录作不同的加权截断引起的频谱差异。实验代码：n=0:1:127;N=128;Ts=1./37