MATLAB讲义第五六章.docx

MATLAB讲义第五六章.docx

《MATLAB讲义第五六章.docx》由会员分享，可在线阅读，更多相关《MATLAB讲义第五六章.docx（17页珍藏版）》请在冰豆网上搜索。

MATLAB讲义第五六章

第五、六章　IIR滤波器的设计

滤波器结构：

一、系统函数的表示法及其转换

1、表示方法

（1）传递函数法

若

　则a=[1a

（2）a（3）…a（N）]

b=[b

（1）b

（2）…b（M）]

（2）零极点增益法

若

则　零点向量　Z=[z1z2zM-1];

极点向量　P=[z1,z2,…,zN-1]

k为系统增益。

（3）部分分式法

　若

则极点向量p=[p

（1）p

（2）…p（n）]

其对应系数向量　r=[r

（1）r

（2）…r（n）]

余数多项式系数向量　k=[k

（1）k

（2）…k（M-N+1）]

（4）二阶分式法

把H（z）划成二阶因式

则其二阶因式为：

　　b01b11b211a11a21

b02b12b221a21a22

Sos=…

b0Nb1Nb2N1a1Na2N

2、各表示方法的转换：

（1）由传递函数转换为零极点增益（tf2zp）

　对应，由零极点增益转换为传递函数（zp2tf）

调用方法：

[z,p,k]=tf2zp（b,a）

[b,a]=zp2tf（z,p,k）

――a,b的长度要相等，不等的话要补零

（2）由零极点增益转换为二次分式（zp2sos）

　对应，由二次分式转换为零极点增益（sos2zp）

　调用方法：

[sos,g]=zp2sos（z,p,k）,g为整个系统的增益，即H（z）=g*H1（z）*H2（z）…*HN（z）

调用方法：

sos2zp（z,p,k）=[sos,g],g为整个系统的增益,默认为1。

（3）二次分式转换为传递函数（sos2tf）

调用方法：

[b,a]=sos2tf（sos）

对应，由传递函数转换为二次分式（tf2sos），调用方式：

[sos,g]=tf2sos（b,a）

例1：

>>b=[1350];

>>a=[1813];

>>[z,p,k]=tf2zp（b,a）

z=

0

-1.5000+1.6583i

-1.5000-1.6583i

p=

-7.9216

-0.0392+0.6141i

-0.0392-0.6141i

k=

1

例2：

差分方程：

16y（n）+12y（n-1）+2y（n-2）-4y（n-3）-y（n-4）=x（n）-3x（n-1）+11x（n-2）-27x（n-3）+18y（n-4）

程序：

>>b=[1-311-2718];

>>a=[16122-4-1];

>>[sos,g]=tf2sos（b,a）

sos=

1.0000-3.00002.00001.0000-0.2500-0.1250

1.00000.00009.00001.00001.00000.5000

g=

0.0625

则其级联结构方程为：

由此可画出其结构图。

例3：

H（z）=1+16.0625z-4+z-8

>>b=[100016.06250001];

>>[sos,G]=tf2sos（b,1）

sos=

1.0000-2.82844.00001.000000

1.00002.82844.00001.000000

1.0000-0.70710.25001.000000

1.00000.70710.25001.000000

G=

1

其传递函数应写成：

H（z）=（1+2.83z-1+4z-2）（1-2.83z-1+4z-2）（1+0.71z-1+0.25z-2）（1－0.71z-1+0.25z-2）

滤波器的设计：

一、模拟滤波器的设计

1、巴氏模拟原型滤波器的设计：

用buttap函数，用于设计N阶归一化（Ωc=1）巴氏模拟低通滤波器。

调用格式：

　[z0,p0,k0]=buttap（N）

其中：

N――为巴氏滤波器阶数

　　　z0,p0,0k――为归一化滤波器的零点、极点、增益

所设计的N阶巴氏模拟原型滤波器的归一化传递函数为：

