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

1、数字信号处理实验指导书数字信号处理上机实验指导书陈纯锴电子与信息工程学院一、引言“数字信号处理”是一门理论和实验密切结合的课程，为了深入地掌握课程内容，应当在学习理论的同时，做习题和上机实验。上机实验不仅可以帮助学生深入地理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。所以，根据本课程的重点要求编写了四个实验。第一章是全书的基础内容，抽样定理、时域离散系统的时域和频域分析以及系统对输入信号的响应是重要的基本内容。由于第一章大部分内容已经在前期信号与系统课程中学习完，所以可通过实验一帮助学生温习以上重要内容，加深学生对“数字信号处理是通过对输入信号的一种运算达到处理目的” 这一重要概念

2、的理解。这样便可以使学生从信号与系统课程顺利的过渡到本课程的学习上来。第二章、三章DFT、FFT是数字信号处理的重要数学工具，它有广泛的使用内容。采用实验二、实验三加深理解DFT的基本概念、基本性质。FFT是它的快速算法，必须学会使用。数字滤波器的基本理论和设计方法是数字信号处理技术的重要内容。学习这一部分时，应重点掌握IIR和FIR两种不同的数字滤波器的基本设计方法。IIR滤波器的单位冲激响应是无限长的，设计方法是先设计模拟滤波器，然后再通过SZ平面转换，求出相应的数字滤波器的系统函数。这里的平面转换有两种方法，即冲激响应不变法和双线性变换法，后者没有频率混叠的缺点，且转换简单，是一种普遍应

3、用的方法。FIR滤波器的单位冲激响应是有限长的，设计滤波器的目的即是求出符合要求的单位冲激响应。窗函数法是一种基本的，也是一种重要的设计方法。学习完第七章后可以进行实验四。二、关于使用计算机语言由于数字信号处理实验的主要目的是验证数字信号处理的有关理论，进一步理解巩固所学理论知识。所以，实现实验用的算法语言可以有许多种，但为了提高实验效率，要求学生用编程效率比C语言高好几倍的MATLAB语言。下面介绍MATLAB的主要特点。（有关MATLAB的启动、程序运行和有关信号处理工具箱函数等内容将放到最后附录中介绍。）MATLAB是一种交互式的以矩阵为基本数据结构的系统。在生成矩阵对象时，不要求明确的

4、维数说明。所谓交互式，是指MATLAB的草稿纸编程环境。即用户每输入一条命令并按回车键，MATLAB系统便解释执行之，并显示执行结果。根据该结果，用户立即知道刚输入的命令的正确性，或利用中间结果进行其他处理等。与C语言或FORTRON语言做科学数值计算的程序设计相比较，利用MATLAB可节省大量的编程时间。将其用于数字信号处理实验，则可大大提高实验效率，在有限的上机时间内，实验内容可增加几倍。例如，C语言FFT子程序有70多行，而用MATLAB只调用一个fft函数即可实现对序列进行FFT计算。另外，MATLAB的工具箱及图形显示（打印）功能，可满足各层次人员直观、方便的进行分析、计算和设计工作

5、，从而可大大节省时间。例如，序列的卷积、滤波，系统函数H(z)的幅频特性和相频特性等计算，均有现成的工具箱函数。而用其它算法语言完成这些计算的编程比较麻烦，且程序较长。由于上述特点，在美国一些大学里，MATLAB已成为辅助教学的有益工具。MATLAB已成功地用于数字信号处理课程中的问题分析、实验、滤波器设计及计算机模拟。附录中所介绍的信号处理工具箱函数及绘图函数基本可满足本教材所要求的上机实验需要。对序列进行谱分析的MATLAB程序及运行结果见附录。实验一 离散信号产生及频谱的绘制一、实验目的: （1）熟悉Matlab环境。（2）掌握 Matlab 中一些基本函数的建立方法（2）通过编程绘制的

