数字信号处理实验4-6

1、实验4 离散系统的变换域分析一、实验目的1、熟悉对离散系统的频率响应分析方法；2、加深对零、极点分布的概念理解。 二、实验原理离散系统的时域方程为其变换域分析方法如下：频域:系统的频率响应为:Z域:系统的转移函数为:分解因式:，其中和称为零、极点。三、预习要求1. 在MATLAB中，熟悉函数tf2zp、zplane、freqz、residuez、zp2sos的使用，其中：z，p，K=tf2zp（num，den）求得有理分式形式的系统转移函数的零、极点；zplane（z，p）绘制零、极点分布图；h=freqz(num,den,w)求系统的单位频率响应；r，p，k=residuez（num，den

2、）完成部分分式展开计算；sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。 2. 阅读扩展练习中的实例，学习频率分析法在MATLAB中的实现； 3. 编程实现系统参数输入，绘出幅度频率响应和相位响应曲线和零、极点分布图。 四、实验内容求系统的零、极点和幅度频率响应和相位响应。解析：【代码】num=0.0528 0.0797 0.1295 0.1295 0.797 0.0528;den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;z,p,k=tf2zp(num,den);disp('零点');disp(z);disp(

3、9;极点');disp(p);disp('增益系数');disp(k);figure(1)zplane(num,den)figure(2)freqz(num,den,128) 【图形】【结果】零点 -1.5870 + 1.4470i -1.5870 - 1.4470i 0.8657 + 1.5779i 0.8657 - 1.5779i -0.0669 + 0.0000i极点 0.2788 + 0.8973i 0.2788 - 0.8973i 0.3811 + 0.6274i 0.3811 - 0.6274i 0.4910 + 0.0000i增益系数 0.0528五、扩展

4、练习例1: 求下列直接型系统函数的零、极点，并将它转换成二阶节形式解:用MATLAB计算程序如下：num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(num,den);m=abs(p);disp('零点');disp(z);disp('极点');disp(p);disp('增益系数');disp(k);sos=zp2sos(z,p,k);disp('二阶节');disp(real(sos);zplane(num,den)输入到“num”和“den”的分别为分子和分

5、母多项式的系数。计算求得零、极点增益系数和二阶节的系数：零点:0.9615-0.5730-0.1443 + 0.5850i-0.1443 - 0.5850i极点:0.5276+0.6997i0.5276-0.6997i-0.5776+0.5635i-0.5776-0.5635i增益系数:1二阶节:1.0000 -0.3885 -0.5509 1.0000 1.15520 0.65111.0000 0.28850 0.36300 1.0000 -1.0552 0.7679系统函数的二阶节形式为：极点图见图:图 系统函数的零、极点图例2: 差分方程所对应的系统的频率响应。解:差分方程所对应的系统函

6、数为：用MATLAB计算的程序如下：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);gridtitle('实部')xlabel('omega/pi');ylabel('幅度')subplot(2,2,2);plot(w/pi,imag(h);gridtitle('虚部')xlabel('omega/pi');ylabel('Am

7、plitude')subplot(2,2,3);plot(w/pi,abs(h);gridtitle('幅度谱')xlabel('omega/pi');ylabel('幅值')subplot(2,2,4);plot(w/pi,angle(h);gridtitle('相位谱')xlabel('omega/pi');ylabel('弧度') 图实验5 有限冲激响应数字滤波器设计一、实验目的： 1、加深对数字滤波器的常用指标理解。 2、学习数字滤波器的设计方法。 二、实验原理：低通滤波器的常用指标

8、：（1）通带边缘频率；（2）阻带边缘频率；（3）通带起伏；（4）通带峰值起伏，（5）阻带起伏，最小阻带衰减。三、预习要求1. 在MATLAB中，熟悉函数fir1、kaiserord 、remezord、remez的使用； B = fir1(n,Wn,'high','noscale')设计滤波器； n,Wn,beta,ftype = kaiserord(f,a,dev)估计滤波器阶数； n,fo,ao,w = remezord (f,a,dev,fs)计算等波纹滤波器阶数n和加权函数w()； B=remez(n,f,a)进行等波纹滤波器的设计。 2. 阅读扩展练习

9、中的实例，学习FIR滤波器的设计方法及其在MATLAB中的实现； 3. 给出FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。 数字滤波器有IIR和FIR两种类型，它们的特点和设计方法不同。四、实验内容：利用MATLAB编程,分别用窗函数法和等波纹滤波器法设计两种FIR数字滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：。阻带边缘频率：，最小阻带衰减：。 l 用窗函数法设计： 【解析】调用函数n,wn,bta,ftype=kaiserord（f，a，dev，fs）参数：f=0.3 0.45 0.65 0.8为对应

10、数字频率a=0 1 0为由f指定的各个频带上的幅值向量，一般只有0和1表示；和f长度关系为（2*a的长度）2=（f的长度） devs=0.01 0.1087 0.01用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等，计算公式：阻带衰减误差=，通带波动衰减误差=fs缺省值为2HZ【代码】n,wn,bta,ftype=kaiserord(0.3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用kaiserord函数估计出滤波器阶数n和beta参数h1=fir1(n,wn,ftyp

11、e,kaiser(n+1,bta),'noscale');hh1,w1=freqz(h1,1,256);figure(1)subplot(2,1,1)plot(w1/pi,20*log10(abs(hh1)gridxlabel('归一化频率w');ylabel('幅度/db')subplot(2,1,2)plot(w1/pi,angle(hh1)gridxlabel('归一化频率w');ylabel('相位/rad');【图形】l 用等波纹法设计： 【解析】调用函数n,fpts,mag,wt=remezo

12、rd(f，a，dev) f=0.3 0.45 0.65 0.8 a=0 1 0 dev=0.01 0.1087 0.01 其含义同函数n,wn,bta,ftype=kaiserord（f，a，dev，fs）中的参数相同。【代码】n,fpts,mag,wt=remezord(0.3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w，使rem

13、ez函数设计出的滤波器满足f、a及dev指定的性能要求h2=remez(n,fpts,mag,wt);%设计出等波纹滤波器hh2,w2=freqz(h2,1,256);figure(2)subplot(2,1,1)plot(w2/pi,20*log10(abs(hh2)gridxlabel('归一化频率w');ylabel('幅度/db')subplot(2,1,2)plot(w2/pi,angle(hh2)gridxlabel('归一化频率w');ylabel('相位/rad');【图形】五、扩展练习例1: 用凯塞窗设计一FIR

14、低通滤波器，通带边界频率，阻带边界频率，阻带衰减不小于50dB。解: 首先由过渡带宽和阻带衰减来决定凯塞窗的N和上图给出了以上设计的频率特性，(a) 为N=30直接截取的频率特性(b)为凯塞窗设计的频率特性。凯塞窗设计对应的MATLAB程序为：wn=kaiser(30,4.55);nn=0:1:29;alfa=(30-1)/2;hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa);h=hd.*wn'h1,w1=freqz(h,1);或者：b = fir1(29,0.4,kaiser(30,4.55);h1,w1=freqz(b,1);plot(w1/pi,20*

15、log10(abs(h1);axis(0,1,-80,10); grid;xlabel('归一化频率/p');ylabel('幅度/dB');还可以使用n,Wn,beta,ftype=kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器：fcuts = 0.3 0.5; %归一化频率omega/pimags = 1 0;devs = 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs); %计算出凯塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n

16、+1,beta),'noscale');freqz(hh);实际中，一般调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W()，调用函数remez可进行等波纹滤波器的设计，直接求出滤波器系数。函数remezord中的数组fedge为通带和阻带边界例2：利用雷米兹交替算法设计等波纹滤波器，设计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动阻带最小衰减At=40dB，采样频率fs=4000Hz。解:在MATLAB中可以用remezord 和remez两个函数设计，其结果如图2，MAT

17、LAB程序如下：fedge=800 1000; mval=1 0; dev=0.0559 0.01; fs=4000; N,fpts,mag,wt=remezord(fedge,mval,dev,fs);b=remez(N,fpts,mag,wt);h,w=freqz(b,1,256);plot(w*2000/pi,20*log10(abs(h);grid; xlabel('频率/Hz');ylabel('幅度/dB');所得图像如下所示：实验6 无限冲激响应数字滤波器设计一、实验目的1、掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方

18、法;2、熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。二、实验原理在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：1)利用buttord和cheb1ord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；2)num,den=butter（N,Wn）(巴特沃斯)和num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以进行滤波器的设计；3）lp2hp，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；4)使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波

19、器的传输函数系数；5)利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。三、预习要求1. 在MATLAB中，熟悉函数butter、cheby1、cheby2的使用，其中：num,den=butter（N,Wn）巴特沃斯滤波器设计；num,den=cheby1（N,Wn）切比雪夫1型滤波器设计；num,den=cheby2（N,Wn）切比雪夫2型滤波器设计。2. 阅读扩展练习中的实例，学习在MATLAB中进行数字滤波器的设计； 3. 给出IIR数字滤波器参数和滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。 四、实验内容利用MATLAB编程，

20、用脉冲响应不变法和双线性变换法设计一个数字带通滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：；阻带边缘频率：，最小阻带衰减： 。【解析】l 巴特沃兹原型1 双线性变换法【代码】ws1=2*8000*tan(0.3*pi/2);ws2=2*8000*tan(0.8*pi/2);wp1=2*8000*tan(0.45*pi/2);wp2=2*8000*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2; Rp=1;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,

21、A=bilinear(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1);plot(f,20*log10(abs(h);axis(0,4000,-60,10);gridxlabel('频率/Hz');ylabel('幅度/dB');subplot(2,1,2)plot(f,angle(h);gridxlabel('频率/Hz');ylabel('相位');【图形】2 脉冲响应不变法【代码】fs=8000;ws1=0.3*pi*fs;ws2=0.8*pi*fs;wp1=0.

22、45*pi*fs;wp2=0.65*pi*fs;ws=ws1 ws2;wp=wp1 wp2;Rp=1;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,A=impinvar(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1)plot(f,20*log10(abs(h);axis(0,4000,-80,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');subplot(

23、2,1,2)plot(f,angle(h)grid;xlabel('频率/Hz');ylabel('相位');【图形】fs=2;ws1=2*fs*tan(0.3*pi/2);ws2=2*fs*tan(0.8*pi/2);wp1=2*fs*tan(0.45*pi/2);wp2=2*fs*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2;Rp=1;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,A=bilinear(num,den,f

24、s);h1,w=freqz(B,A);ws11=0.3*pi*fs;ws22=0.8*pi*fs;wp11=0.45*pi*fs;wp22=0.65*pi*fs;ws0=ws11 ws22;wp0=wp11 wp22; N,Wn=buttord(wp0,ws0,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,A=impinvar(num,den,fs);h2,w=freqz(B,A);figure(1)plot(w/pi,20*log10(abs(h1),'-.',w/pi,20*log10(abs(h2),'

25、-');axis(0,1,-80,0);grid;xlabel('频率/Hz');ylabel('幅度/dB');l 切必雪夫原型1 双线性变换法【代码】ws1=2*8000*tan(0.3*pi/2);ws2=2*8000*tan(0.8*pi/2);wp1=2*8000*tan(0.45*pi/2);wp2=2*8000*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2;Rp=1;Rs=40;N,Wn=cheb1ord(wp,ws,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s&

26、#39;);B,A=bilinear(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1)plot(f,20*log10(abs(h);axis(0,4000,-60,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');subplot(2,1,2)plot(f,angle(h);grid;xlabel('频率/Hz');ylabel('相位');【图形】2 脉冲响应不变法【代码】fs=8000;ws1=0.3*pi*fs;ws2=0.8*pi*

27、fs;wp1=0.45*pi*fs;wp2=0.65*pi*fs;ws=ws1 ws2;wp=wp1 wp2; Rp=1;Rs=40;N,Wn=cheb1ord(wp,ws,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s');B,A=impinvar(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1);plot(f,20*log10(abs(h);axis(0,4000,-90,10);grid;xlabel('频率/Hz');ylabel('幅度/dB&

28、#39;);subplot(2,1,2);plot(f,angle(h);grid;xlabel('频率/Hz');ylabel('相位');【图形】fs=8000;ws1=2*fs*tan(0.3*pi/2);ws2=2*fs*tan(0.8*pi/2);wp1=2*fs*tan(0.45*pi/2);wp2=2*fs*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2; Rp=1;Rs=40;N,Wn=cheb1ord(wp,ws,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s')

29、;B,A=bilinear(num,den,fs);h1,w=freqz(B,A);ws11=0.3*pi*fs;ws22=0.8*pi*fs;wp11=0.45*pi*fs;wp22=0.65*pi*fs;ws0=ws11 ws22;wp0=wp11 wp22;N,Wn=cheb1ord(wp0,ws0,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s');B,A=impinvar(num,den,fs);h2,w=freqz(B,A);figure(1)plot(w/pi,20*log10(abs(h1),'-.',w/

30、pi,20*log10(abs(h2),'-');axis(0,1,-80,10);grid;xlabel('频率/Hz')ylabel('幅度/dB');【结论】 双线性变换法通过将数字频率的取值范围从0到p对应到模拟频率的范围0到的范围，也就对应于模拟域中所有可能的频率值。双线性变换法不会出现频率混叠，但非线性关系却导致数字滤波器的频率响应不能逼真地模仿模拟滤波器的频率响应。脉冲响应不变法通过选择满足设计要求的模拟滤波器冲激响应h（t）的采样值的数字脉冲响应hn得到的被采样的冲激响应将给出与原模拟滤波器非常相近的滤波器形状。由于该方法不可避免

31、的要发生频率混叠现象，所以只适合设计低通和带通滤波器。从实验结果可以看出：双线性变换法所设计的巴特沃兹滤波器最符合设计指标，而用脉冲响应不变法设计的滤波器（无论是巴特沃兹还是切比雪夫）有一定的误差，主要是由于混叠所引起。五、扩展练习例1：设采样周期T=250s（采样频率fs =4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc =1kHz。B,A=butter(3,2*pi*1000,'s');num1,den1=impinvar(B,A,4000);h1,w=freqz(num1,den1);B,A=butter(3,2/0.0002

32、5,'s');num2,den2=bilinear(B,A,4000);h2,w=freqz(num2,den2);f=w/pi*2000; plot(f,abs(h1),'-.',f,abs(h2),'-');grid;xlabel('频率/Hz ')ylabel('幅值/dB')程序中第一个butter的边界频率2×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。同时，也看到双线性变换法，在z=-1即=或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在=处的三阶传输零点通过映射形成的。下图给出了MATLAB计算的结果。 例2: