matlab里impz指令格式,华北电力大学Matla实验指导书.doc

《华北电力大学Matla实验指导书.doc》由会员分享，提供在线免费全文阅读可下载，此文档格式为doc，更多相关《华北电力大学Matla实验指导书.doc》文档请在天天文库搜索。

1、数字信号处理实验指导书(基于MATLAB)尚秋峰 宋文妙华北电力大学二○○六年十一月前 言MATLAB 是一套功能强大的工程计算及数据处理软件，广泛应用于工业，电子，医疗和建筑等众多领域。它是一种面向对象的，交互式程序设计语言，其结构完整，又具有优良的可移植性。它在矩阵运算，数字信号处理方面有强大的功能。另外，MATLAB提供了方便的绘图功能，便于用户直观地输出处理结果。本课程实验要求学生掌握用MATLAB语言仿真数字信号处理算法的基本方法，加深对教学内容的理解。说明：本书给出的为本课程基本实验，教师可根据学生情况安排内容延伸实验。目 录实验一 序列的产生和运算 1实验二 因果性数字系统的时域实现 5实验三 离散傅里叶变换(DFT)及其快速算法(FFT) 10实验四 FFT的典型应用 17实验五 FIR滤波器设计 23实验六 IIR滤波器设计 29附 录 32实验一 序列的产生和运算一、。

2、实验目的：1、熟悉MATLAB环境，掌握基本编程方法。2、熟悉MATLAB中序列产生和运算的基本函数。二、实验内容：1、MATLAB入门(1)了解MATLAB 的工作窗口进入MATLAB环境，选择View中的各个选项，可以打开命令窗口Command Window、变量浏览器Workspace、路径浏览器Current Directory，如下图示：图1－1 Matlab 工作窗口MATLAB命令窗口，是键入指令的地方也是MATLAB显示计算结果的地方。在符号>>之后键入“1+2+3”，或者“x=1+2+3”，显示结果有什么不同？如果在上述的例子结尾加上“;”，则计算结果不会在窗口自动显示，要显示计算结果须键入该变量x。进行以上操作后打开变量浏览器，观察里面的变量名和变量的值。在MATLAB环境选择File 再选择New，即进入程序编辑/调试窗口，如下图示。存档时必须以.m名称储。

3、存。要执行 M-file 可以在命令窗口下直接键入该文件名；或是选择Debug下的Run M-file来执行M-file。如果要修改 M-file 可以选择功能表上的Open M-file ，即可搜寻要修改的 M-file，修改后再存档。尝试编写程序文件，完成以上操作。图1－2 程序编辑/调试窗口练习使用路径浏览器打开特定目录下的M-file。M-file还可以用来定义函数，然后储存起来，就可以和那些内建的函数(如sin, cos,log等)一样的自由使用。举例来说，我们可以定义一函数cirarea来计算圆的面积，以下的 M-file: cirarea.m就是定义这个函数的：% M-file function, cirarea.m % Calculate the area of a circle with raduis r % r can be a scalar or an array f。

4、unction c=cirarea(r) c=pi*r.^2;注意，M-file定义的函数语法上有一些规定： Ø 第一行指令以function这个字做为起头，接着是输出的变量，等号，函数名称，输入的变量是放在括号之内。function out1=userfun(in1)，这行的out1是输出的变量，userfun是函数名称，in1是输入的变量。function [out1, out2]=serfun(in1, in2) 如果输出变量 [out1,out2]和输入变数(in1, in2)不只一个时，则在输出变量部份须加上[ ]。 Ø 上述的输入变量是在调用函数时输入的，而输出的变量即是函数的返回值。 Ø 函数名称的取法的规定与一般变量相同。 Ø 在定义函数之前，最好加上注解行来说明这个函数的特色及如何使用，这样，使用指令如help cirarea，该函数的注解行会出现在指令视窗。尝试编写函数。

