数字滤波器的设计及其MATLAB实现文档格式.docx
《数字滤波器的设计及其MATLAB实现文档格式.docx》由会员分享，可在线阅读，更多相关《数字滤波器的设计及其MATLAB实现文档格式.docx（29页珍藏版）》请在冰豆网上搜索。
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）;
2*K
Arow=poly（Arow）;
A（K+1,:
）=[0real（poly（p（Na）））];
z=cplxpair（roots（b））;
K=floor（Nb/2）;
ifNb==0
B=[00poly（z）];
elseifK*2==Nb
B=zeros（K,3）;
Nb
Brow=z（n:
Brow=poly（Brow）;
B（（fix（n+1）/2）,:
）=real（Brow）;
elseifNb==1
B=[0real（poly（z））];
B=zeros（K+1,3）;
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:
mag=abs（H）;
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:
N-2
br=r（i:
i+1,:
ar=p（i:
[br,ar]=residuez（br,ar,[]）;
B（（fix（i+1）/2）,:
）real（br'
A（（fix（i+1）/2）,:
）real（ar'
[br,ar]=residuez（r（N-1）,p（N-1）,[]）;
B（K,:
）=[real（br'
）0];
A（K,:
）=[real（ar'
N-1
）real（br）;
）real（ar）;
比较两个含同样标量元素但（可能）有不同下标的复数对及其相位留数向量子程序：
functionI=cplxcomp（p1,p2）
I=[];
fori=1:
length（p2）
forj=1:
length（p1）
if（abs（p1（j）-p2（i））<
0.0001）
I=[I,j];
I=I'
双线性变换巴特沃斯低通滤波器设计：
巴特沃思模拟滤波器的设计子程序：
function[b,a]=afd_butt（wp,ws,Rp,rs）
ifwp<
ifws<
=wp
阻带边缘必须大于通带边缘'
if（Rp<
=0）|（Rs<
N=ceil（（log10（（10^（Rp/10）-1）/（10^（Rs/10）-1）））/（2*log10（wp/ws）））;
\n***ButterworthFilterOrder=%2.0f\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）;
k=k*OmegaC^N;
直接型到级联型形式的转换：
function[b0,B,A]=dir2cas（b,a）
b0=b0/a0;
ifN>
M
b=[b,zeros（1,N-M）];
elseifN<
a=[a,zeros（1,M-N）];
NM=0;
k=floor（N/2）;
B=zeros（k,3）;
A=zeros（k,3）;
ifk*2==N
b=[b,0];
a=[a,0];
broots=cplxpair（roots（b））;
aroots=cplxpair（roots（a））;
2*k
br=broots（i:
br=real（polt（br））;
）=br;
ar=aroots（i:
ar=real（polt（ar））;
）=ar;
function[db,mag,pha,grd,w]=freqz_m（b,a）
[h,w]=freqz（b,a,1000,'
h=（h（1:
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
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];
数字通带边缘频率
Rs=40;
Ripple=10^（-Rp/20）;
通带波动
Attn=10^（-Rs/20）;
阻带衰减
[N,wn]=ellipord（wp/pi,ws/pi,Rp,Rs）计算椭圆滤波器参数
0.40000.6000
[b,a]=ellip（N,Rp,Rs,wn）;
数字椭圆滤波器的设计
[b0,B,A]=dir2cas（b,a）级联形式实现
b0=
0.0197
1.00001.50661.0000
1.00000.92681.0000
1.0000-0.92681.0000
1.0000-1.50661.0000
1.00000.59630.9399
1.00000.27740.7929
1.0000-0.27740.7929
1.0000-0.59630.9399
figure
（1）;
subplot（2,2,1）;
gridon;
subplot（2,2,3）;
gridon;
subplot（2,2,2）;
plot（w/pi,pha/pi）;
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];
Rs=30;
[N,wn]=buttord（wp/pi,ws/pi,Rp,Rs）;
stop'
0.0394
1.00000.35590.9994
1.00000.35471.0040
1.00000.35220.9954
1.00000.34991.0046
1.00000.34750.9960
1.00000.34631.0006
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
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）;
pld=[1];
bzord-k-1
pld=conv（pld,Dz）;
bz=bz+bZ（k+1）*conv（pln,pld）;
azord
azord-k-1
az=az+aZ（k+1）*conv（pln,pld）;
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,'
num1=length（num）;
den1=length（den）;
if（num1>
1）
z=roots（num）;
z=0;
if（den1>
p=roots（den）;
p=0;
if（num>
1&
den1>
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>
a_max=max（r_max_z,i_max_z）;
a_max=max（r_max_p,i_max_p）;
axis（[-a_maxa_max-a_maxa_max]）;
plot（[-a_maxa_max],[00],'
plot（[00],[-a_maxa_max],'
plot（[-a_maxa_max],[a_maxa_max],'
plot（[a_maxa_max],[-a_maxa_max],'
Lz=length（z）;
Lz;
plot（real（z（i））,imag（z（i））,'
bo'
Lp=length（p）;
forj=1:
Lp
plot（real（p（j））,imag（p（j））,'
bx'
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:
L];
w=[0:
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;
stem（n,h）;
axis（[-12*L+1aminamax]）;
text（2*L+1.5,amin,'
n'
xlabel（'
h（n）'
脉冲响应'
stem（0:
L,a）;
a（n）'
a（n）系数'
plot（w/pi,Hr）;
text（1.05,-20,'
频率pi'
频率'
Hr'
I型振幅响应'
pzkplot（h,1）;
title（'
零极点分布'
function[hr,w,b,L]=Hr_Type2（h）
L=M/2;
b=2*h（L:
1）;
n=[1:
n=n-0.5;
hr=cos（w*n）*b'
II型线性相位FIR滤波器：
[Hr,w,b,L]=Hr_Type2（h）;
Warning:
Integeroperandsarerequiredforcolonoperatorwhenusedasindex.
InHr_Type2at2
bmax=max（b）+1;
bmin=min（b）-1;
axis（[-12*L+1bminbmax]）;
text（2*L+1.5,bmin,'
stem（1:
L,b）;
b（n）'
b（n）系数'
II型振幅响应'
function[hr,w,c,L]=Hr_Type3（h）
b=2*h（L+1:
hr=cos（w*n）*c'
用MATLAB编程绘制各种窗函数的形状。
窗函数的长度为20：
Nwin=20;
Nwin-1;
figure
（1）;
4
switchi
case1
w=boxcar（Nwin）;
stext='
矩形窗'
case2
w=hanning（Nwin）;
汉宁窗'
case3
w=hamming（Nwin）;
海明窗'
case4
w=bartlett（Nwin）;
布莱克曼窗'
posplot=['
2,2'
int2str（i）];
subplot（posplot）;
stem（n,w）;
holdon;
plot（n,w,'
r'
w（n）'
title（stext）;
holdoff;
用MATLAB编程，采用521个频率点绘制各种窗函数的幅频特性：
clf;
Nf=512;
[y,f]=freqz（w,1,Nf）;
mag=abs（y）;
sub
升级会员 