实验四IIR数字滤波器设计及软件实现实验报告.docx

实验四IIR数字滤波器设计及软件实现实验报告一、实验目的熟悉用双线性变换法设计1IR数字滤波器的原理与方法；(1) 学会调用MATLAB信号处理工具箱中滤波器设计函数(或滤波器设计分析工具 fdatool)设计各种IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。(2) 掌握IIR数字滤波器的MATLAB实现方法。(3)通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。二、实验原理设计11R数字滤波器一般采用间接法(脉冲响应不变法和双线性变换法),应用最广泛 的是双线性变换法。基本设计过程是：先将给定的数字滤波器的指标转换成过渡模拟滤波器 的指标；设计过渡模拟滤波器；将过渡模拟滤波器系统函数转换成数字滤波器的系统函 数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六 章介绍的滤波器设计函数butter、chebyl、cheby2和ellip可以分别被调用来直接设计 巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函 数直接设计IIR数字滤波器。本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给 定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n)。三、实验内容及步骤(1)调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st,该函 数还会自动绘图显示st的时域波形和幅频特性曲线，如图所示。由图可见，三路 信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离， 这就是本实验的目的。(b) s(t)的频谱I4IIIII| - - - -4-iiiii翅；：；：los 0.514- !1-e-IIIII200400600800 1000 1200 U00 1600 1800 2000f/Hz图三路调幅信号st的时域波形和幅频特性曲线要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分 第1页共12页3. 具体步骤：(1) 构造原始信号s(t)(2) 画出s(t)的频谱(3) 设计ellipse数字滤波器(IIR),包括低通，带通，带通，并显示幅频特 性(4) 用得到的滤波器进行滤波，分离出三路信号，观察时域波形和幅频特性(5) 用三路信号si, s2, s3尝试重新合成原始信号问题：为什么重新合成的信号和原信号不相等呢？谁能解释一下？谢谢 程序如下：clearelf% (1)构造原始信号Fs二10000; T二1/Fs; %先设定采样频率t=0:T:0. 1;n=length(t);s=cos (2*pi*250*t). *cos(2*pi*25*t)+cos(2*pi*500*t). *cos(2*pi*50*t) +. cos (2*pi*1000*t) *cos(2*pi*100*t);subplot (2, 1, 1), plot (t, s), axis (0 0. 08 2 3)titleC原始信号s(t)xlabel ( t/s), ylabel ( s (t)% (2)画出s(t)的频谱ft=fftshift(fft(s);i二fix(n/2) ;f=(-i:i)/n*Fs;%貌 似 这 是 公 式。subplot(2, 1, 2), stem(f, abs(ft), Marker, none), xlim(0 1250) title C s (t)的频谱)xlabel (J f/Hz，), ylabel C 幅度)% ( 3 )设计 ellipse数字滤波器(IIR),并显示幅度特性%3. la设计模拟低通滤 波 器 fp=320:fs=400:Ap=0. 1:As=60:wp二2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp= (2*Fs)*tan(wp/2);Ws=(2*Fs)*tan(ws/2);N, Wc= ellipord(Wp , Ws , Ap , As ,s);bLPs , aLPs=ellip(N , Ap, As, Wc, s);H, w=freqs (bLPs, aLPs);db =20*logl0(abs(H):figure, subplot (2, 1, 1), plot (w/2/pi, db);axis (0 1600 -80 5),gridtitle (模拟低通滤波器的幅度特性);xxxxxx C f (Hz);yyyyyy (，dB)；%3. lb将模拟低通滤波器转换为数字低通滤波器bLPz, aLPz = bilinear (bLPs , aLPs , Fs);w 二 linspace (0, pi, 1000);h = freqz(bLPz, aLPz , w);subplot (2,1, 2), plot (w*Fs/2/pi, 20*logl0(abs(h);axis (0 1600 -80 5),grid titlef数字低通滤波器的幅度特性);xxxxxx (f (Hz); yyyyyyCdB)；%3. 2a设计模拟带通滤波器 fp=430 570;fs=330 670;Ap=0. l;As=60;wp二2*pi*fp/Fs;ws二2*pi*fs/Fs;Wp= (2*Fs)*tan(wp/2);Ws=(2*Fs)*tan(ws/2);N, Wc= ellipord(Wp , Ws , Ap , As , s);bBPs , aBPs=ellip (N , Ap, As, Wc, s);H, w =freqs (bBPs, aBPs);db =20*logl0(abs(H): figure, subplot (2, 1, 1), plot (w/2/pi, db);axis (0 1600 -80 5),grid title (模拟带通滤波器的幅度特性)；xxxxxx ( f (Hz);yyyyyy ( dB)%3. 2b将模拟带通滤波器转换为数字带通滤波器bBPz,aBPz = bilinear (bBPs , aBPs , Fs);w = linspace (0, pi, 1000);h = freqz (bBPz, aBPz , w); subplot (2, 1, 2), plot(w*Fs/2/pi, 20*logl0(abs(h);axis(0 1600 -80 5),grid titleC数字低通滤波器的幅度特性)； xxxxxx (f (Hz);yyyyyy (dB)； %3. 3a设计模拟高通滤波器fp=800;fs=700;Ap=0. 1;As=60; wp二2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp=(2*Fs)*tan(wp/2);Ws=(2*Fs)*tan(ws/2);N, Wc= ellipord (Wp , Ws , Ap , As ,s);bHPs , aHPs=ellip (N, Ap, As, Wc , high, s);H, w=freqs (bHPs, aHPs);db 二20*logl0(abs(H); figure, subplot (2, 1, 1), plot (w/2/pi, db);axis(0 1600 -80 5),grid titleC模拟高通滤波器的幅度特性)； xxxxxx ( f (Hz);yyyyyy (，dB)%3. 3b将模拟高通滤波器转换为数字高通滤波器bHPz, aHPz = bilinear(bHPs , aHPs , Fs);w = linspace(0, pi, 1000);h = freqz (bHPz, aHPz , w);subplot (2, 1, 2), plot(w*Fs/2/pi, 20*logl0(abs(h); axis(0 1600 -80 5),grid title (数字高通滤波器的幅度特性)；xxxxxx (f (Hz); yyyyyy(dB)；% (4)用得到的滤波器进行滤波，分离出三路信号，观察时域波形和幅频特 性%4. 1滤波si = filter(bLPz, aLPz, s); s2 二 filter(bBPz, aBPz, s); s3 = filter (bHPz, aHPz, s);%4. 2时域波形figure, subplot (3, 1, 1), plot (t, si) xlabel ( t/s ), ylabel ( si (t) subplot (3, 1, 2), plot (t, s2) xlabel ( t/s), ylabel ( s2 (t) subplot (3, 1, 3), plot (t, s3) xlabel C t/s，), ylabel ( s3 (t) %4. 3幅频特性ftl=fftshift(fft(si):ft2=fftshift(fft(s2) ;ft3=fftshift(fft(s3); figure, subplot (3, 1, 1), stem(f, abs (ftl), Marker5, none), xlim(0 1250) subplot(3, 1, 2), stem(f, abs(ft2), Marker5,none),xlim(0 1250) subplot(3, 1, 3), stem(f, abs(ft3), Marker5, none，), xlim(0 1250) %(5)用三路信号sl,s2, s3尝试重新合成原始信号ss=sl+s2+s3;%合成信号记作ss isequal (s, ss)%判断合成信号是否和原信号相等，结果不相等，why?本帖最后由cwjy于2010-5-318:19编辑离st中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器) 的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0. IdB,阻带最小衰减为60dB。提示：抑制载波单频调幅信号的数学表示式为1s(A) cos(2 ft) cos(2 ft) cos(2 (f f cos(2 (f f)/)0c 2c oc 0其中，cos(2)称为载波，匚为载波频率，cos(2 Q)称为单频调制信号，f为调制正 弦波信号频率，且满足。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频r f和差频f /,这2个频率成分关于载波频c 0C 0率 对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率f对称的2根谱线， f CC其中没有载频成分，故取名为抑制载波单频调幅信号。容易看出，图10. 4.1中三路调幅信 号的载波频率分别为250Hz、500Hz、1000Hzo如果调制信号m(t)具有带限连续频谱，无直流成分，则s(，) m(r)cos(2 ft)就是一般的抑制载波调幅信号。其频谱图是关于载波频率匚对称的2个边带(上下边带)，砖业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅信号，简称双边带(DSB)信号。如果调制信号rn(t)有直流成分，则s(f) m(0cos(2 ffl就是一般的双边带调幅信号。其频谱图是关于载波频率匚对称的2 个边带(上下边带),并包含载频成分。(3) 编程序调用MATLAB滤波器设计函数ellipord和ellip分别设计这三个椭圆滤波 器，并绘图显示其幅频响应特性曲线。(4) 调用滤波器实现函数filter,用三个滤波器分别对信号产生函数mstg产生的信 号st进行滤波，分离出st中的三路不同载波频率的调幅信号yjn)、y2(n)和y3(n),并 绘图显示yl(n)、y2(n)和y3(n)的时域波形，观察分离效果。四、信号产生函数mstg清单function st=mstg%产生信号序列向量st,并显示st的时域波形和频谱 %st=mstg返回三路调幅信号相加形成的混合信号，长度N=1600N=1600%N为信号st的长度。Fs=10000;T=l/Fs;Tp=N*T; %采样频率 Fs=10kHz, Tp 为采样时间 t=O:T: (N-l)T;k=O:N-l;f=k/Tp;fcl=Fs/10;%第1路调幅信号的载波频率fcl=1000Hz,fml=fcl/10;%第1路调幅信号的调制信号频率fml=100Hzfc2二Fs/20;%第2路调幅信号的载波频率fc2=500IIzfm2=fc2/10;%第2路调幅信号的调制信号频率fm2=50Hzfc3=Fs/40;%第3路调幅信号的载波频率fc3=250Hz,fm3=fc3/10;%第3路调幅信号的调制信号频率fm3=25Hzxtl=cos(2*pi*fml*t) *cos(2*pi*fcl*t) ; %产生第 1 路调幅信号 xt2二cos (2*pi*fm2*t). *cos (2*pi*fc2*t) ; %产生第 2 路调幅信号 xt3=cos(2*pi*fm3*t). *cos(2*pi*fc3*t); %产生第 3 路调幅信号st=xtl+xt2+xt3; fxt=fft (st,N);st=xtl+xt2+xt3; fxt=fft (st,N);%三路调幅信号相加 %计算信号st的频谱%=以下为绘图部分，绘制st的时域波形和幅频特性曲线二 subplot (3, 1, 1)plot (t, st) ;grid;xlabel ( t/s，) ;ylabel ( s (t);axis(0, Tp/8, min (st), max (st) ; title (a) s(t)的波形) subplot (3, 1, 2)stem(f, abs (fxt)/max(abs (fxt),);grid; title( (b) s(t)的频谱) axis (0, Fs/5, 0, 1.2); xlabel ( f/Hz) ;ylabel (幅度)五、实验程序框图如图io. 4. 2所示调用函数mstg产生st,自动绘图 显示st的时域波形和幅频特性曲线调用ellipord和ellip分别设计三个椭圆滤波器，并绘图显示其幅频响应特性曲线。出三路不同载波频率的调幅信号y】（n）、y 2（n）和fn）六、滤波器参数及实验程序清单1、滤波器参数选取观察图可知，三路调幅信号的载波 频率分别为250Hz、500Hz. lOOOHzo带宽（也可以由信号产生 函数mstg清单看出）分别为50Hz、100Hz、200Hzo所以，分离混合 信号st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波 器）的指标参数选取如下：对载波频率为250Hz的条幅信号，可以用低通滤波器分离，其指标为带截止频率fp280Hz通带最大衰减p 0.1dBdB;阻带截止频率f 450Hz阻带最小衰减s 60dBdB,对载波频率为500Hz的条幅信号，可以用带通滤波器分离，其指标为带截止频率。/440也，fpu 56C)hz,通带最大衰减p O.ldBdB;阻带截止频率275Hz)f 900Hz hz,阻带最小衰减sis 60dBdB,对载波频率为1000Hz的条幅信号，可以用高通滤波器分离，其指标为带截止频率。89Hz,通带最大衰减p dBdB；阻带截止频率55C)hz,阻带最小衰减、60dBdB,S说明：(1)为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽量使滤波 器过渡带宽尽可能宽。(2) 与信号产生函数mstg相同，采样频率Fs=10kHzo为了滤波器阶数最低，选用椭圆滤波器。按照图所示的程序框图编写的实验程序为exp4.ni。2、实验程序清单%实验4程序% HR数字滤波器设计及软件实现clear all;close allFs = 10000; T=l/Fs;%采样频率%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号stst=mstg;%低通滤波器设计与实现fp=280;fs=450;wp二2*fp/Fs;ws=2*fs/Fs;rp二0. 1 ;rs=60;%DF指标(低通滤波器的通、阻带边界频)N, wp=ellipord(wp, ws, rp, rs) ; %调用ellipord计算椭圆DF阶数N和通带截止频率wpB, Aellip(N, rp, rs,wp);%调用ellip计算椭圆带通DF系统函数系数向量B和Aylt二filter (B, A, st) ;%滤波器软件实现%低通滤波器设计与实现绘图部分figure (2);subplot (3,1,1);myplot(B, A); %调用绘图函数myplot绘制损耗函数曲线yt- y_l(t);subplot (3, 1,2); tplot (ylt, T, yt) ; %调用绘图函数tplot绘制滤波器输出波形%带通滤波器设计与实现fpl=440;fpu=560;fsl=275;fsu=900;wp二2*fpl/Fs, 2*fpu/Fs;ws=2*fsl/Fs,2*fsu/Fs;rp二0.1;rs=60;N, wp=ellipord(wp, ws, rp, rs) ;%调用ellipord计算椭圆DF阶数N和通带截止频率wpB, A=ellip(N, rp, rs, wp); %调用ellip计算椭圆带通DF系统函数系数向量B和Ay2t二filter (B, A, st) ;%滤波器软件实现%带通滤波器设计与实现绘图部分(省略)%高通滤波器设计与实现fp=890;fs=600;wp=2*fp/Fs;ws=2*fs/Fs;rp=0. l;rs=60;%DF指标(低通滤波器的通、阻带边界频)N, wp=ellipord(wp, ws, rp, rs) ;%调用ellipord计算椭圆DF阶数N和通带截止频率wpB, A=ellip(N, rp, rs, wp, high，) ; %调用ellip计算椭圆带通DF系统函数系数向量B和A y3t=filter (B, A, st) ;%滤波器软件实现%高低通滤波器设计与实现绘图部分(省略)七、实验程序运行结果实验4程序exp4.m运行结果如图104.2所示。由图可见，三个分离滤波器指标参数选 取正确，算耗函数曲线达到所给指标。分离出的三路信号yl(n), y2(n)和y3(n)的波形是抑 制载波的单频调幅波。损耗函数曲线(a)低通滤波器损耗函数及其分离出的调幅信号y (t)i损耗函数曲线00.10.203040.506070.80.91(b)带通滤波器损耗函数及其分离出的调幅信号y2(t)O00 010.Q20 030 040 050.060.07t/（c）高通滤波器损耗函数及其分离出的调幅信号y3（t）图104.实验4程序exp4.m运行结果八、思考题简答（1）请阅读信号产生函数mstg,确定三路调幅信号的载波频率和调制信号频率。（2） 信号产生函数nistg中采样点数N=800,对st进行N点FFT可以得到6根理想谱 线。如果取N二1000,可否得到6根理想谱线？为什么？ N二2000呢？请改变函数mstg中采 样点数N的值，观察频谱图验证您的判断是否正确。（3）修改信号产生函数mstg,给每路调幅信号加入载波成分，产生调幅（AM）信号， 重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别。答：分析发现，st的每个频率成分都是25IIz的整数倍。采样频率Fs=10kHz=25X400IIZ, 即在25Hz的正弦波的1个周期中采样400点。所以，当N为400的整数倍时一定为st的 整数个周期。因此，采样点数N=800和N=2000时，对st进行N点FFT可以得到6根理想谱 线。如果取N=1000,不是400的整数倍，不能得到6根理想谱线Mat lab的IIR滤波器描述清楚iir滤波器的设计建模过程，程序注释，以及对不同设计方案的说明比较等。lirl :低通巴特沃斯模拟滤波器设计。通带截至频率3400 Hz ,通带最大衰减3dB阻带截至频率4000 Hz ,阻带最小衰减40dBIir2 :模拟低通滤波器转换为数字低通滤波器,脉冲响应不变法和双线性变换法。Iir3 :切比雪夫二型低通数字滤波器设计通带边界频率0.2兀,通带最大衰减IdB阻带截至频率0.47T ,阻带最小衰减80dBIir4 :椭圆带通数字滤波器设计Iir5 :高通和带通巴特沃思数字滤波器设计双线性变换%低通巴特沃斯模拟滤波器设计clear; close allfp=3400; fs=4000; Rp=3; As=40;N,fc=buttord(fjp,fs,Rp,Ass)B,A=butter(N,fc,s);hf,f =freqs(B, A, 1024);plot(f,20*log 10(abs(hf)/abs(hf( 1)grid, xlabel(f/Hz); ylabel(幅度(dB) axis(0,4000,40,5);line(0,4000,-3,-3);line(3400,3400,-90,5)%用脉冲响应不变法和双线性变换法将模拟滤波器离散化 clear; close allb=1000;a=l,1000;w=0:1000*2*pi;hf,w=freqs(b9a,w);subplot(2,3,l); plot(w/2/pi,abs(hf); grid;xlabel(f(Hz); ylabel(幅度);title(，模拟滤波器频响特性)FsO= 1000,500;for m=l :2Fs=FsO(m)d,c=impinvar(b,a,Fs)f,e=bilinear(b,a,Fs)wd=0:512 *pi/512;hw 1 =freqz(d,c,wd);h w2=freqz(f,e, wd);subplot(2,3,2); plot(wd/pi,abs(hwl)/abs(hwl(l); grid on; hold on title(，脉冲响应不变法)subplot(2,3,3); plot(wd/pi,abs(hw2)/abs(hw2( 1); grid on; hold on title(，双线性变换法)end%切比雪夫II型低通数字滤波器设计clear; close allwp=0.2; ws=0.4; Rp=l; Rs=80;N,wc=cheb2ord(wp,ws,Rp,Rs)B,A=cheby2(N,Rs,wc)freqz(B,A)%直接设计带通数字椭圆滤波器clear; close allWp=0.25,0.45; Ws=O.15055;Rp=0.1; Rs=60;N,wc=ellipord(Wp,Ws,Rp,Rs)b,a=ellip(N,Rp,Rs,wc)hw,w=freqz(b,a);subplot(2,l,l); plot(w/pi,20*logl0(abs(hw); gridaxis(0,l,80,5); xlabel(w/ 兀)；ylabel(幅度(dB)subplot(2,l,2); plot(w/pi,angle(hw); gridaxis(0,l,-pi,pi); xlabel(*w/ 兀)；ylabel(相位(rad)1)%用双线性变换法设计数字高通和带通滤波器 clear; close allT=l; wch=pi/2;wlc=0.35*pi; wuc=0.65*pi;B=1;A=1,2.6131,3.4142,2.6131,1; h,w=freqz(B,A,512);subplot(2,2,l); plot(w,20*logl0(abs(h); grid%axis(0,10,90,0); xlabel(w/ 兀);title(模拟低通幅度(dB)f) % 高通omegach=2*tan(wch/2)/T;B hs, Ahs=lp2hp(B, A,omegach); Bhz,Ahz=bilinear(Bhs,Ahs, 1/T);h,w=freqz(Bhz, Ahz,512);subplot(2,2,3); plot(w/pi,20*logl0(abs(h); gridaxis(0,1,-150,0); xlabel(w/ n ); title(数字高通幅度(dB) %带通omegalc=2*tan(wlc/2)/T; omegauc=2*tan(wuc/2)/T;wo=sqrt(omegalc*omegauc); Bw=omegauc-omegalc; Bbs,Abs=lp2bp(B,A,wo,Bw);Bbz,Abz=bilinear(Bbs,Abs, 1/T);h,w=freqz(Bbz,Abz,512);subplot(2,2,4); plot(w/pi,20*logl0(abs(h); gridaxis(0,1,-150,0); xlabel(w/ 兀)；title(数字带通幅度(dB)数字滤波器设计与应用问题1. 题目：数字滤波器的设计与应用 设计要求：利用Matlab软件，以复合信号分离为例，对“数字信号处理” 课程中的谱分析、数字滤波器设计和信号滤波这三个过程进行了仿真实现，给 出了仿真结果。