5、文件，并在你编写的程序文件中应用此函数。MATLAB会将绘图结果展示在另一个视窗，称为MATLAB Figure Windows，如下图示。如果你看不到此视窗，可以进入Windows再选择Figure。在图形窗口中，可以利用Edit菜单中的选项来改变显示效果以及拷贝图形，尝试这些操作。图1－3 MATLAB 图形视窗(2)寻求帮助在MATLAB系统中相关的线上(on-line)求助方式有三种： ² 利用help指令，如果你已知要找的主题为何的话，直接键入help 。所以即使身旁没有使用手册，也可以使用help指令查询不熟悉的指令或是题材之用法，例如help sqrt, help topic。 ² 利用lookfor指令，它可以从你键入的关键字(key-word)(即使这个关键字并不是MATLAB的指令)列出所有相关的主题，例如lookfor cosine, lookfor sine。 ² 。

6、利用指令视窗的功能选单中的Help，从中选取Table of Contents(目录)或是Index(索引)。² 练习以上寻求帮助的方法。2、波形产生学习附录中的有关知识，完成下列练习。练习一：建立一个MATLAB函数impseq.m，用来生成迟延的单位脉冲序列。其输入参数为：序列起始位置ns，序列终止位置nf，脉冲位置np。练习二：编写MATlAB程序，以产生下图所示的波形。并将序列绘制出来。练习三：设x＝[2 3 4 5]，求将它延拓6个周期所得的序列。练习四：设x＝[2 3 4 5]，实现以下运算并绘制波形练习五： 编写程序实现下列函数，并绘制波形以 t以t=0.01(n=0:N-1)进行取样，N=64三、思考题：算术运算符*和．*之间的区别是什么？ 实验二 因果性数字系统的时域实现一、 实验目的编制实现下列差分方程的程序: y(n)=∑bkx(n-k)+∑aky(n-k)二、 实验。

7、内容1、编制nonrec.m函数文件,实现FIR滤波y(n)=h(n)*x(n)。这里给定h(n)=R8(n)，x(n)= nR16(n)，求y(n)。nonrec.m函数文件:function y=nonrec(x,h)x=[x,zeros(1,length(h)-1)]; %补零w=zeros(1,length(h));for i=1:length(x) for j=length(h):-1:2w(j)=w(j-1); %得到每一延时单元输出变量 end w(1)=x(i); y(i)=w*h.’; %h与w对应相乘end主程序文件:x=0:15;h=ones(1,8);y=nonrec(x,h);n=0:22;stem(n,y);图2－1 卷积运算实现FIR滤波分析：线性卷积y(n)=x(n)*h(n)的长度为16+8-1=23，可利用y(n)=∑h(m)x(n-m)直接计算得 n(。

8、n+1)/2, n≤7y(n)= 4(2n-7), 8≤n≤15 (n+8)(23-n)/2, 16≤n≤22即 y=[ 0 1 3 6 10 15 21 28 36 44 52 60 68 76 84 92 84 75 54 42 29 15] ,与曲线相符。2、编制rec.m函数文件,实现纯递归IIR滤波 y(n)=x(n)+∑aky(n-k)。 这里给定a1=2rcosw0, a2=-r2, r=0.95, w0=π/8, 求单位抽样响应h(n).rec.m函数文件:function y=rec(x,a,n)x=[x,zeros(1,n-length(x))]; %补零到所需长度sum=0;w=zeros(1,length(a));for i=1:n y(i)=sum+x(i); % 递归 for j=length(a):-1:2 w(j)=w(j-1); %延时 end w(1)=。

9、y(i); sum=w*a';end主程序文件:x=[1];a=[2*0.95*cos(pi/8),-0.95^2];h=rec(x,a,75); %取h(n)的长度为75点 n=0:74;stem(n,h);图2－2 纯递规IIR滤波分析计算:由题意, a1=2*0.95*cos(π/8), a2=-0.952, 所以,得到系统函数 H(z)=1/[1-1.9cos(π/8)z-1+0.952z-2],做逆Z变换得 h(n)=0.95ncos(πn/8)+ctg(π/8)*0.95nsin(πn/8),利用MATLAB直接画h(n), 即,使用下列语句 n=0:74;h=0.95.^n.*cos(pi.*n./8)+cot(pi/8).*(0.95.^n).*sin(pi.*n./8);stem(n,h); 得到如图2－3的曲线图2－3 此曲线与用rec函数所画曲线完全一致, 可见, 用。