若要设计未归一化（Ωc≠1）巴氏低通滤波器，则要用Ωc乘以p0、z0，而分子k0乘以Ω

因为极点有N个。

例1：

设计阶数为3，5，10，15的巴氏模拟原型滤波器。

并画出幅频响应曲线。

程序

fori=1:

4;

switchi

case1

N=3;

case2

N=5;

case3

N=10;

case4;

N=15;

end;

[z,p,k]=buttap（N）;

[b,a]=zp2tf（z,p,k）;

[h,w]=freqs（b,a,n）;

Ah=abs（h）;

subplot（2,2,i）,plot（w,Ah）;axis（[0201]）;xlabel（'w/wc'）;ylabel（'|H（jw）|.^2'）;

title（['filerN=',num2str（N）]）;grid;end;

2、巴氏滤波器的最小阶数的求法：

巴氏模拟原型滤波器的阶数越大，响应曲线在通带内越平缓，在阻带内的衰减越大。

如何选择一个合适的阶数，使其刚好能实现系统的最优性能，而实现不复杂。

阶数选择函数buttord：

调用方式：

[N,OmegaC]=buttord（OmegaP,OmegaS,Rp,Rs，’s’）

N――返回的滤波器最小阶数

Wc――3dB频率

OmegaP,OmegaS――通带、阻带截止频率，为归一化频率，范围为[0,1]，对应π。

Rp,Rs――通带、阻带内的衰减

s――表示设计的是模拟滤波器

例1：

设通带、阻带截止频率fp=5kHz、fs=12kHz，通带、阻带最大衰减Rp=1dB,Rs=30dB，要求设计巴氏低通滤波器。

>>OmegaP=2*pi*5000;OmegaS=2*pi*12000;

>>Rp=1;Rs=30;

>>[N,OmegaC]=buttord（OmegaP,OmegaS,Rp,Rs,'s'）%确定阶数N

N=

5

OmegaC=

3.7792e+004

>>[z0,p0,k0]=buttap（N）%确定传递函数

z0=

[]

p0=

-0.3090+0.9511i

-0.3090-0.9511i

-0.8090+0.5878i

-0.8090-0.5878i

-1.0000

k0=

1

>>p=p0*OmegaC;%去归一化

>>k=k0*OmegaC^N

k=

7.7094e+022

>>p=p0*OmegaC

p=

1.0e+004*

-1.1678+3.5943i

-1.1678-3.5943i

-3.0575+2.2214i

-3.0575-2.2214i

-3.7792

>>b=k*real（poly（z））

b=

7.7094e+022

>>a=real（poly（p））

a=

1.0e+022*

0.00000.00000.00000.00000.00077.7094

由结果可见，因为Ωc相当大，Ω

达到1022数量级，致使向量a中的其它系数无法显示。

因此最好仍用归一化系数表示滤波器传递函数。

即归一化传递函数为：

>>b0=real（poly（z0））

b0=

1

>>a0=real（poly（p0））

a0=

1.00003.23615.23615.23613.23611.0000

为检验此滤波器在Ωp、Ωs处的幅频特性，可用freqs函数，调用格式：

　　　　H＝freqs（b,a,w）

b,a――滤波器分子、分母系数向量

　　　　　w――自变量的频率向量，不允许取标量，故至少要同时取两个频点

若b、a取归一化值b0、a0，w也应归一化为w/OmegaC，即

　　　　　H=freqs（b0,a0,w/OmegaC）

本例中，取OmegaP,OmegaS作为w,即w=[OmegaP,OmegaS]，并对其归一化为w/OmegaC

>>b0=real（poly（z0））

b0=

1

>>a0=real（poly（p0））

a0=

1.00003.23615.23615.23613.23611.0000

>>w=[5000,12000]*2*pi/37792;

>>H=freqs（b0,a0,w）;

键入freqs（b,a），就可得出系统的频率特性

