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

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

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

3、波器、高通滤波器的通带截止频率和阻带截止频率。且滤波器的通带最大衰减为0.1dB,阻带最小衰减为60bB。提示：抑制载波单频调幅信号的数学表示式为,1s(t)cos2ftcos2ftcos2fftcos2fft0c2c0c0其中，cos2ft称为载波，f为载波频率，cos2ft称为单频调制信号，cc0f0为调制正弦波信号频率，且满足fc>f0。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频fc+f0和差频fc-f0,这2个频率成分关于载波频率fc对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率fc对称的2根谱线，其中没有载频成分，故取名为抑制

4、载波单频调幅信号。图中三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。3 .编程调用MATLAB滤波器涉及函数ellipord和ellip分别设计这三个椭圆滤波器，并绘图显示其损耗函数曲线；4 .调用滤波器实验函数filter,用三个滤波器分别对信号产生函数mstg产生的信号st进展滤波，别离出st中的三路不同载波频率的调幅信号y1(n卜y2(ny3(n)。滤波器参数的选取：对于载波频率为250Hz的条幅信号，可以选用低通滤波器别离，其指标为：通带截止频率f280Hz,通带最大衰减p0.1dB;p阻带截止频率f450Hz,阻带最小衰减，60dB;ss对于载波频率为500Hz的

5、条幅信号，可以选用带通滤波器别离，其指标为通带截止频率f440Hz,f560Hz，通带最大衰减D0.1dB;plpup阻带截止频率f275Hz,f900Hz,阻带最小衰减60dB;slsus对于载波频率为1000Hz的条幅信号，可以选用高通滤波器别离，其指标为:通带截止频率fp890Hz,通带最大衰减P0.1dB;阻带截止频率f550Hz,阻带最小衰减，60dB;ss说明：1为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原那么是尽可能使滤波器过渡带宽尽可能宽；(2)与信号产生函数mstg一样采用频率Fs=10kHz;(3)为了滤波器阶数最低，选用椭圆滤波器。四.试验程序框图1、阅读信号产生

6、函数mstg,确定三路调幅信号的载波频率和调制信号频率第一路调幅信号的载波频率fc1=1000Hz；第一路调幅信号的调制频率fm1=100Hz；第二路调幅信号的载波频率fc2=500Hz；第二路调幅信号的调制频率fm2=50Hz；第三路调幅信号的载波频率fc3=250Hz；第三路调幅信号的调制频率fm3=25Hz；2、信号产生函数mstg中采样点数N=1600,对st进展N点FFT就可以得到6根理想谱线。如果取N=1800，可否得到6根理想谱线？为什么？N=2000呢？请改变采样点数N的值，观察频谱图验证判断是否正确？答：因为信号st是周期序列，谱分析时要求观察时间为整数倍周期。分析可知，st

7、的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25X400Hz,即在25Hz的正弦波的1个周期中采样400点。所以，当N为400的整数倍时一定为st的整数个周期。因此，采样点数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时:(a)s(t)的波形t/s(b)s(t)的频谱度幅f/HZ(3)N=2000S时:(a)s(t)的波形f/HZ3、修改信号产生函数mstg,给每路调幅信号参加载波成分，产生调幅

8、AM信号,重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差异。AM信号表示式：stAdAm8s2f0tss2fctAdAm取值：Ad10,Am5,结果见附录i六.实验结果程序附录n原信号输出：(a)s(t)的波形3211110-100.0010.0020.0030.0040.0050.0060.0070.0080.0090.01t/s(b)s(t)的频谱f/HZ低通输出:(a) s(t)损耗函数曲线(b) s(t)的波形呷rI1m(c)s(t)的频谱10.500200400600800100012001400160018002000f/Hz带通输出:XDa 1度幅00.010.

9、020.030.040.050.060.07t/ss(t)的频谱mmTfr110.500200400600800100012001400160018002000f/Hz高通输出:度幅附录I：原信号输出:(a)s(t)的波形1度幅0.50 0200400600800100012001400160018002000t/sf/HZ低通输出：XDa 1度幅io0-ioJL1A Ajry! I!1 f A""111 ； jy V Hi,1:J'二,一0.010.020.030.040.050.060.07t/ss(t)的频谱11110.50020040060080010001

10、2001400160018002000f/Hz带通输出：(a)s(t)损耗函数曲线XDa 1度幅1A/iin!1IiL八i11S11ijf八；j1.iJ7jpU'I；iJ17Pii/I,:LJ,.-100-1000.010.020.030.040.050.060.071度幅0.50t/ss(t)的频谱0200400600800100012001400160018002000f/Hz高通输出:XDa l度幅附录n：%主函数%IIR数字滤波器设计及软件实现clearall;%调用信号产生函数mstg产生又三路抑制载波调幅信号相加构成的复合信号symsst;symst;st=mstg;%低通

11、滤波器设计与实现Fs=10000;T=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

12、,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('幅度');%带通滤波器的实现与设计fpl=440;fpu=560;fsl=275;fsu=900;wp=2*fpl/Fs,2*fpu/Fs;ws=2*fsl/Fs,2*fs

13、u/Fs;rp=0.1;rs=60;N,wp0=ellipord(wp,ws,rp,rs);B,A=ellip(N,rp,rs,wp0);y2t=filter(B,A,st);fyt=fft(y2t,n);figure(3);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);x

14、label('f/Hz');ylabel('幅度');%高通滤波器的实现与设计fp=890;fs=600;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;N,wp0=ellipord(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

15、,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%子程序%产生信号程序functionst=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=

16、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*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的频谱.%绘图subp

17、lot(2,1,1)plot(t,st);gridon;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),'.');gridon;title('(b)s(t)的频谱')axis(0,FS/5,0,1.2);xlabel('f/HZ');ylabel('幅度');%损耗输出波形functionmyplot(B,A)H,W=freqz(B,A,1000);m=abs(H);plot(W/pi,20*log10(m/ma