10、MATLAB编制的FIR滤波程序与理论计算所得结果是相同的.3、用nonrec.m和rec.m函数编制dfilter.m 函数文件, 构成完整的一般IIR滤波程序, 并完成下列信号滤波:x(n)=cos(2πn/5)R64(n) 这里给定系统函数 H(z)=(1-2z-1+z-2)/(1-0.5z-1+0.5z-2) 求y(n).dfilter.m函数文件:function y=dfilter(x,b,a,n)y1=nonrec(x,b); y=rec(y1,a,n);主程序文件:n=0:63;x=cos(2*pi/5*n);b=[1,-2,1]; %由H(z)得到系数a,ba=[0.5,-0.5];y=dfilter(x,b,a,64); %取y(n)的长度为64点stem(n,y);图2－44．用help查看内部函数filter.m, 了解调用格式, 重做3, 并与你编写的dfilte。

11、r.m结果进行比较.用help可以看到内部函数为Y = FILTER(B,A,X), 且有 “The filter is a Direct Form II Transposed implementation of the standard difference equation”: a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) -a(2)*y(n-1) - ... - a(na+1)*y(n-na) 因此, 调用内部函数filter时, 要对原系数a做适当变化: a1(1)=1, a1(i)=-a(i-1), i>1.n=0:63;x=cos(2*pi/5*n);b=[1,-2,1];a=[0.5,-0.5];y=dfilter(x,b,a,64); %调用自己编的dfilter函数a1=[1,-0.5,0.5。

12、]; %a变为a1, 用于调用内部函数filtery1=filter(b,a1,x);subplot(2,1,1);stem(n,y);subplot(2,1,2);stem(n,y1);图2－5三、 思考题1、 内部函数filter.m的调用格式是什么?与你编写的dfilter.m调用格式是否一致?差异在何处?2、 实验中如何合理地选取单位抽样响应h(n)的点数?实验内容及要求中第三项的物理概念是什么?给定的H(z)是具有什么性质的滤波器?3、 为实现各实验内容中的滤波器，你编写程序时采用的什么结构？4、 试利用Matlab中的卷积函数conv完成实验内容1,并与本书的例程进行比较。实验三 离散傅里叶变换(DFT)及其快速算法(FFT)一、实验目的1、学习DFT和FFT的原理及编程实现方法。2、熟悉MATLAB编程的特点。二、实验内容1、用三种不同的DFT程序计算x(n)=R8(n)的傅。

13、里叶变换X (ejw)，并比较三种程序计算机运行时间。(1)用for loop 语句的M函数文件dft1.m，用循环变量逐点计算X(k)；(2)编写用MATLAB矩阵运算的M函数文件dft2.m, 完成下列矩阵运算；(3)调用FFT库函数，直接计算X(k)；(4)分别利用上述三种不同方式编写的DFT程序计算序列x(n)的傅立叶变换X(ejw)，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。(一)实验例程M函数文件如下：dft1.m:function[Am,pha]=dft1(x)N=length(x);w=exp(-j*2*pi/N);for k=1:N sum=0; for n=1:N sum=sum+x(n)*w^((k-1)*(n-1)); end Am(k)=abs(sum); pha(k)=angle(sum);enddft2.m:function[Am,pha]=。

14、dft2(x)N=length(x);n=[0:N-1];k=[0:N-1];w=exp(-j*2*pi/N);nk=n'*k;wnk=w.^(nk);Xk=x*wnk;Am=abs(Xk); pha=angle(Xk);dft3.m:function[Am,pha]=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk);源程序及运行结果：(1)x=[ones(1,8),zeros(1,248)];t=cputime;[Am1,pha1]=dft1(x);t1=cputime-tn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(1)subplot(2,1,1), plot(w,Am1,'b'); grid;title('Magnitude part');xlabel('frequency in radians');。

15、ylabel('|X(exp(jw))|');subplot(2,1,2), plot(w,pha1,'r'); grid;title('Phase Part');xlabel('frequency in radians');ylabel('argX[exp(jw)]/radians');图3－1 DFT1记录运行时间：t1=(2)x=[ones(1,8),zeros(1,248)];t=cputime;[Am2,pha2]=dft2(x);t2=cputime-tn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(2)subplot(2,1,1), plot(w,Am2,'b'); grid;title('Magnitude part');xlabel('frequency in radians');ylabel('|X(exp(jw))|')。