二、IIR数字滤波器的设计

1、冲激响应不变法

采用impinvar函数

[bz,az]=impinvar（b,a,fs,tol）

bz,az――转换的分子分母向量

　　fs――采样频率，默认值为1Hz

tol――误差容限，表示转换后的离散系统函数是否有重复的极点。

实现思路：

先把Ha（s）的分子分母向量ba,aa利用residue函数转换为极点留数Ra,pa及直接项Ga,再把模拟滤波器的极点映射为数字极点pd=epa*T，得到H（z）的极点留数形式，再用residuez把H（z）转换为传递函数形式。

function[bd,ad]=impinvar_my（ba,aa,Fs）

[Ra,pa,Ga]=residue（ba,aa）;%把模拟滤波器系数向量变为极点留数形式

T=1/Fs;pd=exp（pa*T）;Rd=Ra*T%将模拟极点留数转换为数字极点留数

[bd,ad]=residuez（Rd,pd,Ga）%将极点留数形式转换为传递函数形式

bd=real（bd）,ad=real（ad）

例1：

利用脉冲响应不变法，把

转换为数字滤波器，T＝0.1。

>>b=[1,1];

>>a=[156];

>>T=0.1;

>>[bz,az]=impinvar（b,a,1/T）

bz=

0.1000-0.0897

得到H（z）为：

az=

1.0000-1.55950.6065

>>subplot（2,1,1）,t=0:

0.1:

3;

>>hat=impulse（b,a,t）;%计算模拟系统的脉冲响应ha（t）,并乘以T

>>plot（t,hat*T,'red'）,holdon

>>hnT=impz（bz,az,31）;%计算数字系统的脉冲响应ha（nT）,并以同样时间轴作图

>>stem（0.1*[0:

30],hnT,'g'）

>>subplot（2,1,2）,w=[0:

0.1:

10]*2*pi;

>>HaS=freqs（b,a,w）;plot（w,abs（HaS）,'red'）,holdon%模拟系统的频率特性Ha（s）

>>Hz=freqz（bz,az,w/Fs）;plot（w/Fs,abs（Hz）,'g'）%数字系统的频率特性H（z）

>>plot（w,abs（Hz）,'blu'）%将数字频率放大Fs倍，与模拟频率的幅频特性相比较（∵w=2π*Fs,范围为0～2π，所以要与模拟的相比，必须乘以Fs）

2、双线性变换法

　　采用bilinear

　　[zd,pd,kd]=bilinear（z,p,k,fs）

――把模拟滤波器的零极点模型转换为数字的零极点模型。

[bz,az]=bilinear（bs,as,fs）

――把模拟滤波器的传递函数转换为数字的

例：

，T＝1和0.1,用双线性变换法设计一个巴氏数字低通滤波器。

>>bs=[11];

>>as=[156];

>>T=0.1;Fs=1/T;

>>[bz,az]=bilinear（bs,as,Fs）

bz=

T＝0.1时

0.04150.0040-0.0375

az=

1.0000-1.55730.6047

T=1;Fs=1/T;

[bz,az]=bilinear（ba,aa,Fs）

bz=

T＝1时

0.15000.1000-0.0500

az=

1.00000.2000-0.0000

>>subplot（2,1,1）,t=0:

0.1:

3;

>>hat=impulse（bs,as,t）;plot（t,hat*T,'red'）,holdon

>>hn=impz（bz,az,31）;stem（0.1*[0:

30],hn,'*','g'）

>>title（'ha（t）andh（n）'）;

>>subplot（2,1,2）,w=[0:

0.1:

10]*2*pi;

>>HaS=freqs（bs,as,w）;plot（w,abs（HaS）,'red'）,holdon

>>Hz=freqz（bz,az,w/Fs）;plot（w/Fs,abs（Hz）,'g'）

>>plot（w,abs（Hz）,'blu'）

>>title（'Ha（s）andH（z）'）;

>>