数字滤波器的设计及其MATLAB实现文档格式.docx
- 文档编号:21317807
- 上传时间:2023-01-29
- 格式:DOCX
- 页数:29
- 大小:252.07KB
数字滤波器的设计及其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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字滤波器 设计 及其 MATLAB 实现
![提示](https://static.bdocx.com/images/bang_tan.gif)