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

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

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

3、以分离 st 中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。且滤波器的通带最大衰减为 0.1dB，阻带最小衰减为 60bB。提示：抑制载波单频调幅信号的数学表示式为 ttttts ffff ccc 000 2os2os12osc)( 其中， 称为载波， 为载波频率， 称为单频调制信号，c2o为调制正弦波信号频率，且满足 。由上式可见，所谓抑制载波单频f0 fc0调幅信号，就是 2 个正弦信号相乘，它有 2 个频率成分：和频 + 和差频fc0- ，这 2 个频率成分关于载波频率 对称。所以，1 路抑制载波单频调fc0 fc幅信号的频谱图

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

5、减 ；Hzfp280BPd1.0阻带截止频率 ，阻带最小衰减 ;s456s 对于载波频率为 500Hz 的条幅信号，可以选用带通滤波器分离，其指标为:通带截止频率 ， ，通带最大衰减 ；zfpl0zfpu50d1.0p阻带截止频率 ， ，阻带最小衰减 ；Hsl275s9B6s 对于载波频率为 1000Hz 的条幅信号，可以选用高通滤波器分离，其指标为：通带截止频率 ，通带最大衰减 ； zfp890BPd1.0阻带截止频率 ，阻带最小衰减 ;Hs56s说明：（1）为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽可能使滤波器过渡带宽尽可能宽；（2）与信号产生函数 mstg 相同采用频率

6、Fs=10kHz；（3）为了滤波器阶数最低，选用椭圆滤波器。4试验程序框图调用函数 mstg 产生 st，自动绘图显示 st 的时域波形和幅频特性曲线调用 ellipord 和 ellip 分别设计三个椭圆滤波器，并绘图显示其幅频响应特性曲线调用 filter，用三个滤波器分别对信号 st 进行滤波，分离出三路不同载波频率的调幅信号 y1(n),y2(n),和 y3(n)绘图显示 y1(n),y2(n)和 y3(n)的时域波形结束25思考题及简答1、阅读信号产生函数 mstg，确定三路调幅信号的载波频率和调制信号频率答：第一路调幅信号的载波频率 fc1=1000Hz； 第一路调幅信号的调制频率

7、 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 的每个频率成分都是 25Hz 的整数倍。采样频率Fs=10kHz=25400

8、Hz，即在 25Hz 的正弦波的 1 个周期中采样 400 点。所以，当N 为 400 的整数倍时一定为 st 的整数个周期。因此，采样点数 N=1600 和N=2000 时，对 st 进行 N 点 FFT 可以得到 6 根理想谱线。如果取 N=1800，不是400 的整数倍，不能得到 6 根理想谱线。（1）N=1600 时：0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02-10123t/ss(t)(a)s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51

9、(b) s(t)函 函 函f/HZ函函（2）N=1800 时：30 0.005 0.01 0.015 0.02-10123t/ss(t)(a)s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(b) s(t)函 函 函f/HZ函函（3）N=2000S 时：0 0.005 0.01 0.015 0.02 0.025-10123t/ss(t)(a)s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(b) s(t)函 函 函f/HZ函函3、修改信号产生函数

10、mstg，给每路调幅信号加入载波成分，产生调幅（AM）信号，重复本实验，观察 AM 信号与抑制载波调幅信号的时域波形及其频谱的差别。4AM 信号表示式： ttts ffAcmd 2oscs0 Amd取值： ， ，结果见（附录）10dA5m6实验结果（程序附录）原信号输出：0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01-10123t/ss(t)(a)s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(b) s(t)函 函 函f/HZ函函低通输出：0 0

11、.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-60-40-200/函函函dB函(a) s(t)函 函 函 函 函 函0 0.01 0.02 0.03 0.04 0.05 0.06 0.07-0.500.51t/sy(n)(b) s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(c) s(t)函 函 函f/Hz函函5带通输出：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-60-40-200/函函函dB函(a) s(t)函 函 函 函 函 函0 0.0

12、1 0.02 0.03 0.04 0.05 0.06 0.07-0.500.51t/sy(n)(b) s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(c) s(t)函 函 函f/Hz函函高通输出：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-60-40-200/函函函dB函(a) s(t)函 函 函 函 函 函0 0.01 0.02 0.03 0.04 0.05 0.06 0.07-0.500.51t/sy(n)(b) s(t)函 函 函0 200 400 600 800 100

13、0 1200 1400 1600 1800 200000.51(c) s(t)函 函 函f/Hz函函6附录（）：原信号输出：0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01-2002040t/ss(t)(a)s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(b) s(t)函 函 函f/HZ函函低通输出：70 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-60-40-200/函函函dB函(a) s(t)函 函 函 函

14、 函 函0 0.01 0.02 0.03 0.04 0.05 0.06 0.07-10010t/sy(n)(b) s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(c) s(t)函 函 函f/Hz函函带通输出：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-60-40-200/函函函dB函(a) s(t)函 函 函 函 函 函0 0.01 0.02 0.03 0.04 0.05 0.06 0.07-10010t/sy(n)(b) s(t)函 函 函0 200 400 600 800

15、1000 1200 1400 1600 1800 200000.51(c) s(t)函 函 函f/Hz函函高通输出：80 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-60-40-200/函函函dB函(a) s(t)函 函 函 函 函 函0 0.01 0.02 0.03 0.04 0.05 0.06 0.07-10010t/sy(n)(b) s(t)函 函 函0 200 400 600 800 1000 1200 1400 1600 1800 200000.51(c) s(t)函 函 函f/Hz函函附录（）：%主函数%IIR 数字滤波器设计及软件实现clea

16、r all;%调用信号产生函数 mstg 产生又三路抑制载波调幅信号相加构成的复合信号syms st;syms t;st=mstg; %低通滤波器设计与实现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

17、 系统函数系数向量 B 和 Ay1t=filter(B,A,st);%滤波器的软件实现fyt=fft(y1t,n);%下面为绘图部分figure(2);subplot(3,1,1);myplot(B,A);9yt=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=90

18、0;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,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/

19、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(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);10subplot(3,1,2);tplot(y1t,T,yt);subplot(3,1,3);stem(f,abs(fyt)/max(a

20、bs(fyt),.);grid;title(c) s(t)的频谱);axis(0,Fs/5,0,1.2);xlabel(f/Hz);ylabel(幅度);clc;clear%子程序%产生信号程序function 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;%第二路调幅信号的调制信号频

21、率为 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 的频谱.%绘图subplot(2,1,1)plot(t,st);grid on;xlabel(t/s);ylabel(s(t);axis(

22、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);11plot(W/pi,20*log10(m/max(m);grid on;xlabel(omega/pi);ylabel(幅度（dB）);axis(0,1,-80,5);title(a) s(t)损耗函数曲线);%滤波器输出波形function tplot(xn,T,yn)n=0:length(xn)-1;t=n*T;plot(t,xn);grid on;xlabel(t/s);ylabel(y(n);axis(0,t(end),min(xn),1.2*max(xn)title(b) s(t)的波形);