6、幅度相位谱加深理解系统的特性二、实验内容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('单位脉冲序列');（2）f(n)=u(n) (-5<n<5)n1=-5;n2=5;n0=0;n=n1:n2;x=n>=n0;stem(n,x,'filled&

7、#39;);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);title('复指数信号的实部');subplot(2,2,3);stem(n,real(x),'filled');title('复指数序列的实部');s

8、ubplot(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) (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);sub

9、plot(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）一个连续的周期性方波信号频率为200Hz，信号幅度在-1+1V之间，要求在图形窗口上显示其两个周期的波形。以4kHz的频率对连续信号进行采样，编写程序生成连续信号和其采样获得的离散信号波形。f=200;nt=2; %显示周期数

10、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);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');box（6）绘制的幅度谱和相位谱。clc;%本语句的作用是清除命令执行界面中所有的输

11、出信息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*j*w)+exp(-6*j*w);mag=abs(hjw);%取模/幅度值ang=angle(hjw);%相位figure(1);%第一个窗口subplot(211)%将第一个窗口分成2行1列的子图，分

12、的方法是从左到右，从上到下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','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','3*pi/4',&#

13、39;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;hjw1=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

14、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,'XtickLabel','0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8

15、9;,'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('相位谱');%设定第2幅子图的标题grid on %可以在同一个子图里面叠绘多张图（7）已知系统函数，用MATLAB绘出8阶系统函数的零极点图、幅频响应和相频响应曲线。b=1 0 0 0 0

16、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');ylabel('|H(ejomega)|'); subplot(1,3,3); plot(w/pi,angle(H); %绘制相频响应曲线xlabel('omega/pi

17、9;);ylabel('phi(omega)');实验二 离散傅立叶变换及谱分析一、 实验目的1掌握离散傅里叶变换的计算机实现方法。2检验实序列傅里叶变换的性质。3掌握计算序列的圆周卷积的方法。4熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。5学习用DFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT。二、 实验内容1实现离散傅里叶变换。2计算序列圆周卷积。3计算实序列傅里叶变换并检验DFT性质。4实现连续信号傅里叶变换以及由不同采样频率采样得到的离散信号的傅里叶变换。5实现补零序列的傅里叶变换。6实现高密度谱

18、和高分辨率谱，并比较二者的不同。三、 实验报告要求 见各程序要求%以下为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=(Xk*WNnk)/N;% (3) 实序列的分解 % 实序列可分解为共扼对称分

19、量 % 和共扼反对称分量 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 zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m

20、,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') xlabel('n');ylabel('

21、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); %求序列的共扼反对称分量的DFTsubplot(2,2,1);stem(n,r

22、eal(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('DFTxec(n)');xlabel('k');su

23、bplot(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)>N error('N must be >= the length of x1')endif length(x2)>N error('N must be >= the length

24、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 本例验证采样定理%令，绘制其傅立叶变换。用不同频率对其进行采样，分别

25、画出离散时间傅立叶变换。已给出采样频率为时的的程序 %实验报告要求：（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'*W)*Dt; %计算连续时间傅立叶变换（利用矩阵运算实现

26、） 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');subplot(2,2,2);plot(W/(2*pi

27、*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); %计算离散时间傅立叶变换（序列的傅立叶变换）X=real(X); w=-fliplr(w

28、),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);plot(w/pi,X); %画出离散时间傅立叶变换

29、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*pi/500).(n'*k); %计算离

30、散时间傅立叶变换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);subplot(3,2,3);stem(n,x);yl

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

