数字信号处理实验指导书(正1)

数字信号处理实验指导书长沙理工大学电气与信息工程学院电子信息工程教研室二零零柒年柒月目录TOC\o"1-3"\h\z实验一：离散时间信号的时域分析 3实验二：离散时间系统的时域分析 6实验三：离散时间信号的频域分析 9实验四：线性时不变离散时间系统的频域分析 13实验五：IIR数字滤波器的设计 17实验六：FIR数字滤波器的设计 24附录AMATLAB系统的常用概念 28附录B信号处理工具箱函数 33

数字信号处理研究数字序列信号的表示方法，并对信号进行运算，以提取包含在其中的特殊信息。近几十年来，由于在研究及应用两方面均取得了进展，数字信号处理领域已日趋成熟。由于计算机的大量使用，从而很容易向学生提供实际环境，以验证所学的概念和算法。本指导书编程语言是MATLAB，它广泛应用于高性能数值计算和可视化。本书假定读者已具备MATLAB基础知识。前面的一些实验帮助学生理解信号处理的重要概念，后面以设计性实验项目为主，有利于加强对原理的理解并且加强对技术的应用。附录中给出了本书中用到的MATALB函数及简要解释。实验一：离散时间信号的时域分析实验目的熟悉MATLAB中产生信号和绘制信号的基本命令。熟悉序列的简单运算，如：加法、标量乘法、时间反转、延时、乘法等。二．实验相关知识准备1．用到的MATLAB命令运算符：：.+-*/;%基本矩阵：ionespirandrandnzeros基本函数：cosexpimagreal数据分析：sum二维图形：axisgridlegendplotsubplotstairsstemtitlexlableylableclf工具箱：sawtoothsquare三．实验内容1．序列的产生(1)程序1－1：单位抽样序列的产生和绘制%ProgramP1_1%GenerationofaUnitSampleSequenceclf;%Generateavectorfrom-10to20n=-10:20;%Generatetheunitsamplesequenceu=[zeros(1,10)1zeros(1,20)];%Plottheunitsamplesequencestem(n,u);xlabel('Timeindexn');ylabel('Amplitude');title('UnitSampleSequence');axis([-102001.2]);(2)程序1－2：正弦序列的产生和绘制%ProgramP1_2%Generationofasinusoidalsequencen=0:40; f=0.1; phase=0; A=1.5; arg=2*pi*f*n-phase;x=A*cos(arg);clf; %Clearoldgraphstem(n,x); %Plotthegeneratedsequenceaxis([040-22]);grid;title('SinusoidalSequence');xlabel('Timeindexn');ylabel('Amplitude');axis;2．序列的运算（信号的平滑）数字信号处理应用的一个常见例子是从被加性噪声污染的信号中移除噪声。假定信号s[n]被噪声d[n]所污染，得到一个含有噪声的信号x[n]=s[n]+d[n]。我们需要对x[n]进行运算，产生一个合理的逼近s[n]，对时刻n的样本求平均，产生输出信号是一种简单有效的方法。如：三点滑动平均的信号。程序1－3实现三点滑动平均的信号运算：y[n]=(x(n-1)+x(n)+x(n+1))/3程序如下:%SignalSmoothingbyAveragingclf;R=51;d=0.8*(rand(R,1)-0.5);%产生随机噪声m=0:R-1;s=2*m.*(0.9.^m);%产生未被污染的信号x=s+d';%产生被噪声污染的信号subplot(2,1,1);plot(m,d','r-',m,s,'g--',m,x,'b-.');xlabel('Timeindexn');ylabel('Amplitude');legend('d[n]','s[n]','x[n]');x1=[00x];x2=[0x0];x3=[x00];y=(x1+x2+x3)/3;subplot(2,1,2);plot(m,y(2:R+1),'r-',m,s,'g--');legend('y[n]','s[n]');xlabel('Timeindexn');ylabel('Amplitude');3．实验问题回答命令clf,axis,title,xlable和ylable的作用是什么？程序1－2中正弦序列的频率是多少？怎样可以改变它？哪个参数控制该序列的振幅？该序列周期是多少？程序1－3中加性噪声d[n]是什么样的形式？语句x=s+d代表什么？信号x1,x2,x3与信号x之间的关系是什么？四．实验报告要求按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。回答实验中提出的问题。总结本次实验结果，按照实验报告格式要求，书写实验报告。五．实验设备PC机，MATLAB软件实验二：离散时间系统的时域分析一．实验目的通过MATLAB仿真一些简单的离散时间系统，并研究它们的时域特性。二．实验相关知识准备1．用到的MATLAB命令语言构造与调试：breakendforifinput基本函数：absnum2str多项式和内插函数：conv工具箱：filterimpz三．实验内容1．线性和非线性系统例2－1设系统为y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2]要求用MATLAB程序仿真系统，输入三个不同的输入序列x1(n),x2(n)和x(n)=a.x1(n)+b.x2(n)，计算并求出相应的输出响应y1[n],y2[n]和y[n]。%Generatetheinputsequencesclf;n=0:40;a=2;b=-3;x1=cos(2*pi*0.1*n);x2=cos(2*pi*0.4*n);x=a*x1+b*x2;num=[2.24032.49082.2403];den=[1-0.40.75];ic=[00];%Setzeroinitialconditionsy1=filter(num,den,x1,ic);%Computetheoutputy1[n]y2=filter(num,den,x2,ic);%Computetheoutputy2[n]y=filter(num,den,x,ic);%Computetheoutputy[n]yt=a*y1+b*y2;d=y-yt;%Computethedifferenceoutputd[n]%Plottheoutputsandthedifferencesignalsubplot(3,1,1)stem(n,y);ylabel('Amplitude');title('OutputDuetoWeightedInput:a\cdotx_{1}[n]+b\cdotx_{2}[n]');subplot(3,1,2)stem(n,yt);ylabel('Amplitude');title('WeightedOutput:a\cdoty_{1}[n]+b\cdoty_{2}[n]');subplot(3,1,3)stem(n,d);xlabel('Timeindexn');ylabel('Amplitude');title('DifferenceSignal');2．时不变系统和时变系统实验内容2－2：因果系统的时不变特性研究编写MATLAB程序仿真例2.1给出的系统，以产生两个不同的输入序列x(n)（在例2－1中x(n)=a.x1(n)+b.x2(n)）和x(n+D)，计算并画出相应的输出序列y1[n],y2[n+D]和y1[n]-y2[n+D]参数：n=0:40;D=10;a=3.0;b=-2;实验内容2－3：线性时不变系统的冲激响应计算设系统为y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2]用y=impz(num,den,N)计算系统的冲激响应的前N个样本。3.实验问题回答修改程序2－1的相关语句，对由加权输入得到的y[n]与相同权系数下输出y1[n]和y2[n]相加得到的yt[n]进行比较，这两个序列是否相等？该系统是线性系统吗？比较实验内容2－2中输出序列y[n]和yd[n-10]，该系统是时不变系统吗？计算并画出实验内容2－3的前40个样本。四．实验报告要求1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。2．回答实验中提出的问题。总结本次实验结果，按照实验报告格式要求，书写实验报告。五．实验设备PC机，MATLAB软件实验三：离散时间信号的频域分析一．实验目的1．在学习了离散时间信号的时域分析的基础上，对这些信号在频域上进行分析，从而进一步研究它们的性质。2．熟悉离散时间序列的3种表示方法：离散时间傅立叶变换（DTFT）,离散傅立叶变换（DFT）和Z变换。二．实验相关知识准备用到的MATLAB命令运算符和特殊字符：<>.*^.^语言构造与调试：errorfunctionpause基本函数：angleconjrem数据分析和傅立叶变换函数：fftifftmaxmin工具箱：freqzimpzresiduezzplane三．实验内容离散傅立叶变换在MATLAB中，使用fft可以很容易地计算有限长序列x[n]的离散傅立叶变换。此函数有两种形式：y=fft(x)y=fft(x,n)求出时域信号x的离散傅立叶变换n为规定的点数，n的默认值为所给x的长度。当n取2的整数幂时变换的速度最快。通常取大于又最靠近x的幂次。（即一般在使用fft函数前用n=2^nextpow2(length(x))得到最合适的n)。当x的长度小于n时，fft函数在x的尾部补0，以构成长为n点数据。当x的长度大于n时，fft函数将序列x截断，取前n点。一般情况下，fft求出的函数多为复数，可用abs及angle分别求其幅度和相位。注意：栅栏效应，截断效应（频谱泄露和谱间干扰），混叠失真例3－1：fft函数最通常的应用是计算信号的频谱。考虑一个由100hz和200hz正弦信号构成的信号，受零均值随机信号的干扰，数据采样频率为1000hz。通过fft函数来分析其信号频率成分。t=0:0.001:1;%采样周期为0.001s，即采样频率为1000hzx=sin(2*pi*100*t)+sin(2*pi*200*t)+1.5*rand(1,length(t));%产生受噪声污染的正弦波信号subplot(2,1,1);plot(x(1:50));%画出时域内的信号y=fft(x,512);%对x进行512点的fftf=1000*(0:256)/512;%设置频率轴（横轴）坐标，1000为采样频率subplot(2,1,2);plot(f,y(1:257));%画出频域内的信号实验内容3－2：混叠现象与频率分辨率假设现有含有三种频率成分的信号x(n)=sin(2nπf1/fs)+sin(2nπf2/fs)+sin(2nπf3/fs),其中f1=20hz,f2=20.5hz,f3=40hz,采样频率fs=100hz.（1）求x(n)在0～128之间的DFT的x(k)；（2）求把（1）中的x(n)以补0方式使其加长到0～512之间后的DFT的x(k)；（3）求x(n)在0～512之间的DFT的x(k)。最后分析三种形式的DFT结果图给出结论。提示：根据采样定理，在信号采样的过程中，如果采样频率小于信号最高频率的两倍，那么信号的各次调制频谱就会相互重叠起来，有些频率部分的幅值就与原始情况不同，因而不能分开和恢复这些部分，取样造成了信号的损失，这种频谱重叠的出现，就成为“混叠现象”。令为采样频率，为信号的最高频率，为了避免混叠现象，则必须使，而频谱分辨率，即在频率轴上所能得到的最小频率间隔F又为。所以，在和两个参数中，保持其中一个不变而增加另一个的唯一方法，就是增加在一记录长度内的点数N。如果与都已给定，则N必须满足一维逆快速傅立叶变换y=ifft(x)y=ifft(x,n)实验内容3-3：频域采样定理的验证产生一个三角波序列x(n)(2)对M=40，计算x(n)的64点DFT，并图示x(n)和X(k)=DFT[x(n)],k=0,1,…,63。(3)对（2）中所得X(k)在[0,2pi]上进行32点抽样得X1(k)=X(2k),k=0,1,…,31。并图示。（4）求X1（k）的32点IDFT,并图示。即x1(n)=IDFT[X1(k)],k=0,1,…,31(5)采用周期延拓的方法绘出的波形，评述与x(n)的关系，并根据频域采样理论加以解释。线性调频Z变换离散傅立叶变换（DFT）可以看作信号在Z域上沿单位圆的均匀采样。但在实际应用中，并非整个单位圆上的频谱都有意义。一些情况下，如对于窄带信号，只希望分析信号所在的一段频带等，采样点的轨迹是一条弧线或圆周。这种需求，就导致了线性调频Z变换（Chirpz变换）的出现。Chirpz变换与DFT计算整个频谱的算法不同，它是一种更为灵活的计算频谱的算法，可以用来计算单位圆上任一段曲线的Z变换，作频谱分析时输入的点数和输出的点数可以不相等，从而达到频域“细化”的目的。y=czt(x,m,w,a)y=czt(x)例3-4：利用Chirpz变换计算滤波器h（且h=fir1(30,125/500,boxcar(31))）在100Hz~200Hz的频率特性，并用图形文件比较CZT函数和FFT函数。h=fir1(30,125/500,boxcar(31));fs=1000;f1=100;f2=200;m=1024;y=fft(h,1024);fy=fs*(0:1023)/1024;subplot(2,1,1);plot(fy,abs(y));axis([0,500,0,1.5]);w=exp(-j*2*pi*(f2-f1)/(m*fs));a=exp(j*2*pi*f1/fs);z=czt(h,m,w,a);%产生1024个z值fz=(f2-f1)*(0:1023)/1024+f1;subplot(2,1,2);plot(fz,abs(z))实验问题回答（1）完成例3－1的程序并观察信号序列x(t)在时域和频域上分析的特点。（2）在实验内容3-2和3-3中按要求完成程序，并回答3-2和3-3中提出的问题。完成例3－4的程序并观察Chirpz变换与fft变换的特点。四．实验报告要求1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。回答实验中提出的问题。总结本次实验结果，按照实验报告格式要求，书写实验报告。五．实验设备PC机，MATLAB软件实验四：线性时不变离散时间系统的频域分析一．实验目的1．在学习了离散时间系统的时域分析的基础上，对这些系统在频域上进行分析，掌握LTI系统的频域分析法。2．熟悉传输函数和频率响应。二．实验相关知识准备1．用到的MATLAB命令MATLAB的工具箱，为同一离散系统的不同系统模型间的多种表达方式的变换，提供了很多功能丰富快捷的MATLAB函数，如下表所示：系统原模型系统转换模型函数传递函数状态空间Tf2ss传递函数零极点增益Tf2zp零极点增益状态空间Zp2ss零极点增益传递函数Zp2tf零极点增益二次分式Zp2sos状态空间传递函数Ss2tf状态空间零极点增益Ss2zp状态空间二次分式Ss2sos二次分式零极点增益Sos2tf二次分式状态空间Sos2zp多项式Lattice结构Sos2ssLattice结构多项式Poly2rc基本函数：absangleimaglog10realzplane可以查看零极点poly2rc可以直接检测稳定性三．实验内容例4－1求截短的理想低通滤波器的冲激响应理想低通滤波器的冲激响应是无限的，因此实际上不可能实现。在实际应用中，将冲激响应截短到一个有限值，可以近似地实现该滤波器。然而，将冲激响应截短后滤波器表现为一个非因果的滤波器。将经过截短后的滤波器的冲激响应向右移动N/2个样本，可得因果上的近似，如下所示：滤波器的长度为N+1。下面这个使用函数sinc的程序，可用来计算上述近似滤波器的冲激响应。clf;fc=0.25;n=[-6.5:1:6.5];y=2*fc*sinc(2*fc*n);k=n+6.5;stem(k,y);title('N=13');axis([013-0.20.6]);xlabel('Timeindexn');ylabel('Amplitude');grid;实验内容4－2=1\*GB3①运行例4－1的程序计算并画出近似理想低通滤波器的冲激响应。=2\*GB3②在例4－1中的低通有限冲激响应滤波器的长度是多少？该程序中哪个语句确定滤波器的长度？哪个参数控制截止频率？=3\*GB3③修改例4－1中的程序，计算长度为20，截止角频率为0.45的有限冲激响应低通滤波器的冲激响应。（h[n]不变）2．例4－3线性相位有限冲激响应滤波器的零点的位置研究。下列程序先产生四类滤波器的冲激响应，然后产生零极点图，最终显示零点的位置。clf;b=[1-8.530.5-63];num1=[b81fliplr(b)];num2=[b8181fliplr(b)];num3=[b0-fliplr(b)];num4=[b81-81-fliplr(b)];n1=0:length(num1)-1;n2=0:length(num2)-1;subplot(2,2,1);stem(n1,num1);xlabel('Timeindexn');ylabel('Amplitude');grid;title('Type1FIRFilter');subplot(2,2,2);stem(n2,num2);xlabel('Timeindexn');ylabel('Amplitude');grid;title('Type2FIRFilter');subplot(2,2,3);stem(n1,num3);xlabel('Timeindexn');ylabel('Amplitude');grid;title('Type3FIRFilter');subplot(2,2,4);stem(n2,num4);xlabel('Timeindexn');ylabel('Amplitude');grid;title('Type4FIRFilter');pausesubplot(2,2,1);zplane(num1,1);title('Type1FIRFilter');subplot(2,2,2);zplane(num2,1);title('Type2FIRFilter');subplot(2,2,3);zplane(num3,1);title('Type3FIRFilter');subplot(2,2,4);zplane(num4,1);title('Type4FIRFilter');disp('ZerosofType1FIRFilterare');disp(roots(num1));disp('ZerosofType2FIRFilterare');disp(roots(num2));disp('ZerosofType3FIRFilterare');disp(roots(num3));disp('ZerosofType4FIRFilterare');disp(roots(num4));实验内容4－4=1\*GB3①运行例4－3的程序，生成每一类线性相位有限冲激响应滤波器的冲激响应。并回答每一个有限冲激响应滤波器的长度是多少？=2\*GB3②在例4－3中验证冲激响应序列的对称性，每一个滤波器的零点位置及它们的线性相位特性。并回答这些滤波器的群延迟是多少？四．实验报告要求按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。回答实验中提出的问题。总结本次实验结果，按照实验报告格式要求，书写实验报告。五．实验设备PC机，MATLAB软件实验五：IIR数字滤波器的设计一．实验目的1．加深对IIR滤波器的理解2．学会用直接法和间接法设计满足响应要求的IIR数字滤波器二．实验相关知识准备1．模拟滤波器的设计函数为了利用模拟滤波器设计数字滤波器，先必须设计对应的模拟滤波器，常用的模拟滤波器有：bessel滤波器，butterworth滤波器，chebyshev=1\*ROMANI型滤波器，chebyshev=2\*ROMANII型滤波器及椭圆函数滤波器。1》．模拟滤波器的设计函数设计bessel模拟低通滤波器[z,p,k]=besselap(n)2）设计butterworth模拟低通滤波器[z,p,k]=buttap(n)3）设计chebyshev=1\*ROMANI型模拟低通滤波器[z,p,k]=cheb1ap(n,Rp)Rp：通带内的波纹系数，单位分贝4）设计chebyshev=2\*ROMANII型模拟低通滤波器[z,p,k]=cheb2ap(n,Rs)Rs：阻带内的波纹系数低于通带Rs分贝5）设计椭圆模拟滤波器[z,p,k]=ellipap(n,Rp.Rs)2》滤波器阶数的选择下列函数除了能选择模拟滤波器的阶数外，同时也能选择数字滤波器的阶数。选择butterworth滤波器阶数数字域：[n,Wn]=buttord(Wp,Ws,Rp,Rs)模拟域：[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)选择chebyshev=1\*ROMANI型滤波器阶数数字域：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)模拟域：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs,’s’)选择chebyshev=2\*ROMANII型滤波器阶数数字域：[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs)模拟域：[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs,’s’)选择椭圆滤波器阶数数字域：[n,Wn]=ellipord(Wp,Ws,Rp,Rs)模拟域：[n,Wn]=ellipord(Wp,Ws,Rp,Rs,’s’)注意：n：返回符合要求性能指标的数字滤波器或模拟滤波器的最小阶数Wn：滤波器的截至频率（即3db频率）Wp：通带的截至频率，Ws：阻带的截至频率，单位rad/s。且均为规一化频率，即。1对应π弧度。频率规一化：信号处理工具箱中使用的频率为莱奎斯特频率，根据香农定理，它为采样频率的一半，在滤波器设计中的截止频率均使用莱奎斯特频率进行规一化。规一化频率转换为角频率，则将规一化频率乘以π。如果将规一化频率转换为Hz，则将规一化频率乘以采样频率的一半。3》．模拟频率变换低通到低通模拟滤波器变换[bt,at]=lp2lp(b,a,Wo)lp2lp函数将截至频率为1rad/s（规一化截至频率）的模拟低通滤波器原型变换成截至频率为Wo的低通滤波器。低通到带通模拟滤波器变换[bt,at]=lp2bp(b,a,Wo,Bw)lp2bp函数可将截止频率为1rad/s的模拟低通滤波器原型变换成具有指定带宽Bw和中心频率Wo的带通滤波器。其中Wo=sqrt(W1*W2)，Bw=W2-W1，W2高端截止频率，W1低端截止频率。低通到高通模拟滤波器变换[bt,at]=lp2hp(b,a,Wo)低通到带阻模拟滤波器变换[bt,at]=lp2bs(b,a,Wo,Bw)5）besself模拟滤波器设计[b,a]=besself(n,Wn)[b,a]=besself(n,Wn,’ftype’)[z,p,k]=besself(n,Wn)[z,p,k]=besself(n,Wn,’ftype’)2．由模拟滤波器变换成等效的数字滤波器方法一：双线性变换法实现模拟滤波器到数字滤波器的变换双线性变换将S域映射成Z域，从而将模拟滤波器变换成等效的数字滤波器[zd,pd,kd]=bilinear(z,p,k,Fs)为零极点增益表示的bilinear函数其中z,p,k为S域传递函数的零点、极点、和增益，Fs为取样频率，zd,pd,kd为双线性变换后Z域传递函数的零点、极点和增益。[numd,dend]=bilinear(num,den,Fs,Fp)为传输函数表示的bilinear函数方法二：冲激响应不变法实现模拟到数字滤波器的变换[bz,az]=impinvar(b,a,Fs)[bz,az]=impinvar(b,a)采用默认值为1hz的Fs。小结：设计IIR数字滤波器的一般步骤：把给出的数字滤波器的性能指标转换为模拟滤波器的性能指标根据转换后的性能指标，通过滤波器阶数选择函数，来确定滤波器的最小阶数N和固有频率Wn由最小阶数N得到低通滤波器原型由固有频率Wn把模拟低通滤波器原型转换为低通、高通、带通、带阻滤波器运用脉冲响应不变法或双线性变换法把模拟滤波器转换成数字滤波器MATLAB信号处理工具箱提供了几个用于直接设计IIR数字滤波器的函数，这些函数把上面这些复杂的步骤融合为一个整体，为设计滤波器带来了极大的方便。3．直接设计IIR数字滤波器1》Butterworth模拟和数字滤波器设计数字域：[b,a]=butter(n,Wn)可设计出截止频率为Wn的n阶butterworth滤波器[b,a]=butter(n,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器[z,p,k]=butter(n,Wn)[zp,k]=buter(n,Wn,’ftype’)[A,B,C,D]=butter(n,Wn)[A,B,C,D]=butter(n,Wn,’ftype’)模拟域：[b,a]=butter(n,Wn,’s’)可设计出截止频率为Wn的n阶模拟butterworth滤波器，其余形式类似于数字域的。2》chebyshev=1\*ROMANI型滤波器（通带等波纹）设计数字域：[b,a]=cheby1(n,Rp,Wn)可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定[b,a]=cheby1(n,Rp,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器[z,p,k]=cheby1(n,Rp,Wn)[zp,k]=cheby1(n,Rp,Wn,’ftype’)[A,B,C,D]=cheby1(n,Rp,Wn)[A,B,C,D]=cheby1(n,Rp,Wn,’ftype’)模拟域：[b,a]=cheby1(n,Rp,Wn,’s’)可设计出截止频率为Wn的n阶chebyshevI型模拟滤波器，其余形式类似于数字域的。3》chebyshev=2\*ROMANII型滤波器（阻带等波纹）设计数字域：[b,a]=cheby2(n,Rs,Wn)可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，阻带内的波纹由Rs确定[b,a]=cheby2(n,Rs,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器[z,p,k]=cheby2(n,Rs,Wn)[zp,k]=cheby2(n,Rs,Wn,’ftype’)[A,B,C,D]=cheby2(n,Rs,Wn)[A,B,C,D]=cheby2(n,Rs,Wn,’ftype’)模拟域：[b,a]=cheby2(n,Rs,Wn,’s’)可设计出截止频率为Wn的n阶chebyshev=2\*ROMANII型模拟滤波器，其余形式类似于数字域的。4》椭圆滤波器设计数字域：[b,a]=ellip(n,Rp,Rs,Wn)可设计出n阶椭圆滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定,阻带内的波纹由Rs确定。[b,a]=ellip(n,Rp,Rs,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器[z,p,k]=ellip(n,Rp,Rs,Wn)[zp,k]=ellip(n,Rp,Rs,Wn,’ftype’)[A,B,C,D]=ellip(n,Rp,Rs,Wn)[A,B,C,D]=ellip(n,Rp,Rs,Wn,’ftype’)模拟域：[b,a]=ellip(n,Rp,Rs,Wn,’s’)可设计出截止频率为Wn的n阶椭圆型模拟滤波器，其余形式类似于数字域的。5》基于幅度特性设计IIR数字滤波器[b,a]=yulewalk(n,f,m)可得到以b,a表示的n阶IIR滤波器，其幅频特性逼近由频率矢量f和幅度矢量m所给定的特性。矢量f和n的长度必须相同。6》基于频率特性设计IIR数字滤波器[b,a]=invfreqz(h,w,nb,na)可得到以b,a表示的IIR滤波器，其幅频特性逼近由频率矢量w和复频率响应矢量h所给定的特性。矢量w和h的长度必须相同，nb,na分别为指定的分子和分母多项式的阶数。三．实验内容1．例5－1设计一butterworth低通滤波器，通带范围0～100hz，通带波纹小于3db，阻带为-30db，并利用最小的阶数来实现。采样频率为1000hz。wp=100/500;ws=200/500;%频率归一化[n,wn]=buttord(wp,ws,3,30);[b,a]=butter(n,wn);freqz(b,a,512,1000)要求与1》相同，但采用椭圆函数滤波器实现。wp=100/500;ws=200/500;%频率归一化[n,wn]=ellipord(wp,ws,3,30);[b,a]=ellip(n,3,30,wn);freqz(b,a,512,1000)应用最小二乘拟合的方法设计一个8阶的低通滤波器，指定的频率点与对应的幅度分别由f和m给出，并画出它希望的频率响应和实际的频率响应。f=[00.60.61];m=[1100];[b,a]=yulewalk(8,f,m);[h,w]=freqz(b,a,128);plot(f,m,w/pi,abs(h),'--');legend('ideal','yulewalkdesigned');title('comparisonoffrequencyresponsemagnitudes')实验内容5－2确定通带截止频率为500HZ，阻带截止频率为600HZ以及通带的最大衰减量为1db，阻带的最小衰减量为50db，采样频率为2000HZ的数字低通巴特沃斯滤波器的阶数，并画出其幅频特性图。实验内容5－3采用直接设计法设计一个通带截止频率为500HZ～800HZ、阻带截止频率为600HZ～700HZ、通带内波纹的最大衰减为0.1db、阻带内波纹的最小衰减为50db、采样频率为2000HZ的带阻椭圆数字滤波器，并画出其幅频响应图。实验内容5－4模拟信号xa（t）=sin(500*pi*t)+cos(50*pi*t)有系统完成xa（t）A/DH(z)D/Aya(t)的处理，采样频率为1000HZ。采用双线性变换法设计一个高通切比雪夫2型数字滤波器，要求所设计的滤波器阶数最小，以小于1db的衰减通过150hz的分量，以至少40db抑制100hz分量。滤波器应有单调的通带和等波纹的阻带。求此滤波器并画出幅频响应图（单位db）产生上述信号xa（t）的150个样本点，用图形文件在时域和频域上分别表示该信号。然后通过上述所设计的滤波器对（2）所产生的信号滤波，并分别在时域和频域上结合滤波后的输出信号进行结果分析。3．实验问题回答（1）完成例5－1的程序并观察IIR滤波器的特点。（2）在实验内容5-2、5-3和5－4中按要求完成程序，并回答提出的问题。四．实验报告要求按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。回答实验中提出的问题。总结本次实验结果，按照实验报告格式要求，书写实验报告。五．实验设备PC机，MATLAB软件实验六：FIR数字滤波器的设计一．实验目的1．加深对FIR滤波器的理解2．学会用窗函数法设计满足响应要求的FIR数字滤波器二．实验相关知识准备1．基于窗函数设计具有标准频率响应的FIR滤波器常用的窗函数矩形序列：w=boxcar(n)可产生一长度为n的矩形窗函数三角窗函数：w=triang(n)可产生n点的三角窗函数bartlett窗：w=Bartlett(n)hamming窗：w=hamming(n)hanning窗：w=hanning(n)blackman窗：w=blackman(n)chebyshev窗：w=chebwin(n,r)可产生n点的chebyshev窗函数，其傅立叶变换后的旁瓣波纹低于主瓣rdb。注意，当n为偶数时，窗函数的长度为n+1kaiser窗：w=Kaiser(n,beta)可产生N点的kaiser窗函数，其中beta为影响窗函数旁瓣的β参数。若阻带最小衰减为As。窗函数选择原则具有较低的旁瓣幅度，尤其是第一旁瓣幅度。旁瓣幅度下降数率要大，以利增加阻带衰减。主瓣的宽度要窄，以获得较陡的过渡带。通常上述几点很难同时满足，当选用主瓣宽度较窄时，虽然得到较陡的过渡带，但通带和阻带的波纹明显增加；当选用最小的旁瓣幅度时，虽然得到匀滑的幅度相应和较小的阻带波动，但过渡带加宽。因此，实际选用的窗函数往往是它们的折中。在主瓣宽度达到一定要求的条件下，适当牺牲主瓣宽度来换取旁瓣波动的减少。3)调用的函数fir1fir1函数实现加窗线性相位FIR数字滤波器设计，可设计出标准的低通、带通、高通和带阻滤波器。b=fir1(n,Wn)可设计n阶低通FIR滤波器，这里Wn为规一化截止频率，0≤Wn≤1,Wn=1相应于0.5fs。所设计的滤波器为加hamming窗的线性相位滤波器。当Wn=[W1W2]时，fir1可得到带通滤波器，其通带为W1<W<W2。b=fir1(n,Wn,’ftype’)当ftype=high时，设计高通FIR滤波器；当ftype=stop时，设计带阻滤波器。其中n必须为偶数，当输入的阶次n为奇数时，fir1函数会自动将阶数加1。b=fir1(n,Wn,window)window长度为n+1，如不指定window参数，则fir1函数采用hamming窗。b=fir1(n,Wn,‘ftype’,window)基于频率采样法和窗函数设计具有任意频率响应的FIR滤波器fir2函数用于设计具有任意频率响应的加窗数字FIR滤波器，其基本思想是基于对滤波器的幅度响应进行频率采样。b=fir(n,f,m)可设计出一n阶的FIR滤波器，其滤波器的频率特性由矢量f和m决定，f为频率采样点矢量，m为与之对应的幅度矢量，有关它们的约定可参见yulewalk函数。b=fir2(n,f,m,window)可设计具有任意频率响应形状的加窗线性相位FIR数字滤波器。其幅频特性由频率点向量f和幅度值向量m给出，0≤f≤1，要求f为单增向量，而且从0开始，以1结束，1表示数字频率w=π。m与f等长度，m(k)表示频点f(k)的幅频响应值。Plot(f,m)可画出期望逼近的幅频响应曲线。如省略window，则自动取为hamming窗。b=fir2(n,f,m,npt)格式中，所希望的频率响应被等间隔地插值，长度变为npt点，默认为512点。相应的b=fir2(n,m,npt,window)格式中，可指定窗函数。b=fir2(n,f,m,npt,lab)格式中，当所指定的f中，如果有2个相继的值相同（例如在转折点的频率），可利用参数lab指定fir2在这2个频率点间插入必要的点数，以保证滤波器的频率响应能平滑过渡且具有必要的陡度，默认值为25点。相应的b=fir2(n,f,m,npt,lab,window)格式中，可指定窗函数。三．实验内容1．例6－1设计一个34阶的高通FIR滤波器，截止频率为0.48rad，并使用具有30db波纹的chebyshev窗。window=chebwin(35,30);b=fir1(34,0.48,'high',window);freqz(b,1,512)设计一个30阶的低通滤波器，使之与期望频率特性相近。f=[00.60.61];m=[1100];b=fir2(30,f,m);[h,w]=freqz(b,1,128);plot(f,m,w/pi,abs(h))2．实验内容6－2设计一个阻带为0.4~0.7，阶数为38，窗函数为切比雪夫窗的带阻滤波器，并与窗函数为默认的hamming窗时的设计结果相比较，要求用图形文件分析比较结果并给出结论。实验内容6－3设计一个具有指定幅频响应的多带FIR滤波器，并与期望的幅频响应的结果进行比较，用图形文件说明比较结果。实验内容6－4设计一个最小阶次的低通FIR数字滤波器，性能指标为：通带0Hz～1500Hz,阻带截止频率2000Hz,通带波动1％，采样频率为8000Hz.用图形文件表示设计出的FIR低通滤波器的幅频响应曲线。3．实验问题回答（1）完成例6－1的程序并观察FIR滤波器的特点。（2）在实验内容6-2、6-3和6－4中按要求完成程序，并回答提出的问题。四．实验报告要求1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。回答实验中提出的问题。总结本次实验结果，按照实验报告格式要求，书写实验报告。五．实验设备PC机，MATLAB软件附录AMATLAB系统的常用概念1、命令窗口在Windows2000下启动MATLAB系统后，Windows2000的工作平台上会弹出一个窗口，如下图所示，这个窗口称为MATLAB的命令窗口。MATLAB的命令窗口是用户与MATLAB解释器进行通信的工作环境，提示符‘>>’表示MATLAB解释器正等待用户输入命令。所有的MATLAB命令、MATLAB函数，以及MATLAB程序都要在这个窗口下运行。在命令窗口，用户可以发出MATLAB命令。例如，为了生成一个3*3的矩阵，可以在提示符下，键入如下的命令：A=[123;456;789]方括号命令表示矩阵，空格或逗号将每行的元素分开，而分号将矩阵的各行数值分开。再键入Enter后，MATLAB将回显如下的矩阵：A=123456789为了求该矩阵的逆矩阵，则只要键入命令？B=inv(A);MATLAB就将计算出相应的结果。如果不想在命令窗口中显示计算结果，只要如上所示，在该命令后多键入一个分号即可。此时，MATLAB系统只完成该命令所要求的计算任务，其计算结果不回显。这项功能在程序设计中是非常必要的。MATLAB系统也可以说是一种新的语言，该语言十分容易掌握，其结构非常类似于数学式子的书写格式，用户花上很少的时间即可掌握MATLAB的大部分命令。2、图形窗口MATLAB系统的强大功能之一是其优秀的图形功能。对于任何作图命令，MATLAB将打开另一个窗口来绘制输出图形，这样的窗口在MATLAB系统中被成为图形窗口。在同一个图形窗口中，可以绘制多个图形，也可以生成多个图形窗口，并选择其中的一个图形窗口，在其中绘制图形。生成图形窗口的方法比较多，在没有图形窗口存在时，每个绘图函数都能自动生成一个图形窗口；也可以用figure命令生成一个新的图形窗口；还可以用命令窗口的File菜单的New子菜单的Figure项来打开一个新的图形窗口。3、搜索路径MATLAB管理着一条搜索路径，它在搜索路径下寻找与命令相关的函数文件。例如，如果在MATLAB提示符下输入example,MATLAB解释器将按照下面的步骤来处理这条字符串：检查example是不是一个变量；如果不是，检查example是不是一个内部函数；如果不是，检查在当前文件夹下是否存在名为example.mex,example.dll或example.m的文件。MEX文件是MATLAB的执行文件，将优先执行；如果不存在，检查在MATLAB的搜索路径的目录下是否存在名为MEX,example.mex,example.dll,或example.m的文件。MEX文件优先执行。用户可以打开路径浏览器（PathBrowse）查看MATLAB系统的当前搜索路径，也可以在其中加入自己的路径。4、文件类型在MATLAB系统中，根据功能可将MATLAB系统所使用的外部文件分为几类，并用不同的扩展名作为其标识，我们用的主要是M文件。M文件以字母m为其扩展名，例如startup.m。一般说来，M文件是以ASCII码文本文件，可以用任何文本编辑器进行编辑。在MATLAB系统中，有两类M文件。一类称为程序M文件，简称M文件；另一类称为函数M文件，或简称为函数，统称为M文件。M文件的内容是由符合MATLAB语法的语句构成的，函数M文件的第一行必须是以关键字function开始的函数说明语句。两类M文件的共同特征是：在MATLAB命令窗口中的命令提示符下键入文件名，来执行M文件中的所有语句规定的计算任务或完成一定的功能。它们的区别在于以下两方面：第一，程序M文件中创建的变量都是MATLAB工作空间中的变量，工作空间中的其他程序或函数可以共享，而函数M文件中创建的所有变量除了全程变量外，均局限于函数运行空间内的局部变量；第二，函数M文件的调用式中可以有输入参数和输出参数，而程序M文件则没有这种功能。5、语言语法要素MATLAB只管理一种对象——矩阵。可以使用下列的任何一种方法在MATLAB环境下创建或输入一个矩阵：显示的输入一个元素序列；用MATLAB的内部函数创建一个矩阵；在M文件中用MATLAB语句创建一个矩阵；从一个外部数据文件中装载并创建一个矩阵。在MATLAB中有两个基本概念：变量和表达式。变量由变量名表示，函数名作为特殊的变量名看待，每个变量名由一个字母后面跟随任意个字母或数字（包括下划线）组成，但MATLAB只能分辨前19个字符。MATLAB能区分组成变量名的大小写字母。MATLAB的语句则是两种形式之一：变量名=表达式或者表达式。在前一种语句形式下，MATLAB将运算的结果赋给“变量名”；而在第二种语句形式下，将运算的结果赋给MATLAB的永久变量ans，每条语句以回车符结束。一般的，运算的结果在命令窗口中显示出来。如果语句的最后一个字符是分号“；”，那末，MATLAB仅执行赋值运算，不再显示运算的结果。与C语言一样，MATLAB将字符串当作数组或矩阵处理。6、MATLAB的基本运算符矩阵运算符矩阵A的转置。A+B,A-B矩阵A与B的和与差。A*B矩阵A与B的乘法。A.*B矩阵A与B的对应元素相乘关系运算符<小于>小于或等于<=大于>=大于或等于==等于~=不等于7、特殊运算符在MATLAB的M文件中，可以加入解释行。解释行的标识符为“%”，该标识符将被作为注解内容。程序执行时，注解被忽略。方括号“[]”用于生成矩阵。特别的，语句A=[]生成空矩阵A。行分隔符“；”用于MATLAB语句后时，表示该语句的执行结果不被回显出来，这可避免显示一些不感性趣的结果。冒号“：“最主要的作用是生成向量，从下面的例子中可以看出它的使用方法：j:k生成向量[j,j+1,j+2,…,k]j:i:k生成向量[j,j+i,j+2*i,…,k]，如果j>k，则生成空矩阵A(:,j)矩阵A的第j列A(I,:)矩阵A的第I行A(j:k)向量A(j),A(j+1),…,A(k)A(:,j:k)从第j列到第k列的矩阵子块换行连接符“…”，有时一条MATLAB语句会很长，在命令窗口的一行内很可能写不下，此时只要在该语句中加入三连点，再回车即可在下一行接着写该语句。8、MATLAB的常用数学函数三角函数sin正弦函数Cos余弦函数Tan正切函数Asin反正弦函数Atan反正切函数Sinh双曲正弦函数Cosh双曲余弦函数Tanh双曲正切函数Asinh反双曲正弦函数Acosh反双曲余弦函数Atanh反双曲正切函数Acos反余弦函数初等函数Abs实数的绝对值、复数的模、字符串的ASIIC值Angle复数的幅角Sqrt方根函数Real复数的实部Imag复数的虚部Conj复共轭运算Round最邻近整数截断（四舍五入）Ceil不大于自变量的最大整数Rem不小于自变量的最小整数Exp自然指数函数（以e为底）Log自然对数函数（以e为底）log10以10为底的对数函数9、程序流控制与其他的程序设计语言一样，MATLAB语言也提供了条件语句。下面分别予以介绍。for循环语句MATLAB也有自己的for循环语句。如果要反复执行的一组语句的循环次数是已知或预定义的，就可以用for循环语句。例如：forI=1:nx(I)=0;end这条语句将向量x的前n个元素赋予零值，这里的变量n必须预先给定。注意：每一个for必须与end配对使用。while循环语句MATLAB提供有while循环语句，它的作用是在一定的逻辑条件的控制下，不断的循环执行一条或一组语句，直到逻辑条件不满足为止。While语句的一般形式是While表达式语句组endif条件语句和bread语句break语句用于退出循环体，if条件语句有两种形式，分别是if表达式语句组1else语句组2end和if表达式1语句组1elseif表达式2语句组2else语句组3end10、MATLAB的在线帮助用户可以随时利用MATLAB的在线帮助查询自己不懂得用法的函数的具体用法，例如：在命令窗口键入helpabs后的显示如下：ABSAbsolutevalueABS(X)istheabsolutevalueoftheelementsofX.WhenXiscomplex,ABS(X)isthecomplexmodulus(m