IIR数字滤波器的设计及软件实现

1、IIR 数字滤波器的设计及软件实现一实验目的（ 1）熟悉用双线性变换法设计IIR 数字滤波器的原理与方法；（ 2）学会用MATLAB 信号处理工具箱中的滤波器设计函数 （或滤波器设计分析工具 FDAtool ）设计各种滤波器，学会根据滤波需求确定滤波器指标参数；（ 3）掌握IIR 数字滤波器的 MATLAB 实现方法；（ 4）通过观察滤波器输入、输出信号的时域波形及其频谱，建立数字滤波的概念。二实验原理设计 IIR 数字滤波器一般采用间接法（脉冲响应不变法和双线性不变法），应用最广泛的是双线性变换法。基本的设计过程是：将给定的数字滤波器指标转换成模拟滤波器的指标；涉及模拟滤波器；将模拟滤波器的

2、系统函数转换成数字滤波器的系统函数。MATLAB 信号处理工具箱中的各种 IIR 数字滤波器设计函数都是采用双线性变换法。 本实验的数字滤波器的 MATLAB 实验是调用 MATLAB 信号处理工具箱的函数filter对给定的输入信号x（n）进行滤波，得到滤波后的输出信号y（n）.三 实验内容及步骤信号处产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号 st ，该函数还会自动回图显示st 的时域波形和幅频特性曲线，由后图可见，三路信号时域混叠无法在时域分离。 但频域是分离的， 所以可通过滤波的方法在 频域分离。.将st中三路调幅信号分离，通过观察 st的幅频特性曲线，分别确定可以