32、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);stem(w1/pi,abs(Y);title(

33、'信号的频谱');实验三 用FFT对信号作频谱分析一、实验目的（1）学习用FFT对连续信号和时域离散信号进行谱分析的方法（2）了解可能出现的分析误差及其原因，以便正确应用FFT。二、实验原理用FFT对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是，因此要求。可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近

34、于连续谱，因此N要适当选择大一些。周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。三、实验步骤及内容（1）对以下序列进行谱分析。 选择FFT的变换区间N为8和16 两种情况进行频谱分析。分别打印其幅频特性曲线。 并进行对比、分析和讨论。（2）对以下周期序列进行谱分析。 选择FFT的变换区间N为8和16 两种情况分别对以上序列进行频谱分析。

35、分别打印其幅频特性曲线。并进行对比、分析和讨论。（3）对模拟周期信号进行谱分析 选择采样频率，变换区间N=16,32,64 三种情况进行谱分析。分别打印其幅频特性，并进行分析和讨论。 四、思考题（1）对于周期序列，如果周期不知道，如何用FFT进行谱分析？（2）如何选择FFT的变换区间？（包括非周期信号和周期信号）（3）当N=8时，和的幅频特性会相同吗？为什么？N=16 呢？五、实验报告要求（1）完成各个实验任务和要求。附上程序清单和有关曲线。（2）简要回答思考题。六、实验程序清单 % 用FFT对信号作频谱分析clear all;close all%实验内容(1)=x1n=ones(1,4);

36、%产生序列向量x1(n)=R4(n)M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb; %产生长度为8的三角波序列x2(n)x3n=xb,xa;X1k8=fft(x1n,8); %计算x1n的8点DFTX1k16=fft(x1n,16); %计算x1n的16点DFTX2k8=fft(x2n,8); %计算x2n的8点DFTX2k16=fft(x2n,16); %计算x2n的16点DFTX3k8=fft(x3n,8); %计算x3n的8点DFTX3k16=fft(x3n,16); %计算x3n的16点DFT%以下绘制幅频特性曲线subplot(2,2,1);mst

37、em(X1k8); %绘制8点DFT的幅频特性图title('(1a) 8点DFTx_1(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k8)subplot(2,2,3);mstem(X1k16); %绘制16点DFT的幅频特性图title('(1b)16点DFTx_1(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k16)figure(2)subplot(2,2,1);m

38、stem(X2k8); %绘制8点DFT的幅频特性图title('(2a) 8点DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k8)subplot(2,2,2);mstem(X2k16); %绘制16点DFT的幅频特性图title('(2b)16点DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k16)subplot(2,2,3);mstem(X3

39、k8); %绘制8点DFT的幅频特性图title('(3a) 8点DFTx_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k8)subplot(2,2,4);mstem(X3k16); %绘制16点DFT的幅频特性图title('(3b)16点DFTx_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k16)%实验内容(2) 周期序列谱分析=N=8;n=0:N-1;

40、%FFT的变换区间N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n,8); %计算x4n的8点DFTX5k8=fft(x5n,8); %计算x5n的8点DFTN=16;n=0:N-1; %FFT的变换区间N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n); %计算x4n的16点DFTX5k16=fft(x5n); %计算x5n的16点DFTfigure(3)subplot(2,2,1);mstem(X4k8); %绘制8点DFT的幅频特性图title(

41、'(4a) 8点DFTx_4(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); %绘制16点DFT的幅频特性图title('(4b)16点DFTx_4(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k16)subplot(2,2,2);mstem(X5k8); %绘制8点DFT的幅频特性图title('(5

42、a) 8点DFTx_5(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X5k8)subplot(2,2,4);mstem(X5k16); %绘制16点DFT的幅频特性图title('(5b)16点DFTx_5(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X5k16)%实验内容(3) 模拟周期信号谱分析=figure(4)Fs=64;T=1/Fs;N=16;n=0:N-1; %FFT的变换

43、区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样X6k16=fft(x6nT); %计算x6nT的16点DFTX6k16=fftshift(X6k16); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %绘制8点DFT的幅频特性图title('(6a) 16点|DFTx_6(nT

44、)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16)N=32;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样X6k32=fft(x6nT); %计算x6nT的32点DFTX6k32=fftshift(X6k32); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点

45、DFT对应的采样点频率（以零频率为中心）subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %绘制8点DFT的幅频特性图title('(6b) 32点|DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32)N=64;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点

46、采样X6k64=fft(x6nT); %计算x6nT的64点DFTX6k64=fftshift(X6k64); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on%绘制8点DFT的幅频特性图title('(6a) 64点|DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-

47、1,N*F/2-1,0,1.2*max(abs(X6k64)%= mstem程序清单=function mstem(Xk)% mstem(Xk)绘制频域采样序列向量Xk的幅频特性图M=length(Xk);k=0:M-1;wk=2*k/M; %产生M点DFT对应的采样点频率(关于归一化值)stem(wk,abs(Xk),'.');box on %绘制M点DFT的幅频特性图xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(Xk)七、实验程序运行结果运行结果如图3所示。图3程序运行结果分析讨论：请读者注意

48、，用DFT（或FFT）分析频谱，绘制频谱图时，最好将X(k)的自变量k换算成对应的频率，作为横坐标便于观察频谱。为了便于读取频率值，最好关于归一化，即以作为横坐标。1、实验内容（1）图（1a）和（1b）说明的8点DFT和16点DFT分别是的频谱函数的8点和16点采样；因为，所以，与的8点DFT的模相等，如图（2a）和（3a）。但是，当N=16时，与不满足循环移位关系，所以图（2b）和（3b）的模不同。2、实验内容（2），对周期序列谱分析的周期为8，所以N=8和N=16均是其周期的整数倍，得到正确的单一频率正弦波的频谱，仅在0.25处有1根单一谱线。如图（4b）和（4b）所示。的周期为16，所以

49、N=8不是其周期的整数倍，得到的频谱不正确，如图（5a）所示。N=16是其一个周期，得到正确的频谱，仅在0.25和0.125处有2根单一谱线, 如图（5b）所示。3、实验内容（3），对模拟周期信号谱分析有3个频率成分，。所以的周期为0.5s。 采样频率。变换区间N=16时，观察时间Tp=16T=0.25s，不是的整数倍周期，所以所得频谱不正确，如图（6a）所示。变换区间N=32,64 时，观察时间Tp=0.5s，1s，是的整数周期，所以所得频谱正确，如图（6b）和（6c）所示。图中3根谱线正好位于处。变换区间N=64 时频谱幅度是变换区间N=32 时2倍，这种结果正好验证了用DFT对中期序列谱

50、分析的理论。注意：（1）用DFT（或FFT）对模拟信号分析频谱时，最好将X(k)的自变量k换算成对应的模拟频率fk，作为横坐标绘图，便于观察频谱。这样，不管变换区间N取信号周期的几倍，画出的频谱图中有效离散谐波谱线所在的频率值不变，如图（6b）和（6c）所示。（2）本程序直接画出采样序列N点DFT的模值，实际上分析频谱时最好画出归一化幅度谱，这样就避免了幅度值随变换区间N变化的缺点。本实验程序这样绘图只要是为了验证了用DFT对中期序列谱分析的理论。实验四 FIR数字滤波器设计与软件实现一、实验目的（1）掌握用窗函数法设计FIR数字滤波器的原理和方法。（2）掌握用等波纹最佳逼近法设计FIR数字滤

51、波器的原理和方法。（3）掌握FIR滤波器的快速卷积实现原理。（4）学会调用MATLAB函数设计与实现FIR滤波器。二、实验原理与方法（1）FIR滤波器的设计在前面的实验中，我们介绍了IIR滤波器的设计方法并实践了其中的双线性变换法，IIR具有许多诱人的特性；但如此同时，也具有一些缺点。例如：若想利用快速傅立叶变换技术进行快速卷积实现滤波器，则要求单位脉冲响应是有限长的。此外，IIR滤波器的优异幅度响应，一般是以相位的非线性为代价的，非线性相位会引起频率色散。FIR滤波器具有严格的相位特性，这对于语音信号处理和数据传输是很重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率取样法和切比雪

52、夫等波纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。本实验中的窗函数法比较简单，可应用现成的窗函数公式，在技术指标要求不高的时候是比较灵活方便的。它是从时域出发，用一个窗函数截取理想的得到，以有限长序列近似理想的；如果从频域出发，用理想的在单位圆上等角度取样得到，根据得到将逼近理想的，这就是频率取样法。（2）窗函数设计法同其它的数字滤波器的设计方法一样，用窗函数设计滤波器也是首先要对滤波器提出性能指标。一般是给定一个理想的频率响应，使所设计的FIR滤波器的频率响应去逼近所要求的理想的滤波器的响应。窗函数法设计的任务在于寻找一个可实现（有限长单位脉冲响应）的传递函

53、数去逼近。一个理想的频率响应的傅立叶反变换所得到的理想单位脉冲响应往往是一个无限长序列。对经过适当的加权、截短处理才能得到一个所需要的有限长脉冲响应序列。对应不同的加权、截短，就有不同的窗函数。所要寻找的滤波器脉冲响应就等于理想脉冲响应和窗函数的乘积，即 由此可见，窗函数的形状就决定了滤波器的性质。例如：窗函数的主瓣宽度决定了滤波器的过渡带宽；窗函数的旁瓣大小决定了滤波器的阻带衰减。1、几种常见的窗函数：（1）矩形窗（Rectangle Window）调用格式：w=boxcar(n)，根据长度n 产生一个矩形窗w。（2）三角窗（Triangular Window）调用格式：w=triang(n

54、) ，根据长度n 产生一个三角窗w。（3）汉宁窗（Hanning Window）调用格式：w=hanning(n) ，根据长度n 产生一个汉宁窗w。（4）海明窗（Hamming Window）调用格式：w=hamming(n) ，根据长度n 产生一个海明窗w。（5）布拉克曼窗（Blackman Window）调用格式：w=blackman(n) ，根据长度n 产生一个布拉克曼窗w。（6）恺撒窗（Kaiser Window）调用格式：w=kaiser(n,beta) ，根据长度n 和影响窗函数旁瓣的参数产生一个恺撒窗w。三、实验内容及步骤（1）认真复习第六章中用窗函数法和等波纹最佳逼近法设计FI

55、R数字滤波器的原理；（2）调用信号产生函数xtg产生具有加性噪声的信号xt，并自动显示xt及其频谱，如图4.1所示；图4.1 具有加性噪声的信号x(t)及其频谱如图（3）请设计低通滤波器，从高频噪声中提取xt中的单频调幅信号，要求信号幅频失真小于0.1dB，将噪声频谱衰减60dB。先观察xt的频谱，确定滤波器指标参数。（4）根据滤波器指标选择合适的窗函数，计算窗函数的长度N，调用MATLAB函数fir1设计一个FIR低通滤波器。并编写程序，调用MATLAB快速卷积函数fftfilt实现对xt的滤波。绘图显示滤波器的频响特性曲线、滤波器输出信号的幅频特性图和时域波形图。（5）重复（3），滤波器指

56、标不变，但改用等波纹最佳逼近法，调用MATLAB函数remezord和remez设计FIR数字滤波器。并比较两种设计方法设计的滤波器阶数。提示：1）、MATLAB函数fir1和fftfilt的功能及其调用格式请查阅本书第6章；2）、采样频率Fs=1000Hz，采样周期T=1/Fs；3）、根据图4.1(b)和实验要求，可选择滤波器指标参数：通带截止频率fp=120Hz，阻带截至频率fs=150Hz，换算成数字频率，通带截止频率，通带最大衰为0.1dB，阻带截至频率，阻带最小衰为60dB。4）、实验程序框图如图4.2所示，供读者参考。Fs=1000，T=1/Fsxt=xtg产生信号xt, 并显示xt及其频谱用窗函数法或等波纹最佳逼近法设计FIR滤波器hn对信号xt滤波：yt=fftfilt(hn,xt)1、计算并绘图显示滤波器损耗函数2、绘图显示滤波器输出信号ytEnd图6.2 实验程序框图四、思考题(1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函数法设计线性相位低通滤波器?请写出设计步骤。(2)如果要求用窗函数法设计带通滤波器,且给定通带上、下截止频率为和，阻带上、下截止频率为和，试求理想带通滤波器的截止频率。（3）解释为什么对同样的技术指