离散信号和系统分析的MATLAB实现

1、1.9 离散信号和系统分析的MATLAB实现1.9.1 利用MATLAB产生离散信号 用MATLAB表示一离散序列xk时，可用两个向量来表示。其中一个向量表示自变量k的取值范围，另一个向量表示序列xk的值。例如序列xk=2,1,1,-1,3,0,2可用MATLAB表示为 K=-2:4;x=2,1,1,-1,3,0,2可用stem(k,f)画出序列波形。当序列是从k=0开始时，可以只用一个向量x来表示序列。由于计算机内寸的限制，MATLAB无法表示一个无穷长的序列。例1-38 利用MATLAB计算单位脉冲序范围内各点的取值。解： %progran 1_1 产生单位脉冲序列 Ks=-4;ke=4;

2、n=2;K=ks:ke;X=(k-n)=0;Stem(k,x):xlabel(k);程序产生的序列波形如图1-49所示。例1-39 利用MATLAB画出信号Xk=10sin(0.02)+nk, 的波形。其中nk表示为均值为0方差为1的Gauss分布随机信号。解：MALAB提供了两个产生（伪）随机序列的函数。Rand(1,N)产生1行N列的0，1均匀分布随机数。Randn(1,N)产生1行N列均值为0方差为1的Gauss分布随机数。%program 1_2 产生受噪声干扰的正弦信号N=100;k=0:N;X=10*sin(0.02*pi*k)+randn(1,N+1);Plot(k,x);Xla

3、bel(k);Ylabel(xk);程序产生序列如图1-50所示。1.9.2 离散卷积的计算 离散卷积是数字信号处理中的一个基本运算，MTLAB提供的计算两个离散序列卷积的函数是conv，其调用方式为 y=conv(x,h)其中调用参数x,h为卷积运算所需的两个序列，返回值y是卷积结果。MATLAB函数conv的返回值y中只有卷积的结果，没有y的取值范围。由离散序列卷积的性质可知，当序列x和h的起始点都为k=0时，y的取值范围为k=0至length(x)+length(h)-2。当序列x或（和）h的起始点不是k=0时，由例1-3知，y的非零值范围可根据例1-4的结论进行计算。例1-40试用MA

4、TLAB函数conv计算例1-2中序列的卷积。解：program 1_3 计算离散卷积x=-0.5,0,0.5,1; %序列x的值kx=-1:2; %序列x的取值范围h=1,1,1;kh=-2:0;y=conv(x,h); %计算卷积k=kx(1)+kh(1):kx(end)+kh(end); %计算y的取值范围stem(k,y);xlabel(k);ylabel(y);程序的运行结果如图1-51所示。1.9.3 离散LTI系统响应MATLAB求解许多离散LTI都可用如下的线性常系数的差分方程描述 (1-151)其中xk、yk分别系统的输入和输出。在已知差分方程的N个初始状态yk,和输入xk,

5、就可由下式迭代计算出系统的输出。 yk=-利用MATLAB提供的filter函数，可方便地计算出上述差分方程的零状态响应。filter函数调用形式为 y=filter(b,a,x)其中 a=a0,a1,aN, b=bo,b1,bM分别表示差分方程系数。X表示输入序列，y表示输出序列。输出序列的长度和序列相同。例1-40 试用M=9点滑动平均系统 yk= 处理例1-39中的受噪声干扰的正弦信号。 解： 由式（1-151）可知，M点滑动平均系统可看成是N=0的差分方程。在调用filter函数时，调用参数a=1,b为有M个元素的向量，b中每个元素的值均为1/M。 %program 1_4 滑动平均去

6、噪 M=9; N=100;k=0:N; s=10*sin(0.02*pi*k); x=s+randn(1,N+1); b=ones(M,1)/M; a=1; y=filter(b,a,x); plot(k,y,s,:); xlabel(k); ylabel（幅度） legend(yk,sk);程序的运行结果如图1-52所示。图中的虚线表示未受噪声干扰的原信号sk,实线为9点滑动的响应yk。比较图1-50的信号可以看出，系统滤出了受干扰信号中的大部分噪声，输出信号相对输入信号有4个样本的延迟。1.9.4 利用MATLAB计算DTET当序列的DTET可写成的有理多项式时，可用MATLAB信号处理工

