完整word版用MATLAB设计滤波器.docx

完整word版用MATLAB设计滤波器.docx

《完整word版用MATLAB设计滤波器.docx》由会员分享，可在线阅读，更多相关《完整word版用MATLAB设计滤波器.docx（13页珍藏版）》请在冰豆网上搜索。

完整word版用MATLAB设计滤波器

用MATLAB设计滤波器

1IIR滤波器的设计

●freqz

功能：

数字滤波器的频率响应。

格式：

[h，w］=freqz（b，a,n）

[h，f]=freqz（b，a，n，Fs）

［h，w］=freqz（b，a,n，’whole'）

［h，f］=freqz（b,a，n，'whole’,Fs）

h=freqz（b，a，w）

h=freqz（b,a，f，Fs）

freqz（b，a）

说明:

freqz用于计算由矢量"和b构成的数字滤波器

H（z）=

=

的复频响应H（jω）.

［h，w]=freqz（b,a，n）可得到数字滤波器的n点的幅频响应,这n个点均匀地分布在

上半单位圆（即0~π），并将这n点频率记录在w中，相应的频率响应记录在h中。

至于n

值的选择没有太多的限制，只要n〉0的整数，但最好能选取2的幂次方，这样就可采用

FFT算法进行快速计算。

如果缺省，则n=512。

[h，f］二freqz（b,a,n，Fs）允许指定采样终止频率Fs（以Hz为单位），也即在0~Fs/2

频率范围内选取n个频率点（记录在f中）,并计算相应的频率响应h。

［h,w]=freqz（b，a，n，’whole’）表示在0~2π之间均匀选取n个点计算频率响应.

［h，f］=freqz（b，a，n，'whole'，Fs）则在O～Fs之间均匀选取n个点计算频率响应.

h=freqz（b，a，w）计算在矢量w中指定的频率处的频率响应，但必须注意，指定的频

率必须介于0和2π之间.

h=freqz（b，a,f，Fs）计算在矢量f中指定的频率处的频率响应,但指定频率必须介于

0和Fs之间。

●butter

功能：

Butterworth（比特沃思）模拟和数字滤波器设计。

格式：

有六种

[b,a］=butter（n，Wn）

[b,a］=butter（n，Wn，’ftype’）

［b,a]=butter（n，Wn，’s’）

［b,a]=butter（n，Wn，'ftyPe'，’s'）

[z，p，k］=butter（…）

［A，B,C，D]=butter（…）

说明:

butter函数可设计低通、带通、高通和带阻的数字和模拟滤波器，其特性为使通带内

的幅度响应最大限度地平坦，这会损失截止频率处的下降斜度.在期望通带平滑的情况

下,可使用butter函数,但在期望下降斜度大的场合,应使用椭圆和Chebyshev（切比雪夫）

滤波器。

butter函数可设计出数字域和模拟域的Butterworth滤波器。

（1）数字域

[b，a]=butter（n，Wn）可设计出截止频率为Wn的n阶低通Butterworth滤波器，其

滤波器为

H（z）=

=

截止频率是滤波器幅度响应下降至1/

处的频率，Wn∈[0，1］，其中1相应于0。

5fs（取样

频率，即Nyquist频率）。

当Wn=［WlW2]（Wl〈W2）时，butter函数产生一2n阶的数字带通滤波器，其通带

为Wl<ω

［b，a]=butter（n，Wn，'ftype’）可设计出高通或带阻滤波器:

·当ftype=high时,可设计出截止频率为Wn的高通滤波器；

·当ftype=stop时，可设计出带阻滤波器,这时Wn=［WlW2］，且阻带为Wl〈ω〈W2

利用输出变量个数的不同，可得到滤波器的另外两种表示：

零极点增益和状态方程。

·[z，p,k]=butter（n，Wn）或［z,p，k]=butter（n，Wn,'ftype'）可得到滤波器的零

极点增益表示；

·[A,B，C，D］=butter（n,Wn）或［A，B,C，D]=butter（n，Wn，'ftype’）可得到滤

波器的状态空间表示。

（2）模拟域

[b，a]=butter（n,Wn，’s’）可设计截止频率为Wn的n阶低通模拟Butterworth滤波

器为

H（s）=

=

其中截止频率Wn〉0。

