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