16、;subplot(2,1,2), plot(w,pha2,'r'); grid;title('Phase Part');xlabel('frequency in radians');ylabel('argX[exp(jw)]/radians');图3－2 DFT2 记录运行时间：t2=(3)x=[ones(1,8),zeros(1,248)];t=cputime;[Am3,pha3]=dft3(x);t3=cputime-t;n=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(3)subplot(2,1,1), plot(w,Am3,'b'); grid;title('Magnitude part');xlabel('frequency in radians');ylabel('|X(exp(jw))|');subplot(2,1,2), plo。

17、t(w,pha3,'r'); grid;title('Phase Part');xlabel('frequency in radians');ylabel('argX[exp(jw)]/radians');图3－2 DFT2记录运行时间：t3=从以上运行结果可以看出，调用FFT库函数直接计算X(k)速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。因此FFT算法大大提高了DFT的实用性。(二)练习改变输入信号类型和点数，重做例程中的内容，并整理总结实验结果，给出结论。2、实现序列的内插和抽取所对应的傅里叶变换。给定序列x(n)=[cos(π/36×n)+cos(1.5π/36×n)]R128(n)，做N=128点的傅里叶变换，并求x1(n)=x(4n) 和 对应的傅里叶变换(N=128点)。(一)实验例程源程序如下：N=32; %抽取32点t=0:N-1;x2=cos(pi/9*t)。