模拟域的butter函数说明与数字域的完全相同，它也有六种形式。

例1：

设数据采样频率为900Hz，现要设计一9阶的高通Butterworth滤波器,截止频率为300Hz。

％ex104.m

[b，a］=butter（9,300/500，'high’）；

freqz（b,a，128，1000）

例2:

设计一10阶的带通Butterworth滤波器，带通为100--200Hz，并画出滤波器的冲击响应.

％ex105.m

n=5；

Wn=［100200］/500；

[b，a］=butter（n，Wn）；

figure

（1）；

freqz（b，a，128,1000）

figure

（2）；

impz（b,a,101）

例3利用双线性变换设计数字Butterworth滤波器

源程序:

Wp=0。

2＊pi；

Ws=0。

3＊pi；

Rp=1；

As=15；

T=1;

Fs=1/T;

OmegaP=（2/T）＊tan（Wp/2）;

OmegaS=（2/T）*tan（Ws/2）；

ep=sqrt（10^（Rp/10）—1）；

Ripple=sqrt（1/（1+ep*ep））;

Attn=1/10^（As/20）;

[cs,ds]=adf_butt（OmegaP,OmegaS,Rp，As）；

［b,a］=bilinear（cs,ds，T）;

[C,B,A]=sdir2cas（b，a）

［cs，ds］=s521（OmegaP,OmegaS,Rp,As）；％调用1

[b,a]=bilinear（cs,ds，T）；

［b0,B,A]=s522（b,a）；%调用2

%

％plot

figure

（1）;

[db,mag,pha,grd，w]=freqz_m（b，a）；

subplot（2，2，1）；plot（w/pi，mag）；

title（'MagnitudeResponse'）；

xlabel（’frequencyinpiunits'）;ylabel（’|H|'）;axis（［0，1，0，1］）；

set（gca，’XTickMode’，'manual’，'XTick’，［0，0。

2,0。

3,1］）；

set（gca,’YTickMode’，’manual'，'YTick’，［0,Attn,Ripple,1］）;grid；

subplot（2，2，3）;plot（w/pi，db）；

title（’MagnitudeindB’）；

xlabel（’frequencyinpiunits’）；ylabel（'decibels'）；

axis（［0,1，-40,5］）；

set（gca，’XTickMode','manual’，'XTick'，[0，0.2，0。

3，1]）;

set（gca,'YTickMode'，’manual','YTick’，［-50,-15，—1,0]）；grid；

%set（gca,'YTickLabelMode','manual’，'YTickLabels’,['50’；'15’;'1’；'0’]）；

subplot（2,2，2）;plot（w/pi，pha/pi）；title（’PhaseResponse’）

xlabel（’frequencyinpiunits’）；ylabel（’piunits'）;axis（［0，1，-1，1]）;

set（gca,’XTickMode’，'manual’，’XTick',［0，0.2,0。

3,1]）;

set（gca，’YTickMode’,’manual’,'YTick',［—1，0,1］）;grid

subplot（2,2,4）;plot（w/pi，grd）；title（'GroupDelay'）；

xlabel（’frequencyinpiunits'）；ylabel（'samples’）;axis（［0，1，0，10］）;

set（gca,’XTickMode','manual'，'XTick’,[0，0。

2，0.3,1]）;

set（gca,’YTickMode’,’manual'，’YTick',［0：

2：

10］）;grid

子程序1:

freqz的修正

%freqz_m.m

function［db,mag，pha,grd，w]=freqz_m（b，a）;

[H,w］=freqz（b,a,1000,'whole'）；

H=（H（1:

501））'；w=（w（1:

501））’；

mag=abs（H）;

db=20＊log10（（mag+eps）/max（mag））；

pha=angle（H）;

grd=grpdelay（b,a，w）;

子程序2

变滤波器直接形式为级联形式

%sdir2cas。

m

function［C,B,A]=sdir2cas（b，a）

％DIRECT—formtoCASCADE-formconversinins-plane

%-——-————-———--———---———--——-—————-————--—

%［C，B，A］=sdir2cas（b,a）

%C=gaincoefficient

%B=Kby3matrixofrealcoefficientscontainingbk's

%A=Kby3matrixofrealcoefficientscontainingak’s

%b=numberatorpolynomialcoefficientsofDIRECTform

