FIR高通/低通/带通/带阻滤波器设计

wp=0.2*pi;wr=0.3*pi;Ap=1;Ar=15; T=1;Omegap=(2/T)*tan(wp/2);Omegar=(2/T)*tan(wr/2); %function [b,a]=afd_butt(Omegap,Omegar,Ap,Ar); if Omegap<=0 error('通带边缘必须大于0') end if (Omegar<=Omegap) error('阻带边缘必须大于通带边缘必须') end N=ceil((log10((10^(0.1*Ap)-1)/(10^(0.1*Ar)-1)))/(2*log10(Omegap/Omegar))); fprintf('\n***巴特沃斯模拟低通滤波器阶数=%2.0f\n',N) OmegaC=Omegap/((10^(0.1*Ap)-1)^(1/(2*N))) z=[]; p=exp(i*(pi*(1:2:N-1)/(2*N)+pi/2)); p=[p;conj(p)]; p=p(:); if rem(N,2)==1 p=[p;-1]; end k=real(-p); %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[A,B,C]=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); end elseif Na==1 A=[0 real(poly(Arow))]; 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)))]; end z=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); end elseif Nb==1 A=[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 %function [db,mag,pha,Omega]=freqs_m(b,a,Omega_max); Omega_max=0:50; Omega=Omega_max.*pi/50; H=freqs(b,a,Omega); mag=abs(H); dB=20*log10((mag+eps)/max(mag)); pha=angle(H); subplot(224);plot(Omega/pi,dB); title('模拟滤波器幅度响应|H|(j/Omegar)'); [b,a]=bilinear(b,a,T); %function[C,B,A]=dir2cas(b,a); b0=b(1);b=b/b0; a0=a(1);a=a/a0; C=b0/a0; M=length(a)-1;N=length(b)-1; p=cplxpair(roots(a));k=floor(Na/2); if N>M b=[b zeros(1,N-M)]; elseif M>N a=[a zeros(1,M-N)];N=M; else NM=0; end k=floor(N/2);B=zeros(k,3);A=zeros(k,3); if k*2==N; b=[b 0]; a=[a 0]; end broots=cplxpair(roots(b)); aroots=cplxpair(roots(a)); for i=1:2:2*k Brow=broots(i:1:i+1,:);Brow=real(poly(Brow)); B(fix((i+1)/2),:)=Brow; Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow)); A(fix((i+1)/2),:)=Arow; end %function [db,mag,pha,w]=freqz_m(b,a); [H,w]=freqz(b,a,1000,'whole'); H=(H(1:1:501))'; w=(w(1:1:501))'; mag=abs(H) db=20*log10((mag+eps)/max(mag)) pha=angle(H) subplot(221);plot(w/pi,mag);title('数字滤波器幅度响应|H|(ej/Omegar)'); subplot(222);plot(w/pi,db);title('数字滤波器幅度响应(db)'); subplot(223);plot(w/pi,pha/pi);title('数字滤波器相位响应'); delta_w=2*pi/1000 Ap=-(min(db(1:1:wp/delta_w+1))) Ar=-round(max(db(delta_w+1:1:501)))