数字滤波器的设计及其MATLAB实现.docx
《数字滤波器的设计及其MATLAB实现.docx》由会员分享，可在线阅读，更多相关《数字滤波器的设计及其MATLAB实现.docx（29页珍藏版）》请在冰豆网上搜索。
数字滤波器的设计及其MATLAB实现
设计低通数字滤波器，要求在通带内频率低于0.2pirad时，允许幅度误差在1dB以内，
在频率0.3pirad～pirad之间的阻带衰减大于15dB，用脉冲响应不变法设计数字滤波器，T=1:
切比雪夫I型模拟滤波器的设计子程序：
function[b,a]=afd_chb1（Omegap,Omegar,Ar）
ifOmegap<=0
error（'通带边缘必须大于0'）
end
if（Dt<=0）|（Ar<0）
error（'通带波动或阻带衰减必须大于0'）;
end
ep=sqrt（10^（Dt/10）-1）;
A=10^（Ar/20）;
OmegaC=Omegap;
OmegaR=Omegar/Omegap;
g=sqrt（A*A-1）/ep;
N=ceil（log10（g+sqrt（g*g-1））/log10（OmegaR+sqrt（OmegaR*OmegaR-1）））;
fprintf（'\n***切比雪夫I型模拟低通滤波器阶数=%2.0f\n',N）;
[b,a]=u_chblap（N,Dt,OmegaC）;
设计非归一化切比雪夫I型模拟低通滤波器原型程序：
function[b,a]=u_chblap（N,Dt,OmegaC）
[z,p,k]=cheb1ap（N,Dt）;
a=real（poly（p））;
aNn=a（N+1）;
p=p*OmegaC;
a=real（poly（p））;
aNu=a（N+1）;
k=k*aNu/aNn;
b0=k;
B=real（poly（z））;
b=k*B;
直接形式转换成级联形式子程序：
function[C,B,A]=sdir2cas（b,a）
Na=length（a）-1;
Nb=length（b）-1;
b0=b
（1）;b=b/b0;
a0=a
（1）;a=a/a0;
C=b0/a0;
p=cplxpair（roots（a））;K=floor（Na/2）;
ifK*2==Na
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
A=[0real（poly（p））];
else
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
z=cplxpair（roots（b））;K=floor（Nb/2）;
ifNb==0
B=[00poly（z）];
elseifK*2==Nb
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
B=[0real（poly（z））];
else
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
计算系统函数的幅度响应和相位响应子程序：
function[db,mag,pha,w]=freqs_m（b,a,wmax）
w1=0:
500;
w=w1*wmax/500;
h=freqs（b,a,w）;
mag=abs（h）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（h）;
脉冲响应不变法程序：
function[b,a]=imp_invr（c,d,T）
[R,p,k]=residue（c,d）;
p=exp（p*T）;
[b,a]=residuez（R,p,k）;
b=real（b）.*T;
a=real（a）;
数字滤波器响应子程序：
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）;
直接转换成并联型子程序：
function[C,B,A]=dir2par（b,a）
M=length（b）;
N=length（a）;
[r1,p1,C]=residuez（b,a）;
p=cplxpair（p1,10000000*eps）;
x=cplxcomp（p1,p）;
r=r1（x）;
K=floor（N/2）;
B=zeros（K,2）;
A=zeros（K,3）;
ifK*2==N
fori=1:
2:
N-2
br=r（i:
1:
i+1,:
）;
ar=p（i:
1:
i+1,:
）;
[br,ar]=residuez（br,ar,[]）;
B（（fix（i+1）/2）,:
）real（br'）;
A（（fix（i+1）/2）,:
）real（ar'）;
end
[br,ar]=residuez（r（N-1）,p（N-1）,[]）;
B（K,:
）=[real（br'）0];
A（K,:
）=[real（ar'）0];
else
fori=1:
2:
N-1
br=r（i:
1:
i+1,:
）;
ar=p（i:
1:
i+1,:
）;
[br,ar]=residuez（br,ar,[]）;
B（（fix（i+1）/2）,:
）real（br）;
A（（fix（i+1）/2）,:
）real（ar）;
end
End
比较两个含同样标量元素但（可能）有不同下标的复数对及其相位留数向量子程序：
functionI=cplxcomp（p1,p2）
I=[];
fori=1:
length（p2）
forj=1:
length（p1）
if（abs（p1（j）-p2（i））<0.0001）
I=[I,j];
end
end
end
I=I';
双线性变换巴特沃斯低通滤波器设计：
巴特沃思模拟滤波器的设计子程序：
function[b,a]=afd_butt（wp,ws,Rp,rs）
ifwp<=0
error（'通带边缘必须大于0'）;
end
ifws<=wp
error（'阻带边缘必须大于通带边缘'）;
end
if（Rp<=0）|（Rs<0）
error（'通带波动或阻带衰减必须大于0'）;
end
N=ceil（（log10（（10^（Rp/10）-1）/（10^（Rs/10）-1）））/（2*log10（wp/ws）））;
fprintf（'\n***ButterworthFilterOrder=%2.0f\n',N）;
OmegaC=wp/（（10^（Rp/10）-1）^（1/（2*N）））;
[b,a]=u_buttap（N,OmegaC）
设计非归一化巴特沃思模拟低通滤波器原型子程序：
function[b,a]=u_buttap（N,OmegaC）
[z,p,k]=buttap（N）;
p=p*OmegaC;
k=k*OmegaC^N;
B=real（poly（z））;
b0=k;
b=k*B;
a=real（poly（p））;
直接型到级联型形式的转换：
function[b0,B,A]=dir2cas（b,a）
b0=b
（1）;
b=b/b0;
a0=a
（1）;
a=a/a0;
b0=b0/a0;
M=length（b）;
N=length（a）;
ifN>M
b=[b,zeros（1,N-M）];
elseifNa=[a,zeros（1,M-N）];
else
NM=0;
end
k=floor（N/2）;
B=zeros（k,3）;
A=zeros（k,3）;
ifk*2==N
b=[b,0];
a=[a,0];
end
broots=cplxpair（roots（b））;
aroots=cplxpair（roots（a））;
fori=1:
2:
2*k
br=broots（i:
1:
i+1,:
）;
br=real（polt（br））;
B（（fix（i+1）/2）,:
）=br;
ar=aroots（i:
1:
i+1,:
）;
ar=real（polt（ar））;
A（（fix（i+1）/2）,:
）=ar;
End
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=grdelay（b,a,w）;
设计一个巴特沃思高通滤波器，要求通带截止频率为0.6pi，通带内衰减不大于1dB，阻带·起始频率为0.4pi，阻带内衰减不小于15dB，T=1：
>>wp=0.6*pi;ws=0.4*pi;
>>Rp=1;Rs=15;T=1;
>>[N,wn]=buttord（wp/pi,ws/pi,Rp,Rs）计算巴特沃思滤波器阶数和截止频率
N=
4
wn=
0.5344
>>[b,a]=butter（N,wn,'high'）;频率变换法计算巴特沃思高通滤波器
>>[C,B,A]=dir2cas（b,a）
C=
0.0751
B=
1.0000-2.00001.0000
1.0000-2.00001.0000
A=
1.00000.15620.4488
1.00000.11240.0425
>>[db,mag,pha,grd,w]=freqz_m（b,a）;
>>subplot（2,1,1）;plot（w/pi,mag）;
>>subplot（2,1,2）;plot（w/pi,db）;
椭圆带通滤波器的设计--ellip函数的应用：
>>ws=[0.3*pi0.75*pi];数字阻带边缘频率
>>wp=[0.4*pi0.6*pi];数字通带边缘频率
>>Rp=1;Rs=40;
>>Ripple=10^（-Rp/20）;通带波动
>>Attn=10^（-Rs/20）;阻带衰减
>>[N,wn]=ellipord（wp/pi,ws/pi,Rp,Rs）计算椭圆滤波器参数
N=
4
wn=
0.40000.6000
>>[b,a]=ellip（N,Rp,Rs,wn）;数字椭圆滤波器的设计
>>[b0,B,A]=dir2cas（b,a）级联形式实现
b0=
0.0197
B=
1.00001.50661.0000
1.00000.92681.0000
1.0000-0.92681.0000
1.0000-1.50661.0000
A=
1.00000.59630.9399
1.00000.27740.7929
1.0000-0.27740.7929
1.0000-0.59630.9399
>>figure
（1）;
>>[db,mag,pha,grd,w]=freqz_m（b,a）;
>>subplot（2,2,1）;plot（w/pi,mag）;
>>gridon;
>>subplot（2,2,3）;plot（w/pi,db）;gridon;
>>subplot（2,2,2）;plot（w/pi,pha/pi）;gridon;
>>subplot（2,2,4）;plot（w/pi,grd）;
设计一个巴特沃思带阻滤波器，要求通带上下截止频率为0.8pi、0.2pi，通带内衰减不大于1dB，阻带上起始频率为0.7pi、0.4pi，阻带内衰减不小于30dB：
>>wp=[0.2*pi0.8*pi];
>>ws=[0.4*pi0.7*pi];
>>Rp=1;Rs=30;
>>[N,wn]=buttord（wp/pi,ws/pi,Rp,Rs）;
>>[b,a]=butter（N,wn,'stop'）;
>>[C,B,A]=dir2cas（b,a）
C=
0.0394
B=
1.00000.35590.9994
1.00000.35471.0040
1.00000.35220.9954
1.00000.34991.0046
1.00000.34750.9960
1.00000.34631.0006
A=
1.00001.35680.7928
1.00001.03300.4633
1.00000.61800.1775
1.0000-0.24930.1113
1.0000-0.66170.3755
1.0000-0.97820.7446
>>[db,mag,pha,grd,w]=freqz_m（b,a）;
>>subplot（2,1,1）;plot（w/pi,mag）;
>>subplot（2,1,2）;plot（w/pi）;
数字低通---数字带阻：
function[bz,az]=zmapping（bZ,aZ,Nz,Dz）
bzord=（length（bZ）-1）*（length（Nz）-1）;
azord=（length（aZ）-1）*（length（Dz）-1）;
bz=zeros（1,bzord+1）;
fork=0:
bzord
pln=[1];
fori=0:
k-1
pln=conv（pln,Nz）;
end
pld=[1];
fori=0:
bzord-k-1
pld=conv（pld,Dz）;
end
bz=bz+bZ（k+1）*conv（pln,pld）;
end
fork=0:
azord
pln=[1];
fori=0:
k-1
pln=conv（pln,Nz）;
end
pld=[1];
fori=0:
azord-k-1
pld=conv（pld,Dz）;
end
az=az+aZ（k+1）*conv（pln,pld）;
end
all=az
（1）;az=az/az1;
bz=bz/az1;
线性相位FIR滤波器的幅度特性：
functionpzkplot（num,den）
holdon;
axis（'square'）;
x=-1:
0.01:
1;
y=（1-x.^2）.^0.5;
y1=-（1-x.^2）.^0.5;
plot（x,y,'b',x,y1,'b'）;
num1=length（num）;
den1=length（den）;
if（num1>1）
z=roots（num）;
else
z=0;
end
if（den1>1）
p=roots（den）;
else
p=0;
end
if（num>1&den1>1）
r_max_z=max（abs（real（z）））;
i_max_z=max（abs（imag（z）））;
a_max_z=max（r_max_z,i_max_z）;
r_max_p=max（abs（real（p）））;
i_max_p=max（abs（imag（p）））;
a_max_p=max（r_max_p,i_max_p）;
a_max=max（a_max_z,a_max_p）;
elseif（num1>1）
r_max_z=max（abs（real（z）））;
i_max_z=max（abs（imag（z）））;
a_max=max（r_max_z,i_max_z）;
else
r_max_p=max（abs（real（p）））;
i_max_p=max（abs（imag（p）））;
a_max=max（r_max_p,i_max_p）;
end
axis（[-a_maxa_max-a_maxa_max]）;
plot（[-a_maxa_max],[00],'b'）;
plot（[00],[-a_maxa_max],'b'）;
plot（[-a_maxa_max],[a_maxa_max],'b'）;
plot（[a_maxa_max],[-a_maxa_max],'b'）;
Lz=length（z）;
fori=1:
Lz;
plot（real（z（i））,imag（z（i））,'bo'）;
end
Lp=length（p）;
forj=1:
Lp
plot（real（p（j））,imag（p（j））,'bx'）;
end
title（'Thezeros-poleplot'）;
xlabel（'虚部'）;
ylabel（'实部'）;
function[Hr,w,a,L]=Hr_Type1（h）
M=length（h）;
L=（M-1）/2;
a=[h（L+1）2*h（L:
-1:
1）];
n=[0:
1:
L];
w=[0:
1:
500]'*pi/500;
Hr=cos（w*n）*a';
设计I型线性相位FIR滤波器：
>>h=[-41-1-2565-2-11-4];
>>M=length（h）;n=0:
M-1;
>>[Hr,w,a,L]=Hr_Type1（h）;
>>amax=max（a）+1;
>>amin=min（a）-1;
>>subplot（2,2,1）;stem（n,h）;
>>axis（[-12*L+1aminamax]）;text（2*L+1.5,amin,'n'）;
>>xlabel（'n'）;ylabel（'h（n）'）;title（'脉冲响应'）;
>>subplot（2,2,3）;stem（0:
L,a）;
>>axis（[-12*L+1aminamax]）;
>>xlabel（'n'）;ylabel（'a（n）'）;title（'a（n）系数'）;
>>subplot（2,2,2）;plot（w/pi,Hr）;
>>gridon;text（1.05,-20,'频率pi'）;
>>xlabel（'频率'）;ylabel（'Hr'）;title（'I型振幅响应'）;
>>subplot（2,2,4）;pzkplot（h,1）;
>>title（'零极点分布'）;
function[hr,w,b,L]=Hr_Type2（h）
M=length（h）;
L=M/2;
b=2*h（L:
-1:
1）;
n=[1:
1:
L];
n=n-0.5;
w=[0:
1:
500]'*pi/500;
hr=cos（w*n）*b';
II型线性相位FIR滤波器：
>>h=[-41-1-2565-2-11-4];
>>M=length（h）;n=0:
M-1;
>>[Hr,w,b,L]=Hr_Type2（h）;
Warning:
Integeroperandsarerequiredforcolonoperatorwhenusedasindex.
>InHr_Type2at2
>>bmax=max（b）+1;
bmin=min（b）-1;
>>subplot（2,2,1）;stem（n,h）;
axis（[-12*L+1bminbmax]）;text（2*L+1.5,bmin,'n'）;
xlabel（'n'）;ylabel（'h（n）'）;title（'脉冲响应'）;
>>subplot（2,2,3）;stem（1:
L,b）;
axis（[-12*L+1bminbmax]）;
xlabel（'n'）;ylabel（'b（n）'）;title（'b（n）系数'）;
>>subplot（2,2,2）;plot（w/pi,Hr）;
gridon;text（1.05,-20,'频率pi'）;
xlabel（'频率'）;ylabel（'Hr'）;title（'II型振幅响应'）;
>>subplot（2,2,4）;pzkplot（h,1）;
title（'零极点分布'）;
function[hr,w,c,L]=Hr_Type3（h）
M=length（h）;
L=（M-1）/2;
b=2*h（L+1:
-1:
1）;
n=[1:
1:
L];
w=[0:
1:
500]'*pi/500;
hr=cos（w*n）*c';
用MATLAB编程绘制各种窗函数的形状。
窗函数的长度为20：
>>Nwin=20;
n=0:
Nwin-1;
figure
（1）;
>>fori=1:
4
switchi
case1
w=boxcar（Nwin）;
stext='矩形窗';
case2
w=hanning（Nwin）;
stext='汉宁窗';
case3
w=hamming（Nwin）;
stext='海明窗';
case4
w=bartlett（Nwin）;
stext='布莱克曼窗';
end
posplot=['2,2',int2str（i）];
subplot（posplot）;
stem（n,w）;
holdon;
plot（n,w,'r'）;
xlabel（'n'）;ylabel（'w（n）'）;title（stext）;
holdoff;gridon;
End
用MATLAB编程，采用521个频率点绘制各种窗函数的幅频特性：
>>clf;
>>Nf=512;
>>Nwin=20;
>>figure
（1）;
>>fori=1:
4
switchi
case1
w=boxcar（Nwin）;
stext='矩形窗';
case2
w=hanning（Nwin）;
stext='汉宁窗';
case3
w=hamming（Nwin）;
stext='海明窗';
case4
w=bartlett（Nwin）;
stext='布莱克曼窗';
end
[y,f]=freqz（w,1,Nf）;
mag=abs（y）;
posplot=['2,2',int2str（i）];
sub
升级会员 