％a=denominatorpolynomialcoefficientsofDIRECTform

％

Na=length（a）-1；Nb=length（b）-1；

%computegaincoefficientC

b0=b

（1）;b=b/b0;

a0=a

（1）；a=a/a0;

C=b0/a0；

％

％Debnominatorsecond—ordersections:

p=cplxpair（roots（a））；K=floor（Na/2）；

ifK＊2==Na%ComputationwhenNaiseven

A=zeros（K,3）；

forn=1：

2:

Na

Arow=p（n：

1：

n+1，：

）；

Arow=poly（Arow）；

A（fix（（n+1）/2），：

）=real（Arow）;

end

elseifNa==1%ComputationwhenNa=1

A=[0real（poly（p））];

else%ComputationwhenNaisoddand>1

A=zeros（K+1,3）;

forn=1:

2:

2*K

Arow=p（n:

1：

n+1，：

）;

Arow=poly（Arow）;

A（fix（n+1）/2,：

）=real（Arow）；

end

A（K+1，:

）=[0real（poly（p（Na）））]；

end

%numeratorsecond—ordersections：

z=cplxpair（roots（b））;K=floor（Nb/2）；

ifNb==0%computationwhenNb=0

B=[00poly（z）］；

elseifK＊2==Nb%computationwhenNbiseven

B=zeros（K，3）；

forn=1:

2:

Nb

Brow=z（n：

1:

n+1,:

）

Brow=poly（Brow）;

B（fix（（n+1）/2），:

）=real（Brow）；

end

elseifNb==1％computationwhenNb=1

B=［0real（poly（z））］;

else%ComputationwhenNbisoddand>1

B=zeros（K+1，3）；

forn=1:

2：

2*K

Brow=z（n:

1：

n+1,:

）

Brow=poly（Brow）;

B（fix（（n+1）/2）,：

）=real（Brow）；

end

B（K+1，：

）=［0real（poly（z（Nb）））］；

end

●chebyl

功能：

Chebyshev（切比雪夫）I型滤波器设计（通带等波纹）。

格式：

[b，b]=chebyl（n，Rp，Wn）

［b，a］=chebyl（n，Rp,Wn，’ftype'）

［b，a］=chebyl（n,Rp，Wn,'s’）

[b，a］=chebyl（n，Rp,Wn,'ftype'，'s’）

［z,p，k］=chebyl（…）

［A，B，C,D]=chebyl（…）

说明:

chebyl函数可设计低通、带通、高通和带阻的数字和模拟ChebyshevI型滤波器，其

通带内为等波纹，阻带内为单调.ChebyshevI型滤波器的下降斜度比Ⅱ型大，但其代价

是在通带内波纹较大。

与butter函数类似,chebyI函数可设计出数字域和模拟域的ChebyshevI型滤波器.

（1）数字域

[b，a］=chebyI（n，Rp,Wn）可设计出n阶低通数字ChebyshevI型滤波器.其截止频

率由Wn确定,通带内的波纹由Rp（分贝）确定,滤波器为

H（z）=

=

截止频率是滤波器幅度下降至—Rp分贝处的频率，对chebyI函数,Wn∈［0，1］，Wn=l

时相应于0。

5fs。

通带波纹Rp越小,可得到更宽的变换宽度.

当Wn=［WlW2]时，chebyI函数可产生一2n阶的数字带通滤波器，其通带为Wl<ω〈W2

[b,a]=chebyI（n,Rp，Wn，'ftype’）可设计高通或带阻滤波器：

·当ftype=high时，可设计出截止频率为Wn的高通滤波器;

·当ftype=stop时，可设计出带阻滤波器，这时Wn=［WlW2］，且阻带为Wl〈ω〈W2。

利用输出变量个数的不同，可得到滤波器的另外两种表示：

零极点增益和状态方程。

·［z，p，k]=chebyI（n,Rp，Wn）或［z，p,k］=chebyI（n，Rp，Wn，’ftype’）可得到滤

波器的零极点增益表示;