18、+cos(1.5*pi/9*t);Xk2=fft(x2,128); %做128点傅里叶变换Am2=abs(Xk2);n=[0:127];w=(2*pi/128)*n;figure(2)plot(w,Am2,'b'); title('Magnitude part');xlabel('frequency in radians');ylabel('|X(exp(jw))|');x3=zeros(1,128); %初始化数组为零N=32;for i=0:N-1 %n=4i时赋值 x3(1,4*i+1)=cos(pi/36*i)+cos(1.5*pi/36*i);endXk3=fft(x3); %做128点fftAm3=abs(Xk3);n=[0:127];w=(2*pi/128)*n;figure(3)plot(w,Am3,'b');title('Magnitude part');xlabel('。

19、frequency in radians');ylabel('|X(exp(jw))|');抽取后的离散傅里叶变换：图3－4 抽取的离散傅里叶变换内插后的离散傅里叶变换：图3－3 内插的离散傅里叶变换(二)练习仿照例程编写给定序列x(n)=[cos(π/36×n)+cos(1.5π/36×n)]R128(n)的128点傅里叶变换程序，比较这三个计算结果得到的幅频特性图，分析其差别产生的原因。三、思考题实验内容4中，分别求X1(ejww)和X2(ejww)与原X(ejww)的关系，并说明为什么抽取后能分辨出两个频率分量？有无任何条件的限制？实验四 FFT的典型应用一、实验目的1、学习FFT的应用于卷积运算和谱分析时的编程实现方法。2、熟悉MATLAB编程的特点。二、实验内容1、利用FFT实现两序列的卷积运算，并研究FFT点数与混叠的关系。给定x(n)=R16(n)，h(n)=R8(n)，用F。

20、FT和IFFT分别求线性卷积和混叠结果输出，并用函数stem(n,y)画出相应图形。(一)实验例程源程序如下：N=16;x=[0:N-1];h=ones(1,8);线性卷积：Xk1=fft(x,23); %做23点fftHk1=fft(h,23);Yk1=Xk1.*Hk1;y1=ifft(Yk1);n=0:22;figure(1)stem(n,y1)运行结果：线性卷积：图3－1 线性卷积(二)练习改写上面例程做18点fft,由程序运行结果总结前多少个点为混叠结果，后多少个点是线性卷积的结果。如果改做30点fft，又会出现什么结果？依据实验结果总结：要用DFT来做线性卷积，DFT的点数必须满足什么要求？2、研究高密度频谱与高分辨率频谱(一)实验例程设有连续信号Xa(t)=cos(2π×6.5×10^3t)+cos(2π×7×10^3t)+cos(2π×9×10^3t)以采样频率fs=32kH。

21、z对信号Xa(t)采样,分析下列三种情况的幅频特性。(1) 采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。(2) 采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。(3) 采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。源程序如下：fs=32000; %采样频率N=16; %采集16点n=0:N-1; %做16点DFTxa=cos(2*pi*6.5*10^3*n/fs)+cos(2*pi*7*10^3*n/fs)+cos(2*pi*9*10^3*n/fs);Xk1=fft (xa);Am1=abs(Xk1);n=[0:(length(xa)-1)];w=(2*pi/length(xa))*n;figure(1)plot(w,Am1,'b'); grid;title('。

22、Magnitude part');xlabel('frequency in radians');ylabel('|X(exp(jw))|');x=[xa(1:16),zeros(1,240)]; %补零到256Xk2=fft(x); %做256点DFTAm2=abs(Xk2);n=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(2)plot(w,Am2,'b'); grid;title('Magnitude part');xlabel('frequency in radians');ylabel('|X(exp(jw))|');N0=256; 。效的数组方法来解给定的问题时，应避免用For循环。两种方法得出同样的结果，而后者执行更快，更直观，要求较少的输入。5. 为了得到最大的速度，在For循环(While循环)被执行之前，应预先分配数组。九、滤。

23、波器设计1、滤波器结构IIR滤波器结构IIR滤波器可写成即可用滤波器系数和表示。当时，滤波器的阶为N。IIR滤波器的差分方程表示为实现IIR滤波器可采用三种不同的结构；直接形式、级联形式和并联形式。1．直接形式滤波器差分方程直接用延迟线、乘法器和加法器实现，例如M=4，N=4时，得到如图3所示的结构，其差分方程为这种结构称为直接形式1。将这种结构进行简化，可得到直接形式2，如图4所示。 在MATLAB中，可利用DSP工具箱函数filter实现IIR的直接形式。图3 IIR滤波器直接形式1结构(N=4) 图4 IIR滤波器直接形式2结构(N=4)2．级联形式下面先假设N为偶数，则IIR可简化其中K：N／2，每个二阶形式都可用直接形式2得以实现，例如当N=4时，其级联结构如图1．5所示。图5 IIR滤波器的级联结构(N=4)在MATLAB中，给定直接形式滤波器的系数{， }，可计算出，，这可由。

24、扩展函数dir2cas实现。然后利用扩展函数casfiltr实现滤波器的级联形式。当然，在给定级联形式时，也可得到滤波器的直接形式，这由扩展函数cas2dir实现。3．并联形式在级联形式的基础上，利用部分分式展开，可进一步得到(设M≥N)其中K=N／2，每一部分均为二阶，它们之间是并联关系。同样以M=4，N=4为例，画出其结构如图6所示。图6 IIR滤波器的并联结构在MATLAB下，我们设计了将滤波器直接形式转化为并联形式的扩展函数dir2par，其中用到了复共轭对比较函数cplxcomp。这样我们就可利用parfiltr函数实现滤波器的并联形式。最后还设计出了将并联形式转换成直接形式的扩展函数par2dir。图7 FIR滤波器的直接结构形式其中K=M/2.以M=7为例，其结构如图8所示。图8 FIR滤波器级联形式(M=7)在MATLAB中，可采用dir2cas和cas2dir函数，在直接。

25、形式和级联形式之间转换，滤波器实现可通过casfiltr函数实现。3．线性相位形式线性相位条件其中b=O或b=±p／2，a为常数。对因果FIR滤波器，设冲激响应处于[O，M-1]区间上，则上述条件就表示下列对称性(对称冲激响应): (反对称冲激响应): 以M=7或M=6为例，对称型的结构如图1．9所示。在MATLAB中，实现FIR线性相位形式的函数等同于直接形式，即采用filter函数。图9 FIR滤波器线性相位形式(M=7或M=6)4．频率取样形式冲激响应h(n)的M点DFT为H(k)(O≤k≤M一1)，则有利用内插公式可得:这表明在这种结构中采用了离散傅里叶变换H(k)，而不是冲激响应h(n)。这种结构分成两部分：第一部分由(构成，第二部分由M个并联而成。在MATLAB中，给定h(n)或H(k)，则可借助于扩展函数dir2fs得到频率取样形式，2、滤波器分析与实现函数1．abs 功能：。

26、求绝对值(幅值)。 格式： y=abs(x) 说明：y=abs(x)用于计算x的绝对值。当x为复数时，得到的是复数模(幅值)，即当x为字符串时，abs(x)得到字符串的各个字符的ASCII码，例如x=’123 ’，则abs(x)得到49 50 51。 2．angle功能：求相角。格式： p=angle(h) ，说明： p=angle(h)用于求取复矢量或复矩阵h的相角(以弧度为单位)，相角介于-π和π之间。例如对复数h可用两种方法表示 则幅值m和相角p可由x+i y格式求出m=abs(h)p=angle(h)当然，由m和p也可求取x+i y格式 h=m．*exp(i*p) x=real(h) y=imag(h)3．conv 功能：求卷积。 格式： C=CONV(a，b) 说明： conv(a，b)用于求取矢量a和b的卷积，即其中N为矢量a和b的最大长度。例如，当a=[1 2 3]，b：[4。

27、 5 6]时，则 c=conv(a，b) c= 4 13 28 27 184.filter功能：利用递归滤波器(IIR)或非递归滤波器(FIR)对数据进行滤波。格式： y=filter(b，a，x) [y，zf]=filter(b，a，x) y=filter(b，a，x，zi)说明： filter采用数字滤波器对数据进行滤波，其实现采用移位直接型Ⅱ结构，因而适用于IIR和FIR滤波器.滤波器的z域表示为这里a(1)=1，如果输入的a(1)≠1，则MATLAB将自动进行归一化系数的操作；如果a(1)=O，则给出出错信息。y=filter(b，a，x)利用给定的矢量a和b，对x中的数据进行滤波，结果放入y矢量，长度取max(，)。[y，zf]=filter(b，a，x)除得到结果矢量y外，还得到x的最终状态矢量zf。y=filter(b，a，x，zi)可在zi中指定x的初始状态。 5. fil。

28、tic 功能：为filter函数选择初始条件。 格式： z=filtic(b，a，y，x) z=filtic(b，a，y) 说明： filtic函数在给定过去输出y和输入x的前提下，为移位直接形式Ⅱ滤波器实现中的延迟找到初始条件z。矢量x和y分别表示过去的输入和输出x={x(-1)，x(-2)，…，x(-nb)，…}y={y(-1)，y(-2)，…，y(-na)，…}其中nb=length(b)-l(分子阶数)，na=length(a)-1(分母阶数)。矢量x和y的长度应分别等同于nb和na，如果长度不够，则补O；如果长度超出，则截断处理，使x和y的长度正好为na和nb。输出矢量z的长度为max(na，nb)。 z=filtic(b，a，y，x)可得到给定输入x和输出y时的初始状态；z=filtic(b，a，y)则假定输入x=O。6.freqs 功能：模拟滤波器的频率响应。 格式: h=f。

29、reqs(b，a，w) l [h，w]=freqs(b，a) I [h，w]=freqs(b，a，n) 说明： freqs用于计算由矢量a和b构成的模拟滤波器 : 的复频响应。 H=freqs(b，a，w)用于计算模拟滤波器的复频响应，其中实矢量w用于指定频率值，即freqs沿虚轴计算频率响应。 [h，w]=freqs(b，a)自动设定200个频率点来计算频率响应，这200个频率值记录在w中。[h，w]=freqs(b，a，n)设定n个频率点计算频率响应。 不带输出变量的freqs函数，将在当前图形窗口中绘制出幅频和相频曲线。 7.freqz功能：数字滤波器的频率响应。 格式： [h，w]=freqz(b，a，n) [h，f]=freqz(b，a，n，Fs) [h，w]=freqz(b，a，n，’whole’) [h，f]=freqz(b，a，n，’whole’，Fs) h=freqz(b。

30、，a，w) h=freqz(b，a，f，Fs) freqz(b，a)说明： freqz用于计算由矢量a和b构成的数字滤波器的复频响应。[h，w]=freqz(b，a，n)可得到数字滤波器的n点的复频响应，这n个点均匀地分布在上半单位圆(即0～p)，并将这n点频率记录在w中，相应的频率响应记录在h中。至于n值的选择没有太多的限制，只要n>0的整数，但最好能选取2的幂次方，这样就可采用 FFT算法进行快速计算。如果缺省，则n=512。[h，f]=freqz(b，a，n，Fs)允许指定采样终止频率Fs(以Hz为单位)，也即在O～Fs／2频率范围内选取n个频率点(记录在f中)，并计算相应的频率响应h。[h，w]=freqz(b，a，n，’whole’)表示在O～2p之间均匀选取n个点计算频率响应。[h，f]=freqz(b，a，n，’whole’，Fs)则在O～Fs之间均匀选取n个点计算频率。

31、响应。h=freqz(b，a，w)计算在矢量w中指定的频率处的频率响应，但必须注意，指定的频率必须介于O和2p之间。h=freqz(b, a，f，Fs)计算在矢量f中指定的频率处的频率响应，但指定频率必须介于 O和Fs之间。8.impz 功能：数字滤波器的冲激响应。 格式： [h，t]=impz(b，a) [h，t]=impz(b，a，n) [h，t]=impz(b，a，n, Fs) impz(b，a) 说明： 由矢量a和b构成数字滤波器，即[h，t]=impz(b，a)计算出滤波器的冲激响应h，取样点数n由impz函数自动选取，并记录在矢量t中。[h，t]=impz(b，a,n)可由用户指定取样点数或取样时刻。当n为标量时，即在0～n-1时刻计算冲激响应，0时刻表示滤波器的起始点；当n为矢量(其值应为整数)，则表示t=n,即在这些指定的时刻计算冲激响应。9.residuez 功能：Z变换。

32、部分分式展开或留数计算。 格式: [r，p，k]=residuez(b，a) [b，a]=residuez(r，p，k)residuez函数可将以多项式之比表示的离散时间系统转化成部分分式展开或留数形式的系统，而且可进行相反转换。[r，p，k]=residuez(b，a)可将以多项式之比b(z)／a(z)转化成留数(r)、极点(p)和直接项(k)的部分分式展开。设 如果不存在重根，且m>n一1，则有因此函数[r，p，k]=residuez(b，a)得到列矢量r为留数，列矢量p为极点位置，行矢量为直接项。极点数为n=length(a)-1=length(r)-length(p)如果length(b)

33、含 [b，a]=residuez(r，p，k)可将部分分式转化成多项式系数a和b。注意，residuez函数与MATLAB环境中的residue函数非常类似，只不过，residue函数适用于连续系统，而residuez函数适用于离散系统。 9.butter功能：Butterworth(巴特沃思)模拟和数字滤波器设计。格式： [b，a]=butter(n，Wn) [b，a]=butter(n，Wn，’ftypet’) [b，a]=butter(n，Wn，’s’) [b，a]=butter(n，Wn，’ftype’，’s ’) [z，p，k]=butter(…) [A，B，C，D]=butter(…)说明：butter函数可设计低通、带通、高通和带阻的数字和模拟滤波器，其特性为使通带内的幅度响应最大限度地平坦，这会损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数，但在期。

34、望下降斜度大的场合，应使用椭圆和Chebyshev(切比雪夫)滤波器。butter函数可设计出数字域和模拟域的Butterworth滤波器。(1)数字域[b，a]=butter(n，Wn)可设计出截止频率为Wn的n阶低通Butterworth滤波器，其滤波器为截止频率是滤波器幅度响应下降至1／处的频率，但与Besself函数不同，Wn∈[O，1]，其中1相应于O.5fs。(取样颇率,即Nyquist频率)．当Wn=[Wl W2](Wl

35、阶的选择1．buttord功能：Butterworth滤波器阶的选择。格式[n，Wn]= buttord(wp，Ws，Rp，Rs) [n，Wn]=buttord(Wp，ws，Rp，Rs，’s’)buttord可在给定滤波器性能情况下，选择模拟或数字Butterworth滤波器最小的阶。其中，Wp和Ws分别是通带和阻带的拐角频率(截止频率)，其值为0≤wp(或Ws)≤1，当其值为1时表示0.5fs。Rp和Rs分别是通带和阻带区的波纹系数。(1)数字域[n,Wn]=butter(Wp，Ws，Rp，Rs)可得到数字Butterworth滤波器的最小阶n，并使在通带(0, Wp)内波纹系数小于Rp，在阻带(Ws，1)内衰减系数大于Rs。buttord还得到了截止频率Wn，这样利用butter函数可产生满足指定性能的滤波器。利用buttord函数，还可以得到高通、带通和带阻滤波器的阶。当Wp>。

36、Ws时，这时为高通滤波器；当Wp，Ws为二元矢量时，若Wp

37、系数包含在b中，这可表示成 b(z)=b(1)+b(2)z-1+…+b(n+1)z-n这是一个截止频率为Wn的Hamming(汉明)加窗线性相位滤波器，0≤Wn≤1，Wn=l相应于O.5fs 当Wn=[Wl W2]时，firl函数可得到带通滤波器，其通带为w1

38、w长度为n+1。如果不指定Window参数，则firl函数采用Hamming窗。 b=firl(n，Wn，’ftype’，Window)可利用ftype和Window参数，设计各种加窗的滤波器。 由firl函数设计的FIR滤波器的群延迟为n／2。5、变换1．fft功能：一维快速傅里叶变换(FFT)。格式： y=fft(x) y=fft(x，n)说明： fft函数用于计算矢量或矩阵的离散傅里叶变换，这可通过实现，其中N＝length(x)，。y＝fft(x)为利用FFT算法计算矢量x的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。当x的长度为2的幂次方时，则fft函数采用基2的FFT算法，否则采用稍慢的混合基算法。y＝fft(x，n)采用n点FFT。当x的长度小于n时，fft函数在x的尾部补零，以构成n点数据；当x的长度大于n时，fft函数会截断序列x。当x为矩阵时，fft函数按类。

39、似的方式处理列长度。2．IFFT功能：一维逆快速傅里叶变换(IFFT)。格式： y=ifft(x) y=ifft(x，n)说明： ifft函数用于计算矢量或矩阵的逆傅里叶变换，即其中N=length(x)，y=ifft(x)用于计算矢量x的IFFT。当x为矩阵时，计算所得的y为矩阵x中每一列的IFFT。y=ifft(x，n)采用n点IFFT。当length(x)n时，将x截断，使length(x)=n。ifft函数是fft函数的逆，其相应的M文件中采用的算法类似。6．数字滤波器结构MATLAB的扩展函数有11种。1．变直接形式为级联形式 dir2cas．m function[bO，B，A]=dir2cas(b，a)； ％DIRECT-form to CASCADE-form conversion(cplxpair version) ％============================。

40、========================== ％[bO，B，A]=dir2cas(b，a) ％b0一gain coefficient ％B=K by 3 matrix of real coefficients containing bk's ％A=K by 3 matrix 0f real coefficients containing ak's ％b=nomirator polynomial coefficients of DIRECT form ％a=denominator polynomial coefficients of DIRECT form ％compute gain coefficient b0 b0=b(1)；b=b／bO； a0=a(1)；a=a／a0； bO=b0／a0； ％ M=length(b)；N＝length(a)； if N>M b=[b ze。

41、ros(1，N-M)]； elseif M>N a=[a zeros(1，M-N)]：N=M； elseNM=O： end ％ K=floor(N／2)；B=zeros(K，3)；A=zeros(K，3)； if K*2==N； b=[b O]； a=[a O]； end ％ broots=cplxpair(roots(b))； aroots=cplxpair(roots(a))； for i=1：2：2*K Brow=broots(i：1：i+1，：)； Brow=real(poly(Brow))； B(fix(i+1)／2，：)=Brow； Arow=aroots(i：1：i+1，：)； Arow=real(poly(Arow))； A(fix(i+1)／2，：)=Arow； end2、bilinearBilinear transformation with optional f。

42、requency prewarping. [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer function specified by Z, P, and K to a z-transform discrete equivalent obtained from the bilinear transformation: H(z) = H(s) |s = 2*Fs*(z-1)/(z+1) where column vectors Z and P specify the zeros and poles, scalar K specifies the gain, and Fs is the sample frequency in Hz. [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where N。

43、UM and DEN are row vectors containing numerator and denominator transfer function coefficients, NUM(s)/DEN(s), in descending powers of s, transforms to z-transform coefficients NUMd(z)/DENd(z). [Ad,Bd,Cd,Dd] = BILINEAR(A,B,C,D,Fs) is a state-space version. Each of the above three forms of BILINEAR accepts an optional additional input argument that specifies prewarping. For example, [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs,Fp) applies prewarping before the bilinear transformation so that the frequency responses before and after mapping match exactly at frequency point Fp (match point Fp is specified in Hz).。