7、具箱提供的freqz函数计算DTFT的抽样值。另外，可用MATLAB提供的abs、angle、real、imag等基本函数计算 DTFT的幅度、相位、实部、虚部。若X（）可表示为 X（）=则freqz的调用形式为 X=freqz(b,a,w) (1-153)式（1-153）中的b和 a分别是表示式（1-152）中分子多项式和分母多项式系数的向量，即 b=b0,b1,bM a=a0,a1,aNw为抽样的频率点，在以式（1-153）形式调用freqz函数时，向量w的长度至少为2。返回值X就是DTFT在抽样点w上的值。注意一般情况下，函数freqz的返回值X是复数。例1-41 已知离散系统的H（z）

8、为 H（z）= 试画出该系统的幅度响应。解： %program 1_5 计算系统幅度响应 b1=0.5009 -1.0019 0.5009; b2=0.5320 1.0640 0.5320; a1=1.0000 -0.8519 0.4167; a2=1.0000 0.8519 0.4167; b=conv(b1,b2);%计算H（z）的分子多项式 a=conv(a1,a2);%计算H（z）的分母多项式 w=linspace(0,pi,512); H=freqz(b,a,w); plot(w/pi,abs(H); ylabel(幅度); xlabel(Normalized frequency);

9、系统幅度响应如图1-531.9.5 部分分式法的MATLAB实现MATLAB的信号处理工具箱提供了一个计算X（z）的部分分式展开的函数，它的调用形式如下 r,p,k=resifuez(b,a)其中调用参数b,a分别表示用z表示X（z）的分子和分母多项式。如果X（z）的部分分式展开为 X(z)=则residuez的返回参数r,p,k分别为 r=r r r r p=p1 p2 p3 p3 k=k1 k2同一极点p3在向量p中出现了两次，表示p3是一个二阶的重极点。Residuez也用于由r,p,k计算z表示的X（z）的分子和分母多项式，其调用形式为 b,a=residuez(r,p,k)例1-43

10、 试用MATLAB计算X（z）=的部分分式展开。 解：计算部分分式展开的MATLAB程序如下 %program 1_6 部分分式展开 b=1.5,0.98,-2.608,1.2,-0.144; a=1,-1.4,0.6,-0.072; r,p,k=residuez(b,a); disp(留书)；disp(r) disp(极点)；disp(p) disp(常数)；disp(k);程序运行结果为留数0.7000+0.0000i 0.5000-0.0000i 0.3000极点0.6000+0.0000i 0.6000-0.0000i 0.2000常数0 2部分分式展开的结果为 X(z)=z反变换的结

11、果为 xk=(0.7x0.6+0.5x(k+1)0.6+0.3x0.2)uk+2k-11.9.6 用MATLAB计算系统的零极点MATLAB信号处理工具箱提供的tf2zp和zp2tf函数可进行系统函数的不同表示形式的转换。z正幂有理多项式表示的系统函数为 (1-154)用零点、极点和常数表示的一阶因子形式的系统函数为 （1-155）MATLAB函数z,p,k=tf2zp(b,a)将有理多项式表示的系统函数转换为一阶因子形式的系统函数。b,a=zp2tf(z,p,k)将一阶因子形式的系统函数转换为有理多项式表示的系统函数。 例1-44试将下面的系统函数表示为一阶因子形式。 （1-156） 解：用

12、z正幂有理多项式表示的系统函数为 %program 1_7 确定一阶因子形式的系统函数b=1 0 0.04 0;a=1 -0.8 0.16 -0.128;z,p,k=tf2zp(b,a);disp(零点);disp(z);disp(极点)；disp(p);disp(常数)；disp(k);程序运行结果为零点 0 0+0.2000i 0-0.2000i极点 0.8000 0.0000+0.4000i 0.0000-0.4000i常数 1系统函数的一阶因子形式为 （1-157）为了在H（z）的表达式中不出现复数，也可将式（1-157）等价的写成二阶因子的形式 （1-158）当H（z）是表示为式（1

13、-154）的形式时，可用z,p,k=tf2zp(b,a)求出系统的零极点，从而可根据极点的位置判断系统的稳定性，也可用函数zplane(b,a)画出z平面的零极点分布来判断系统的稳定性。例-45已知一因果的系统的函数为试用函数zplane画出系统的零极点分布，并判断系统的稳定性。解：%program 1_8 系统零极点分布 b=1 2 1; a=1 -0.5 -0.005 0.3; zplane(b,a);图画出了系统零极点分布。图中符号表示零点。符号旁的数字表示零点的阶数，符号x表示极点，图中的虚线画的是单位圆。由图可知，该系统的极点全在单位圆内，故系统是稳定的。.简单数字滤波器的设计例1-

14、46已知信号xk中含有角频率分别为和的离散正弦信号。试设计一高通滤波器，滤出信号xk中低频分量，保留信号xk高频分量。解：为了简单起见，假设数字滤波器是一个具有如下形式的长度为的系统h0=h2=,h1=系统的查分方程+系统的频率响应为由上式可知系统的群延迟为为了滤除低频分量，系统需满足可用命令求解上面的的方程组得：下面的MATLAB程序实现了滤波器的设计及信号的滤波。progam 1_9 简单数字滤波器的设计确定滤波器系数w1=0.015*pi;W2=0.4*pi;N=50;A=2*cos(w1) 1;2*cos(w1) 1;B=0;1x=A/B;a=1;b=x(1) x(2) x(1);%产

15、生两个余弦信号k=0:N;x1=cos(w1*k);x2=cos(w2*k);信号滤波y=filter(b,a,x1+x2);plot(k,y,r,k,x2,b-,x1+x2,k:);ylabel(幅度);xlabel(k)legend(yk,x2k,x1k+x2k);图画出了程序运行结果。由图可以看出在k=0，时，输出响应有波动。这是因为系统响应中存在瞬态响应。当时，系统进入稳态响应。滤波后的信号相对与原信号有一个样本的延迟，着是因为在时，.语音信号的滤波例1-47已知一段语音信号中混入了一频率f=1200Hz正弦型干扰信号。语音信号的抽样频率f=22050Hz。试用二阶带阻滤波器滤除语音信

16、号中的噪音信号。解：干扰信号的数字频率为由式（1-110）可得取，由式（1-111）可得解上述方程得的值为1.0619和0.9417。当的值为1.0619时，系统有一个极点在单位圆外，当的值为0.941时，系统的二个极点在单位圆内。为了保证系统的稳定，取0.9417。由和的值可得系统函数为%program 1_10 语音信号滤波Fn=1200;w=0.06;%读入语音信号x,Fs=wavread(sample);%播放语音信号sound(x,Fs);pause;%设计滤波器w0=2*pi*Fn/Fs;beta=cos(w0);alpha=min(roots(1-2/cos(w) 1);a=1,

17、-beta*(1+alpha),alpha;b=1 -2 *beta 1*(1+alpha)/2 ;%信号滤波y=filter(b,a,x);%播放滤波后语音信号sound(y,Fs);图画出了滤波前后语音信号的频谱。滤波前的语音信号的频谱在z有一高峰，这是干扰所造成的。由滤波后语音信号的频谱可看出，滤波器成功地制造了语音信号中的干扰。2.6 利用MATLAB实现信号DFT的计算2.6.1 利用MATLAB计算信号DFT在MATLAB信号处理工具箱中，函数dftmtx(N)可用来产生NN的DFT矩正D。NN的IDFT矩正D可用函数conj(dfmtx(N)/N来确定。此外，MATLAB提供了4

18、个内部函数用于计算DFT和IDFT，它们分别是： fft(x), fft(x,N), ifft(X), ifft(X,N)fft(x) 计算M点的DFT。M是序列x的长度，即M=length(x)。fft(x,N) 计算N点的DFT。若M>N，则将原序列截短为N点序列，再计算其N点DFT；若M<N，则将原序列补零至N点，然后计算其N点DFT。ifft(X) 计算M点的IDFT。M是序列X的长度。ifft(X,N) 计算N点IDFT。若M>N，则将原序列截短为N点序列，再计算其N点IDFT；若M<N，则将原序列补零至N点，然后计算其N点IDFT。为了提高fft与ifft的

19、运算效率，应尽量使序列长度M=2，或对序列补零使N=2。实现例2-7的程序如下：%program 2_1 计算有限长余弦信号频谱N=30；%数据长度L=512；%DFT的点数f1=100; f2=120; fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;f=cos(2*pi*f1*t)+cos(2*pi*f2*t);F=fftshift(fft(f,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(F);ylabel(幅度谱)实现例2-8的程序为： %program 2_2 利用Hamming计算有限长余弦信号频谱 N=50；%数

20、据长度 L=512；%DFT的点数 f1=100; f2=120; fs=600; %N=25；T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;f=cos(2*pi*f1*t)+cos(2*pi*f2*t); wh=(hamming(N); f=f.*wh; F=fftshift(fft(f,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(F);ylabel(幅度谱)例2-11 已知一长度为16点的有限序列 试用MATLAB计算序列xk的16点和512点DFT。 解： %progam 2_3 Numerical Computation of

21、DTFT Using FFT k=0:15; f=cos(2*pi*k*4./16); F_16=fft(f); F_512=fft(f,512); L=0:511; plot(L/512,abs(F_512); hold on; plot(k/16,abs(F_16),0); set(gca,xtick,0,0.25,0.5,0.75,1); set(gca,ytick,0,2,4,6,8); grid on; xlabel(Normalized frequency); ylabel(Magnirude); hold off 从图2-19可见，对长度为16的序列xk进行512点DFT，可以获

22、得频谱函数X（）更多的细节，因为序列N点的离散傅立叶变换就是序列X（）一个周期的N个等间隔抽样，尽管从信息表达的角度，由序列xk16点DFT Xm足以恢复原序列xk。2.6.2 利用MATLAB实现由DFT计算线性卷积例2-12 试利用MATLAB实现由DFT计算有限序列线性卷积，并与线性卷积直接计算的结果进行比较。 解：在MALAB中，序列线性卷积的直接结果是通过函数conv(x,h)来实现的。 %program 2_4 linear Convolution through DFT x=1 2 0 1; h=2 2 1 1; %Determine the length of the resu

23、lt of convolution L=length(x)+length(h)-1; %Compute the DFT by zero-padding XE=fft(x,L); HE=fft(h,L); %Determine the IDFT of the product y1=ifft(XE.*HE); %plot the sequence generated by circular convolution %and the error from direct linear convolution k=0:L-1; subplot(2,1,1); setm(k,real(y1);axis(0

24、 6 0 7); title(Result of Linear Convolution); xlabel(Time index k);ylabel(Amplitude); y2=conv(x,h); error=y1-y2;subplot(2,1,2); stem(k,abs(error); xlabel(Time index k);ylabel(Amplitude); title(Error Magnitude);由图2-20（b）可见，由DFT计算的线性卷积的误差非常小，这些误差主要是计算机在计算DFT时，由计算机有限字长而产生的舍入误差与计算误差。此外，MATLAB信号处理工具箱中提供了

25、fftfilt函数，该函数是基于重叠相加法的原理，可以实现长序列与短序列之间的线性卷积，其调用形式为 y=fftfilt(h,x,n)其中h:短序列；x：长序列；n：DFT的点数。一般选择n为2正整次幂，以提高DFT运算的效率。3.6 线性调频z变换算法例3-3 已知序列xk=,(),试计算序列xk在单位圆上的CZT，其中。解：由于是计算序列xk在单位圆上的CZT，故A0=1，W0=1。计算序列czt的MATLAB程序如下：%program 3_1% Compute CZT of sequence xkN=13; %length of xkn=0:N-1;xn=(0.8).n;sita=pi/

26、4; %phase of start pointfai=pi/6; %phase intervalA=exp(j*sita); %complex starting pointW=exp(-j*fai); %complex ratio between pointsM=8; %length of czty=czt(xn,M,W,A); %call czt functionsubplot(2,1,1)stem(n,xn);xlabel(k);title(sequence xk);m=0”M-1;subplot(2,1,2);stem(m,abs(y);xlabel(m);title(CZT of x

27、k);运行结果如图3-22所示。 4.5 用MATLAB实现滤波器设计4.5.1 用MATLAB实现模拟低通滤波器的设计MATLAB信号处理工具箱提供了常用的设计模拟低通滤波器的MATLAB函数。在实现应用中，可方便地调用这些函数完成模拟滤波器的设计。关于这些函数现分别介绍如下：Butterworth滤波器 N，wc=buttord(wp,ws,Ap,As,s) num,den=butter(N,wc,s)根据BW型滤波器的设计指标，利用MATLAB函数buttord获得BW型滤波器参数N和 wc。函数buttord的输入参数wp和ws(rad/s)分别表示滤波器的通带和阻带截频，Ap和As(

28、dB)表示滤波器的通带和阻带衰减。s表示所设计的是模拟滤波器。函数buttord返回参数N为BW滤波器的阶数，wc(rad/s)等于BW滤波器3 dB截频。由于wc由阻带方程确定，故由参数N、wc得出的滤波器在阻带刚好满足设计指标，在通带将超过设计指标。当BW型滤波器的参数N、wc确定后，可用MATLAB函数butter获得BW型滤波器系统函数的分子多项式（num）和分母多项式（den）。Chebyshev I型滤波器 N，wc=cheblord(wp,ws,Ap,As,s) num,den=cheby1(N,Ap,wc,s)MATLAB函数cheblord返回参数N表示CB I型滤波器的阶数

29、，wc(rad/s)等于wp。参数N、wc的取值可使cheby1设计出的CB I型滤波器在通带刚好满足设计指标。cheby1函数利用参数N、wc和Ap确定CB I型滤波器系统函数的分子多项式（num）和分母多项式（den）。Chebyshev II型滤波器 N，wc=cheb2ord(wp,ws,Ap,As,s) num,den=cheby2(N,As,wc,s)MATLAB函数cheb2ord返回参数N表示CBII型滤波器的阶数，wc(rad/s)的取值可使cheby2设计出的滤波器在通带刚好满足设计指标。cheby2函数利用参数N、wc和As确定CBII滤波器系统函数的分子多项式（num）

30、和分母多项式（den）。椭圆滤波器 N，wc=ellipord(wp,ws,Ap,As,s) num,den=ellip(N,Ap,wc,s)MATLAB函数ellipord返回参数N表示椭圆滤波器的阶数，wc=wp。MATLAB函数ellip确定阶数为N，通带衰减为Ap(dB),阻带衰减为As(dB),通带截频为wc的椭圆滤波器系统函数的分子多项式和分母多项式。在滤波器的阶数确定后，采用了4.1.3节中方案A重新计算滤波器的参数k1。故设计出的椭圆滤波器阻带衰减超过设计指标，而其他指标刚好满足设计要求。例 4-15 设计一个满足下列指标的模拟BW低频滤波器。 解：%program 4_1 设

31、计模拟Butterworth低频滤波器%滤波器指标wp=2*pi*1000;ws=2*pi*3000;Ap=1;As=50;%设计滤波器N，wc=buttord(wp,ws,Ap,As,s)num,den=butter(N,wc,s)fprintf(Order of the fiter=%.Ofn,N)disp(Numerator polynomial);fprintf(%.4en,num);disp(Denominator polynomial);fprintf(%.4en,den);%计算Ap,As及增益响应omega1=linspace(0,wp,500);omega2=linspace

32、(wp,ws,200);omega3=linspace(ws,5*1000*pi*2,500);h1=20*log10(abs(freqs(num,den,omega1);h2=20*log10(abs(freqs(num,den,omega2);h3=20*log10(abs(freqs(num,den,omega3);fprintf(Ap=%.4fn,max(-h1);fprintf(As=%.4fn,min(-h3);plot(omega1 omega2 omega3/(2*pi),h1 h2 h3);grid;xlabel(Frequency in Hz);ylabel(Gain in

33、 dB);程序运行后的结果如下：order of the filter=6 Numerator polynomial0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+0000.0000e+000 0.0000e+000 1.4184e+023Denominator polynomial1.000e+000 2.7902e+004 3.8927e+008 3.4429e+0122.0301e+016 7.5889e+019 1.4184e+023Ap=0.7488 As=50.0000图4-25画出了该滤波器的增益响应。例4-16 利用模拟椭圆低通滤波器重新

34、设计例4-15。解：只需将例4-15程序的第6行和第7行改为N，wc=ellipord(wp,ws,Ap,As,s) num,den=ellip(N,Ap,wc,s)即可。程序的运行结果如下： Order of the filter=4 Numerator polynomial3.1615e-003 1.8190e-012 3.2366e+006 1.5259e-004 4.4736e+014Denominator polynomial1.0000e+000 5.9359e+003 5.8689e+007 1.9320e+011 5.0195e+014Ap=1.0000 As=51.3315图

35、4-26画出了该滤波器的增益响应。比较例4-15和例4-16可知，实现同样设计指标的滤波器，椭圆滤波器所需阶数较低。 4.5.2 用MATLAB实现模拟域频率变换MATLAB信号处理工具箱提供了实现四种模拟域频率变换的函数，它们分为：低通到低通的变换 numt,dent=1p21p(num,den,w0)低通到高通的变换 numt,dent=1p2hp(num,den,w0)低通到带通的变换 numt,dent=1p2bp(num,den,w0,B)低通到带阻的变换 numt,dent=1p2bs(num,den,w0,B)其中num、den分别表示变换前模拟滤波器系统函数的分子多项式和分母多

36、项式。numt、dent分别表示变换后模拟滤波器系统函数的分子多项式和分母多项式。w0和B为变换中的参数。 例4-47 试设计一个满足下列指标的模拟BW型带阻滤波器 解： （1）由式（4-77）确定低通到带阻变换中的参数 B=2， （2）由式（4-79） 由式（4-80）得原型低通滤波器的通带截频wp_bar为 wp_bar=0.3714 （3）由 N，wc=buttord(wp_bar,1,Ap,As,s); num,den=butter(N,wc,s);得原型低通滤波器为 （4）由 numt,dent=1p2bp(num,den,w0,B)；得带阻低通滤波器为 图4-27画出了该滤波器的增

37、益响应曲线。该滤波器的通带和阻带衰假分别为Ap1=0.05dB, Ap2=0.68dB, As1=As2=10dB。 4.5.3 脉冲响应不变的MATLAB实现 MATLAB信号处理提供的impinvar(num,den,Fs)函数，可实现用脉冲响应不变法将模拟滤波器转化为数字滤波器，起调用形式为： numd,dend=impinvar(num,den,Fs)式中num和den分别表示模拟滤波器系数函数H（s）的分子多项式和分母多项式，Fs是脉冲响应不变法中的抽样频率，单位是Hz。输出变量numd和dend分别表示数字滤波器的系统函数H(z)分子多项式和分母多项式。若某个因果模拟滤波器的系统函

38、数为 取T=1，则由 numd,dend=impinvar(2 14,1 2 5,1)得 numd=2.0000 2.3133 dend=1.0000 0.3062 0.1353所以由脉冲响应不变法获得的数字滤波器为 例4-18 利用Butterworth低通滤波器及脉冲响应不变法设计满足下列指标的数字滤波器 解：设计满足上述指标数字滤波器的MATLAB程序如下： %program 4_2 脉冲响应不变法设计BW型LP DF %DF BW LP 指标Wp=0.1*pi;Ws=0.4*pi;Ap=1;As=25; Fs=1;%抽样频率（Hz）%确定模拟BW 指标Wp=Wp*Fs;Ws=Ws*Fs

39、;%确定AF阶数N=buttord(Wp,Ws,Ap,As,s);%由通带指标确定3-dB截频Wc=Wp/(10(0.1*Ap)-1)(1/2/N);%确定BW AFnuma,dena=butter(N,Wc,s);%确定 DFnumd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);%幅度归一化DF的幅度响应norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm);grid;xlabel(Normalized frequency);y

40、label(Gain,dB);disp(Numerator polynomial);fprintf(%.4en,numd);disp(Denominator polynomial);fprintf(%.4en,dend);%计算Ap和Asw=Wp Ws;h=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log10(abs(h(1);fprintf(As=%.4fn,-20*log10(abs(h(2);运行结果为Numerator polynomial-5.5515e-0172.3231e-0021.7880e-0020.0000e+000Denominato

41、r polynomial1.0000e+000-2.2230e+0001.7193e+000-4.5520e-001Ap=0.9985As=30.3240Ap=0.9985As=30.3240图4-28画出了该滤波器的增益响应 4.5.4 双线性变换法的MATLAB实现 MATLAB信号处理工具箱提供的（num,den,Fs）函数可以用来实现双线性变换，其调用形式为： numd,dend=bibinear(num,den,Fs)式中num和den分别表示模拟滤波器系统函数H（s）的分子多项式和分母多项式，Fs=1/T.输出变量numd和dend分别数字滤波器系统函数H（s）的分子多项式和分母多

42、项式。例4-19 三阶归一化Butterworth滤波器的系统函数为 解：取T=2，则由 den=conv(1 1,1 1 1); numd,dend=bilinear(1,den,0.5)得经双线性变换后的数字滤波器的分子多项式和分母多项式为 numd=0.1667 0.5000 0.5000 0.1667 dend=1.0000 0.0000 0.3333 0.0000注意0.1667是1/6的近似值，0.3333是1/3的近似值，所以双线性变换后的数字滤波器的系统函数H（z）为 例4-20 利用CB I型模拟滤波器和双线性变换法，设计一满足下列指标的数字低通滤波器。 解：取T=2，由得模

43、拟滤波器的频率指标 由 N,wc=cheb1ord(0.1854,0.7265,1,10,s) num,den=cheby1(N,1,wc,s)得模拟滤波器的分子多项式和分母多项式为 num= 0 0 0.0246 den=1.0000 0.1739 0.0277即 再由 numd,dend=bilinear(num,den,0.5)得双线性变换后的数字滤波器的分子多项式和分母多项式 numd=0.0205 0.0410 0.0205 dend=1.0000 -1.6185 0.7106即 上式表示的数字滤波器的增益响应如图4-29所示。该滤波器的通带和阻带衰减分别为。 4.5.5 用MATL

44、AB实现数字滤波器设计MATLAB也提供了直接设计IIR数字滤波器的函数。设计的基本思想是用式（4-103）将数字滤波器的频率指标转化为模拟滤波器的指标，然后设计模拟滤波器，最后用双线性变换把模拟滤波器转换为数字滤波器。由于在设计中用的模拟滤波器的类型不同，所以给出了以下四组不同的函数BW型数字低通滤波器 N,Wc=buttord(Wp,Ws,Ap,As) num,den=butter(N,Wc)在函数buttord函数中，调用参数Wp、Ws是数字低通滤波器的归一化的通带和阻带截频。例如要求数字滤波器的，则Wp=0.1,Ws=0.4。Ap和As为滤波器的通带和阻带的衰减（dB）。返回的参数为和

45、Wc，其中N为滤波器阶数。函数butter中，调用参数N和Wc由函数buttord确定。返回参数num和den分别是数字滤波器的分子多项式和分母多项式的系数。CB I型数字低通滤波器 N,Wc=cheb1ord(Wp,Ws,Ap,As) num,den=cheby1(N,Ap,Wc)CB II 型数字低通滤波器 N,Wc=cheb1ord(Wp,Ws,Ap,As) num,den=cheby1(N,Ap,Wc)椭圆型数字低通滤波器 N,Wc=ellipord(Wp,Ws,Ap,As) num,den=ellip(N,Ap,As,Wc)利用上述四组MATLAB函数，也可设计数字高通、带通和带阻滤

46、波器。调用方法可参见MATLAB手册。 例4-21 设计满足下列指标的BW型和椭圆型数字低通滤波器。 解： %program 4_3 设计数字低通滤波器 %DF 指标 Wp=0.2;Ws=0.4; Ap=0.5;As=20; %设计BW型DF LP N,Wc=buttord(Wp,Ws,Ap,As); b1,a1=butter(N,Wc); %设计C型DF LP N,Wc=ellipord(Wp,Ws,Ap,As); b2,a2=ellip(N,Ap,As,Wc); disp(BW型分子多项式)； fprintf(%.4en,b1); disp(BW型分母多项式)； fprintf(%.4en

47、,a1); disp(C型分子多项式)； fprintf(%.4en,b2); disp(C型分母多项式)； fprintf(%.4en,a2); w=linspace(0,0.6*pi,200); h1=freqz(b1,a1,w); h2=freqz(b2,a2,w); plot(w/pi,20*log10(abs(h1),w/pi,20*log10(abs(h2); legend(BW型，C型)； xlabel(Normalized frequency); ylabel(Gain,dB); axis(0 0.6 -50 1);程序运行的结果为 BW型分子多项式 4.7792e-003 2

48、.3896e-002 4.7792e-002 4.7792e-002 2.3896e-002 4.7792e-003 BW型分母多项式 1.0000e+000 -2.2360e+000 2.4061e+000 -1.3906e+000 4.2865e-001 -5.5151e-002 C型分子多项式 9.4379e-002 -1.5627e-002 -1.5627e-002 9.4379e-002 C型分母多项式 1.0000e+000 -1.9853e+000 1.6032e+000 -4.6040e-001 图4-30画出了所设计的滤波器的增益响应。由运行结果知对同样的设计指标，用BW型模拟滤波器设计出的数字滤波器是五阶的，而用C型模拟滤波器设计出的数字滤波器只有三阶。 5.5 利用MATLAB实现FIR滤波器设计5.5.1 窗函数法的MATLAB实现窗函数的计算 MATLAB提供了许多常用的窗函数，其中部分窗函数的调用形式为 w=hanning(N) w=hamming(N) w=blackman(N) w=Kaiser(N,beta)其中N是窗函数的长度,