巴特沃兹模拟滤波器的MATLAB实现.docx

巴特沃兹模拟滤波器的MATLAB实现设计一个模拟巴特沃特低通滤波器，它在30rad/s处具有1dB或更好的波动，在50rad/s处具有至少30dB的衰减。求出级联形式的系统函数，画出滤波器的幅度响应、对数幅度响应、相位响应和脉冲响应图。解：程序清单如下Wp=30;Ws=50;Rp=1;As=30; %技术指标Ripple=10(-Rp/20);Attn=10(-As/20);b,a=afd_butt(Wp,Ws,Rp,As) %巴特沃兹低通滤波器C,B,A=sdir2cas(b,a) %计算二阶节系数,级联型实现db,mag,pha,w=freqs_m(b,a,50); %计算幅频响应ha,x,t=impulse(b,a); %计算模拟滤波器的单位脉冲响应figure(1);clf;subplot(2,2,1);plot(w,mag);title(Magnitude Response);xlabel(Analog frequency in rad/s);ylabel(H);%axis(0,50,0,1.1)set(gca,XTickMode,manual,XTick,0,30,40,50);set(gca,YTickMode,manual,YTick,0,Attn,Ripple,1);gridsubplot(2,2,2);plot(w,db);title(Magnitude in dB);xlabel(Analog frequency in rad/s);ylabel(decibels);%axis(0,50,-40,5)set(gca,XTickMode,manual,XTick,0,30,40,50);set(gca,YTickMode,manual,YTick,-40,-As,-Rp,0);gridsubplot(2,2,3);plot(w,pha/pi);title(Phase Response);xlabel(Analog frequency in rad/s);ylabel(radians);%axis(0,50,-1.1,1.1)set(gca,XTickMode,manual,XTick,0,30,40,50);set(gca,YTickMode,manual,YTick,-1,-0.5,0,0.5,1);gridsubplot(2,2,4);plot(t,ha);title(Impulse Response);xlabel(time in seconds);ylabel(ha(t);axis(0,max(t)+0.05,min(ha),max(ha)+0.025);set(gca,XTickMode,manual,XTick,0,0.1,max(t);set(gca,YTickMode,manual,YTick,0,0.1,max(ha);grid%巴特沃兹模拟滤波器的设计子程序functionb,a=afd_butt(Wp,Ws,Rp,As);if Wp=0 error(Passband edge must be larger than 0)endif Ws=Wp error(Stopband edge must be larger than Passed edge)endif (Rp=0)|(As0) error(PB ripple and /0r SB attenuation must be larger than 0)end N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);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*OmegaCN;B=real(poly(z);b0=k;b=k*B;a=real(poly(p);%计算系统函数的幅度响应和相位响应子程序function db,mag,pha,w=freqs_m(b,a,wmax);w=0:1:500*wmax/500;H=freqs(b,a,w);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);%直接形式转换成级联形式子程序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);if K*2=Na A=zeros(K,3); for n=1:2:Na Arow=p(n:1:n+1,:);Arow=poly(Arow); A(fix(n+1)/2),:)=real(Arow); endelseif Na=1 A=0 real(poly(p);else A=zeros(K+1,3); for n=1:2:2*K Arow=p(n:1:n+1,:);Arow=poly(Arow); A(fix(n+1)/2),:)=real(Arow); end A(K+1,:)=0 real(poly(p(Na);endz=cplxpair(roots(b);K=floor(Nb/2);if Nb=0 B=0 0 poly(z);elseif K*2=Nb B=zeros(K,3); for n=1:2:Nb Brow=z(n:1:n+1,:);Brow=poly(Brow); B(fix(n+1)/2),:)=real(Brow); endelseif Nb=1 B=0 real(poly(z);else B=zeros(K+1,3); for n=1:2:2*K Brow=z(n:1:n+1,:);Brow=poly(Brow); B(fix(n+1)/2),:)=real(Brow); end B(K+1,:)=0 real(poly(z(Nb);end 运行结果如下（如图6.10所示）。图6.10 巴特沃兹模拟低通滤波器例b = 3.8682e+013a = 1.0e+013 * Columns 1 through 10 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0036 0.0613 0.6888 3