无限长单位脉冲响应IIR滤波器的设计方法MATLAB

  • 格式:doc
  • 大小:1.06 MB
  • 文档页数:18

如有帮助,欢迎支持。

1 数字信号处理实验指导

实验四、 无限长单位脉冲响应(IIR)滤波器的设计方法

(一) 实验目的

加深对无限冲激响应( IIR )数字滤波器的常用指标和设计过程的理解。

(二) 实验内容

常用函数介绍:

1、Matlab信号处理工具箱中提供了设计巴特沃思模拟滤波器的函数buttord、buttap和butter,格式如下:

(1)[,](,,,,CPSPSNWbuttordWWRRs

用于计算巴特沃思模拟低通滤波器的阶N和3dB截止频率Wc (即本书中的符号c)。其中,Wp 和Ws分别是滤波器的通带截止频率p和阻止截止频率s,单位为rad/s;Rp和Rs分别是通带最大衰减系数p和阻带最小衰减系数s,单位为dB。

(2)[,,]()zpGbuttapN

用于计算N阶巴特沃思归一化(c=1)模拟低通滤波器系统函数的零、极点和增益因子,返回长度为N的向量z和p分别给出N个零点和极点,G是滤波器增益。得到的滤波器系统函数形式如下:

1212()()()()()()()()()aNaaNQsszszszHsGPsspspsp

其中,kz和kp分别是向量z和p的第k个元素。如果要从零、极点得到系统函数的分子和分母多项式系数向量B和A,可以调用结构转换函数

(3)[,]2(,,)BAzptfzpG,结构转换后系统函数的形式为

111111()()()MMMaNNNbsbsbBsHsAsasasa

其中,M 是向量B的长度,N是向量A的长度,kkba和分别是向量B和A的第k个元素。 如有帮助,欢迎支持。

2 (3)[,](,,,)CBAbutterNWftypes

用于计算巴特沃思模拟滤波器系统函数中分子和分母多项式系数向量B和A,其中N和CW分别是滤波器的阶和3dB截止频率c,返回向量B和A中的元素ka和kb分别是上面的()aHs表示式中的分母和分子系数。ftype缺省时,设计低通滤波器;ftype=high时,设计高通滤波器;ftype=stop时,设计带阻滤波器,此时CW为向量[1,]ccuWW,且ftype缺省时,设计带通滤波器,带通的频率区间为1ccuWW。S缺省时,设计数字滤波器。

例如:设计一个满足下列指标要求的巴特沃思模拟低通滤波器。指标

sraddBsraddsspp/4,30,/2,1B

Matlab程序:

2;4;1;30[,](,,,,[,](,,)00.0016[](pspsCPSPSCWWRRNWbuttordWWRRsBAbutterNWsWHWfreqs; %滤波器参数; %计算滤波器的阶和3dB截止频率; %计算滤波器系统函数分子、分母多项式系数::; %在0~6rad/s频率范围内取1001个频率点,,,);20*log10(());(,),(//);(/)BAWWHabsHBplotWHgridxlabelradsylabeldB %计算频率向量上的滤波器频率响应 %纵坐标的单位取d on %绘制幅频响应曲线频率(幅度

运行结果如下:

N=6

Wc=2.2496

B= 0 0 0 0 0 0 129.5917

A=1.0000 8.6916; 37.7720; 104.0667;191.1447 ;222.5973 ; 129.5917

即:

65432129.5917()8.691637.7720104.0667191.1447222.5793129.5917aHsssssss幅频响应曲线如图6.3.4所示。

2、Matlab中设计数字滤波器的函数都是采用双线性变换法,将模拟滤波器转换为数字滤波器。这些函数及其凋用格式如下(巴特沃思数字滤波器):

(1)[,](,,,)CPSPSNWbuttordWWRR

该格式用于计算巴特沃思数字滤波器的阶N和3dB截止频率的归一化值CW如有帮助,欢迎支持。

3 (关于归一化)。调用参数Wp 和Ws分别是数字滤波器的通带截止频率和阻带截止频率的归一化值(关于归一化),要求01pW和01sW,其中1表示数字频率(对应模拟频率/2sf,sf为采样频率)。pR和sR分别是通带最大衰减和阻带最小衰减,单位dB。当spWW时,设计高通滤波器;当pW和sW是二元向量时,设计带通(spWW)或带阻(spWW)滤波器,这时返回参数CW也是二元向量。

(2)[,](,,)CBAbutterNWftype;

该格式用于计算巴特沃思模拟滤波器系统函数中分子和分母多项式系数向量B和A。调用参数N 和WC分别是巴特沃思数字滤波器的阶和3dB截止频率归一化值(关于归一化)。当ftype缺省时,设计低通滤波器;当ftype=high时,设计高通滤波器;ftype=stop时,设计带阻滤波器,此时CW为二元向量[1,]ccuWW,1cW和cuW分别是带阻滤波器的通带3dB下,上截止频率的归一化值(关于归一化);CW为向量[1,]ccuWW,且ftype缺省时,设计带通滤波器,带通的频率区间为1ccuWW。注意设计出的带通和带阻数字滤波器是2N阶的,这是因为带通滤波器可表示为一个N阶低通滤波器与一个N阶高通滤波器的级联。

由函数的返回向量B和A可写出数字滤波器的系统函数为:

1112111121()NNNNNNbbzbzbzHzaazazaz

其中,kkba和分别是向量B和A的第k个元素。

例如:用双线性变换法设计巴特沃思数字低通滤波器,指标要求

通带截止频率0.2pwrad,通带最大衰减1pdB;

阻带截止频率0.35swrad,阻带最小衰减10sdB

解 设计步骤如下:

(1)给出数字滤波器的指标。

0.2;0.35;10;PSSWWR

(2)计算巴特沃思数字滤波器的阶N和3dB截止频率。

[,](,,,)CPSPSNWbuttordWWRR; 如有帮助,欢迎支持。

4 (3)用双线性变换法设计巴特沃思数字低通滤波器。

[,](,)ZZCBAbutterNW;

运行结果如下:

ZB=0.0335 0.1006 0.1006 0.0335

ZA=1.000 -1.4245 0.8827 -1.1900

设计出的数字滤波器系统函数为

1231230.03350.10060.10060.0335()1.0001.42450.88270.1900zzzHzzzz

程序如下:

0.2;0.35;10;10;[,](,,,);[,](,);PSPSCPSPSZZCWWRRNWbuttordWWRRBAbutterNW %滤波器参数 %计算数字滤波器的阶和3dB截止频率 %计算数字滤波器系统函数分子和分母多项 式的系数向量B和A%省去绘图部分程序

例:利用脉冲响应不变法,把系统函数为651)(2ssssHa的模拟滤波器变换成等价的数字滤波器,采样间隔ST1.0

手算:2132651)(2ssssssHa

21112136065.05595.1108966.01.0112)(zzzzeTzeTzHTT

Matlab:

B=[1,1];

A=[1,5,6];

T=0.1;

Fs=1/T;

[Bz,Az]=impinvar(B,A,Fs);%用脉冲响应不变法将模拟滤波器变换成数字滤波器

运行结果: 如有帮助,欢迎支持。

5 Bz=0.1000 -0.0897

Az=1.000 -1.5595 0.6065

1、 设采样周期T=250μs(采样频率fs =4kHz),用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器,其3dB边界频率为fc =1kHz。

[B,A]=butter(3,2*pi*1000,'s');

[num1,den1]=impinvar(B,A,4000);

[h1,w]=freqz(num1,den1);

[B,A]=butter(3,2/0.00025,'s');

[num2,den2]=bilinear(B,A,4000);

[h2,w]=freqz(num2,den2);

f=w/pi*2000;

plot(f,abs(h1),'-.',f,abs(h2),'-');

grid;

xlabel('频率/Hz ')

ylabel('幅值/dB')

程序中第一个butter的边界频率2π×1000,为脉冲响应不变法原型低通滤波器的边界频率;第二个butter的边界频率2/T=2/0.00025,为双线性变换法原型低通滤波器的边界频率.图3.1给出了这两种设计方法所得到的频响,虚线为脉冲响应不变法的结果;实线为双线性变换法的结果。脉冲响应不变法由于混叠效应,使得过渡带和阻带的衰减特性变差,并且不存在传输零点。同时,也看到双线性变换法,在z=-1即ω=π或f=2000Hz处有一个三阶传输零点,这个三阶零点正是模拟滤波器在Ω=∞处的三阶传输零点通过映射形成的。 如有帮助,欢迎支持。

6

2、 设计一数字高通滤波器,它的通带为400~500Hz,通带内容许有0.5dB的波动,阻带内衰减在小于317Hz的频带内至少为19dB,采样频率为1,000Hz。

wc=2*1000*tan(2*pi*400/(2*1000));

wt=2*1000*tan(2*pi*317/(2*1000));

[N,wn]=cheb1ord(wc,wt,0.5,19,'s');

[B,A]=cheby1(N,0.5,wn,'high','s');

[num,den]=bilinear(B,A,1000);

[h,w]=freqz(num,den);

f=w/pi*500;

plot(f,20*log10(abs(h)));

axis([0,500,-80,10]);

grid;