·[A，B，C，D］=chebyI（n，Rp,Wn）或[A，B，C,D］=chebyI（n,Rp,Wn，’ftype'）

可得到滤波器的状态空间表示。

（2）模拟域

［b，a]=chebyI（n，Rp，Wn，'s’）可设计出截止频率为Wn的n阶低通模拟ChebyshevI型滤波器

H（s）=

=

其中截止频率Wn>0。

模拟域的chebyl函数说明与数字域完全相同，它也有六种形式。

例4：

设数据采样频率为900Hz，现要设计一9阶的低通Butterworth滤波器，截止频率为300Hz.

并设定波纹系数为0.5dB。

其程序如下

%ex106.m

［b，a］=butter（9，300/500）;％butterworth

figure

（1）;

freqz（b,a，128，1000）

[b，a]=cheby1（9，0.5,300/500）;％chebyshev

figure

（2）;

freqz（b，a，512，1000）

2椭圆滤波器的设计

若信号由5Hz、15Hz、30Hz三个正弦频率成分构成.设计一个椭圆滤波器,滤除5Hz和30Hz频率成分。

％hz15.m

Fs=100;

t=（0：

99）/Fs；

s1=sin（2*pi＊t＊5）；s2=sin（2*pi*t*15）；s3=sin（2＊pi＊t*30）;

s=s1+s2+s3；

[b,a]=ellip（4，0.1,40，[1020]＊2/Fs）；

［H,w]=freqz（b,a,512）;

S=fft（s,512）;

sf=filter（b，a，s）；

SF=fft（sf，512）；

subplot（3,2,1）；plot（t，s）;title（’Timedomain’）;

grid

ylabel（'originalsignal’）；

subplot（3,2,5）；plot（t，sf）;

xlabel（'Time（s）’）；

ylabel（'filtedsignal'）；

w1=（0：

255）/256*（Fs/2）;

subplot（3,2，2）;plot（w1,abs（S（1：

256）））;title（'Frequencydomain'）；

grid

ylabel（'originalspectrum’）;

subplot（3,2,4）;plot（w*Fs/（2＊pi），abs（H））;

ylabel（’filter’）;

grid;

subplot（3,2，6）;plot（w1,abs（SF（1:

256）））;

xlabel（’Frequency（Hz）’）;

ylabel（'filtedspectrum’）;

grid

3利用kaiser窗函数法设计FIRDF

%利用kaiser窗函数,设计长度为45的带阻滤波器,使阻带衰减为60dB.

％Hd=10<=w

％Hd=0pi/3〈=w<=2＊pi/3

%aa。

m

As=60;

M=45;％或M=ceil（（As—7。

95）/（2.286＊deltaw））+1

n=[0:

1：

M-1];

beta=0。

1102*（As-8.7）％若验证不满足指标要求，要增大beta或M,修正设计

w_kai=（kaiser（M,beta））'；

wc1=pi/3；wc2=2＊pi/3；

hd=ideal_lp（wc1,M）+ideal_lp（pi，M）-ideal_lp（wc2,M）;

h=hd。

＊w_kai；

［db，mag，pha，grd，w]=afreq（h，[1］）;

figure

（1）;

subplot（2，2,1）;stem（n，hd）;title（'IdealImpulseResponse'）;

axis（［-1M—0.20.8］）;xlabel（’n'）；ylabel（'hd（n）’）;

subplot（2，2,2）；stem（n，w_kai）;title（'KaiserWindow'）;

axis（[—1M01。

1]）;xlabel（’n’）;ylabel（'w（n）'）;

subplot（2,2，3）;stem（n，h）；title（’ActualImpulseResponse’）；

axis（［—1M-0.20.8]）；xlabel（’n'）;ylabel（’h（n）'）;

subplot（2,2,4）;plot（w/pi，db）;title（’MagnitudeResponseindB’）；grid;

xlabel（'frequencyinpiunits’）;ylabel（'Eecibels'）；

axis（［01-8010］）;

％子程序1-ideal_lp

functionhd=ideal_lp（wc，M）;

alpha=（M-1）/2;

n=0：

1：

M-1；

m=n-alpha+eps;

hd=sin（wc＊m）。

/（pi*m）；

%子程序2-afreq。

m

function［db,mag，pha，grd，w]=afreq（b，a）;

［H，w]=freqz（b，a，1000,’whole'）；

H=（H（1：

1：

501））’；w=（w（1:

1：

501））';

mag=abs（H）;

db=20＊log10（（mag+eps）/max（mag））;

pha=angle（H）;

grd=grpdelay（b，a，w）；