IIR数字滤波器的设计

一、目的和要求

1.掌握脉冲响应不变法和双线性变换法设计IR数字滤波器的原理和方法:
2.观察双线性变换法和脉冲响应不变法设计的滤波器的频域特性,了解双线性变换法和脉冲响应不变法的特点和区别。

二、原理和内容

使用MATILAB编写程序,实现IR数字滤波器的设计。涉及脉冲响应不变法和双线性变换法设计IR数字滤波器的方法、不同设计方法得到的IR滤波器频域特性异同等知识点。

1.脉冲响应不变法
所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h（n）等干模拟滤波器 的单位冲激响应和ha（t）的采样值，即:h(n)=ha(t)|t=nT=ha(nt),其中，T为采样周期在脉冲响应不变法中,模拟角频率和数字角频率的变换关系为: w=ΩT,可见, Ω和w 之间的变换关系为线性的.
在MATLAB中,可用函数impivar实现从模拟滤波器到数字滤波器的脉冲响应不变映射,调用格式为:
[B,A]=impinvar (b, a, fs1)
[B,A]=impinvar (b, a)
其中,b、a分别为模拟滤波器的分子和分母多项式系数向量; fs1为采样频率(H2),缺省值fs-IHz; B、 A分别为数字滤波器分子和分母多项式系数向量。
2.双线性变换法: 由于s平面和2平面的单值双线性映射关系为s-21其中工为采材周T142
期,因此,若已知模拟滤波器的传递函数,将上式代人即可得到数字地波器的系统百散Hz).
在双线性变换中,模拟和频率和数字角频率的变换关系为: 2-31见, Q和 之间的变换关系为非线性的。
在MATLAB中,可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射,调用格式为: [B.A]-bilinearb.a. 511
3.数字滤波器设计
(1)定技术指标转换为模拟滤波器设计性能指标。
(2)估计满足性能指标的模拟相应滤波器性能阶数和截止频率。
利用MATLAB中buttord, cheblord, cheb2ord, ellipord等函数,
调用格式如: [N,9c] = buttord(Sp, 2s, cp, as,'s)
其中, Sp为通带边界频率, rads: 2s为阻带边界频率, rad/s: a p为带通波动, dB;as为阻带衰减, dB: 's'表示为模拟滤波器;函数返回值N为模拟滤波器的最小阶数; 3c为模拟滤波器的截止频率(3dB频率), rads。函数适用低通、高通、带通、带阻滤波器。
(3)设计模拟滤波器。
MATLAB信号处理工具箱提供了模拟滤波器设计的完全工具函数: butterchebyl, chebyz2、 ellip. besself.用户只需一次调用就可完成低通、高通、带通.
带阻滤波器设计
调用格式如: [ba]-btter(N. .'ftype,'s),其中, 'type
为滤波器类型:
high'表示高通滤波器,截止频率为8c
stop表示带阻波器,Q.c-[21 a2] (Q1< 2)
"tper缺省时表示为低通或带通滤波器
(4)利用脉冲响应不变法或双线性不变法,实现模拟滤波器到数字滤波器的映射

三、步骤

在“开始一程序”菜单中,找到MATLAB程序,运行启动;
进入MATLAB后,在Command Window中输入自己编写的主程序,并执行;记录运行结果图形,作分析对比。
具体步骤如下
1、查看帮助文件,了解相关函数的调用格式。
2、设计模拟巴特斯低通滤波器, fp=skHz, ap=2dB, fs=12kHz, a s=30dB.
3、设计模拟巴特沃斯高通滤波器, fp=200Hz, fs=100Hz,幅度特性单调下降, fp处最大衰减为34B,阻带最小衰减as=15dB.
4. Fp=0. 2kHZ, a p=ldB, fs=0. 4kHZ, a s=30daB, 采样间隔T=1ms:分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器(1)观察所设计数字滤波器的幅频特性曲线,记录带宽和衰减量;
(2)比较两种方法的优缺点:
常,酸t
5、设计一数字高通滤波器,它的通带为700~ 100012,通带内容许有2dB的波动,阻带内衰减在小于500Hz的频带内至少为30dB,采样频率为2000Hz
6、利用MATLAB编程设计一个Butterworth数字带通滤波器,指标要求如下:通带边缘频率: wpl=0.5irad, wp2=0.7mrad,通带峰值起伏: ap<=2dB;
阻带边缘频率: w

ws1=0.35Πrad, ws2=0.95Πrad, 最小阻带衰减: a s>=35dB.分别用脉冲响应不变法和双线性变换法进行IR数字滤波器的设计。T=0.2ms

四、结果与分析

1.

ap=2;
as=30;
fp=5000;
fs=12000;
omegap=2*pi*fp;
omegas=2*pi*fs;
[N,omegac]=buttord(omegap,omegas,ap,as,'s');
[b,a]=butter(N,omegac,'s');
[H,omega]=freqs(b,a);
plot(omega/2/pi,20*log10(abs(H)));
axis([0,20000,-50,0]);
xlabel('频率Hz’)
ylabel('H幅值dB‘)
title(‘巴特沃斯低通滤波器幅频特性‘) 

2.

ap=3;
as=15;
fp=200;
fs=100;
omegap=2*pi*fp;
omegas=2*pi*fs;
[N,omegac]=buttord(omegap,omegas,ap,as,'s');
[b,a]=butter(N,omegac,'high','s');
[H,omega]=freqs(b,a);
plot(omega/2/pi,20*log10(abs(H)));
axis([0,1000,-50,0]);
xlabel('频率Hz‘)
ylabel('H幅值dB’)
title('巴特沃斯高通滤波器幅频特性‘) 

3.

 ap=1;
as=30;
fp=200;
fs=400；
omegap=2*pi*fp;
omegas=2*pi*fs;
fs1=1000;
[N,omegac]=buttord(omegap,omegas,ap,as,'s');
[b,a]=butter(N,omegac,'s');
[B,A]=impinvar(b,a,fs1)
[H1,w]=freqz(B,A,'whole');
subplot(211);plot(w*fs1/2/pi,20*log10(abs(H1)));
axis([0,1000,-100,0]);
xlabel('频率Hz‘)
ylabel('H1幅值dB‘)
ap=1;
as=30;
fp=200;
fs=400;
fs1=1000;
wp=2*pi*fp/fs1;
ws=2*pi*fs/fs1;
omegap=2*fs1*tan(wp/2);
omegas=2*fs1*tan(ws/2);
[N,omegac]=buttord(omegap,omegas,ap,as,'s');
[b,a]=butter(N,omegac,'s');
[B,A]=bilinear(b,a,fs1);
[H2,w]=freqz(B,A,'whole');
subplot(212);plot(w*fs1/2/pi,20*log10(abs(H2)));
axis([0,1000,-100,0]);
xlabel('频率Hz‘)
ylabel('H2幅值dB‘) 

4.

omegap=2*2000*tan(2*pi*700/(2*2000));
omegas=2*2000*tan(2*pi*500/(2*2000));
[N,omegac]=cheb1ord(omegap,omegas,2,30,'s');
[b,a]=cheby1(N,0.5,omegac,'high','s');
[B,A]=bilinear(b,a,1000);
[H,w]=freqz(B,A);
f=w/pi*500;
plot(f,20*log10(abs(H)));
axis([0,1000,-80,10]);
grid;
xlabel('频率Hz‘)；
ylabel('幅值/dB‘)；

5.

fs=5000;
omegas1=2*fs*tan(0.35*pi/2);
omegas2=2*fs*tan(0.95*pi/2);
omegap1=2*fs*tan(0.5*pi/2);
omegap2=2*fs*tan(0.7*pi/2);
omegas=[omegas1,omegas2];
omegap=[omegap1,omegap2];
ap=2;as=35;
[N,omegac]=buttord(omegap,omegas,ap,as,'s');
[b,a]=butter(N,omegac,'s');
[B,A]=bilinear(b,a,fs);
[H,w]=freqz(B,A);
f=w/2/pi*fs;
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('相位‘)；
figure(2)
fs=5000;
omegas1=0.35*pi*fs;
omegas2=0.95*pi*fs;
omegap1=0.5*pi*fs;
omegap2=0.7*pi*fs;
omegas=[omegas1,omegas2];
omegap=[omegap1,omegap2];
ap=2;as=35;
[N,omegac]=buttord(omegap,omegas,ap,as,'s');
[b,a]=butter(N,omegac,'s');
[B,A]=impinvar(b,a,5000);
[H,w]=freqz(B,A);
f=w/pi*2500;
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('相位‘)；