3、分离st中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器)的通带截止频率和阻带截止频率。且滤波器的通带最大衰减为0.1dB,阻带最小衰减为 60bB。提示：抑制载波单频调幅信号的数学表示式为s(t) = cos2二 f ot cos2二 f t =g Cos2二 f c - f 0 t 1 cos2二 f c f。，其中，cos(2nf t端为载波，f为载波频率，cos2nf t麻为单频调制信号， cc0f 0为调制正弦波信号频率，且满足 f c f 0。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频f c+ f 0和差频 f c-

4、 f ,这2个频率成分关于载波频率f c对称。所以，1路抑制载波单频调幅 信号的频谱图是关于载波频率 fc对称的2根谱线，其中没有载频成分，故取名为抑制载波单频调幅信号。图中三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。.编程调用MATLAB滤波器涉及函数ellipord和ellip分别设计这三个椭圆滤波器，并绘图显示其损耗函数曲线；.调用滤波器实验函数filter,用三个滤波器分别对信号产生函数 mstg产生的信号st进行滤波，分离出st中的三路不同载波频率的调幅信号 y1(n)、y2(n)、y3(n)。滤波器参数的选取：对于载波频率为250Hz的条幅信号，可以选用低通滤

5、波器分离，其指标为：通带截止频率f = 280Hz ,通带最大衰减Fp = 0.1dB ； p阻带截止频率f = 450Hz,阻带最小衰减6 6 6 60dB ; ss对于载波频率为500Hz的条幅信号，可以选用带通滤波器分离，其指标为 TOC o 1-5 h z 通带截止频率f = 440Hz , f =560Hz,通带最大衰减0 n = 0.1dB ; plpup阻带截止频率f =275Hz, f =900Hz ,阻带最小衰减a s = 60dB ; slsus对于载波频率为1000Hz的条幅信号，可以选用高通滤波器分离，其指标为:通带截止频率f p= 890Hz ,通带最大衰减eP =

6、0.1dB ;阻带截止频率f = 550Hz,阻带最小衰减 氤=60dB; ss说明：(1)为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽可能使滤波器过渡带宽尽可能宽；(2)与信号产生函数mstg相同采用频率Fs=10kHz ;(3)为了滤波器阶数最低，选用椭圆滤波器。四.试验程序框图五思考题及简答1 、阅读信号产生函数mstg ，确定三路调幅信号的载波频率和调制信号频率答：第一路调幅信号的载波频率fc1=1000Hz ；第一路调幅信号的调制频率fm1=100Hz ； TOC o 1-5 h z 第二路调幅信号的载波频率fc2=500Hz；第二路调幅信号的调制频率fm2=50Hz；

7、第三路调幅信号的载波频率fc3=250Hz；第三路调幅信号的调制频率fm3=25Hz；2、信号产生函数mstg中采样点数N=1600 ,对st进彳T N点FFT就可以得到6根理想谱线。如果取N=1800 ，可否得到 6 根理想谱线？为什么？ N=2000 呢？请改变采样点数N 的值，观察频谱图验证判断是否正确？答：因为信号st 是周期序列，谱分析时要求观察时间为整数倍周期。分析可知，st的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25 400Hz ,即在 25Hz 的正弦波的 1 个周期中采样400 点。所以，当 N 为 400 的整数倍时一定为 st 的整数个周期。因此，采

8、样点数N=1600 和 N=2000 时，对 st 进行 N点 FFT 可以得到 6 根理想谱线。如果取N=1800 ，不是 400 的整数倍，不能得到6根理想谱线。(1) N=1600 时:(a)s(t)的波形t/s(b) s(t)的频谱f/HZ(2) N=1800 时:3(a)s(t)的波形00.0050.010.0150.02t/s(b) s(t)的频谱(3) N=2000S 时:(a)s(t)的波形t/s(b) s(t)的频谱f/HZ3、修改信号产生函数 mstg,给每路调幅信号加入载波成分，产生调幅(AM)信号，重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差 别。A

9、M 信号表示式：st )= Ad + Amcos2 兀 f 0tLs2K f ct) Ad 2 Am取化Ad =10 , Am =5,结果见(附录i)六.实验结果(程序附录n)原信号输出：1JJ1*1/ 11.(a)s(t)的波形210-100.0010.0020.0030.0040.0050.0060.0070.0080.0090.01t/s(b) s(t)的频谱f/HZ低通输出:B -20(-40区-60幅-80 C=0.10.20.30.40.50.60.70.80.91“二(b) s(t)的波形(c) s(t)的频谱带通输出:(a) s(t)损耗函数曲线1 0.50 -0.5Lfl 1

10、 1f1 / 11 1 ./ 100.010.020.030.040.050.060.07t/ss(t)的频谱f/Hz高通输出:B -20(-40 度-60 幅-8001 0.50-0.51Y1_L0.10.20.30.40.(b) s(t)白50.6勺波形0.70.80.9/fl 1 1f工1 J1 1 1 f二1 r? 11U00.010.020.030.040.050.060.07t/s s(t)的频谱1度幅0.50012004006008001000120014001600180020f/Hz0074-4 IM z/l 录B 附原(a)s(t)的波形t/s(b) s(t)的频谱f/HZ

11、低通输出:(a) s(t)损耗函数曲线1度幅111A A:八; J-A l l 1. i/X F a八八1 .1 J 7，/ 1 I J 1, I1i1串i f i /j* /1 j 1 j-10 0-1000.010.020.030.040.050.060.07t/ss(t)的频谱f/Hz度 幅带通输出：1度幅t/ss(t)的频谱高通输出:3 1度幅度幅附录(n)：% 主函数%IIR 数字滤波器设计及软件实现clear all;% 调用信号产生函数mstg 产生又三路抑制载波调幅信号相加构成的复合信号syms st;syms t;st=mstg; %低通滤波器设计与实现Fs=10000;T=

12、1/Fs;n=800;Tp=n*T;k=0:n-1;f=k/Tp;fp=280;fs=450;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF 指标 ;(低通滤波器的通阻带边界频率)N,wp0=ellipord(wp,ws,rp,rs);% 调用 ellipod 计算椭圆 DF 阶数 N 和通带截止频率wpB,A=ellip(N,rp,rs,wp0);% 调用 ellip 计算椭圆带通DF 系统函数系数向量B 和 Ay1t=filter(B,A,st);% 滤波器的软件实现fyt=fft(y1t,n);% 下面为绘图部分figure(2);subplot(3,1,

13、1);myplot(B,A);yt=y_1(t);tplot(y1t,T,yt);subplot(3,1,3);stem(f,abs(fyt)/max(abs(fyt),.);grid;title(c) s(t) 的频谱 );axis(0,Fs/5,0,1.2);xlabel(f/Hz);ylabel( 幅度 );% 带通滤波器的实现与设计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,wp0=ellipord(wp,ws,rp,rs);B,A=ellip(N,rp

14、,rs,wp0);y2t=filter(B,A,st);fyt=fft(y2t,n);figure(3);subplot(3,1,1);myplot(B,A);yt=y_1(t);tplot(y1t,T,yt);subplot(3,1,3);stem(f,abs(fyt)/max(abs(fyt),.);grid;title(c) s(t) 的频谱 );axis(0,Fs/5,0,1.2);xlabel(f/Hz);ylabel( 幅度 );% 高通滤波器的实现与设计fp=890;fs=600;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;N,wp0=ellipord(

15、wp,ws,rp,rs);B,A=ellip(N,rp,rs,wp0,high);y3t=filter(B,A,st);fyt=fft(y3t,n);figure(4);subplot(3,1,1);myplot(B,A);yt=y_1(t);subplot(3,1,2);tplot(y1t,T,yt);subplot(3,1,3);stem(f,abs(fyt)/max(abs(fyt),.);grid;title(c) s(t) 的频谱 );axis(0,Fs/5,0,1.2);xlabel(f/Hz);ylabel( 幅度 );clc;clear% 子程序% 产生信号程序function

16、 st=mstgN=800FS=10000;T=1/FS;TP=N*T;t=0:T:(N-1)*T;K=0:N-1;f=K/TP;fc1=FS/10;%第一路调幅信号的载波频率fc1=1000HZfm1=fc1/10;%第一路调幅信号的调制信号频率为fm1=100hz.fc2=FS/20;%第二路调幅信号的载波频率fc2=500HZfm2=fc2/10;%第二路调幅信号的调制信号频率为fm2=50hz.fc3=FS/40;%第三路调幅信号的载波频率fc3=250HZfm3=fc3/10;%第三路调幅信号的调制信号频率为fm3=25hz.xt1=cos(2*pi*fm1*t).*cos(2*pi

17、*fc1*t);xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);st=xt1+xt2+xt3;fxt=fft(st,N);% 计算信号 st 的频谱 .%绘图subplot(2,1,1)plot(t,st);grid on;xlabel(t/s);ylabel(s(t);axis(0,TP/8,min(st),max(st);title(a)s(t) 的波形 )subplot(2,1,2)stem(f,abs(fxt)/max(abs(fxt),.);grid on;title(b) s(t) 的频谱 )axis(0,FS/5,0,1.2);xlabel(f/HZ);ylabel( 幅度 );% 损耗输出波形function myplot(B,A)H,W=freqz(B,A,1000);m=abs(H);plot(W/pi,20*log10(m/