zhang_butt.m

给出了MATLAB语言实现IIR滤波器的实现

function [b,a] = zhang_butt(WP,WS,AP,AS,FS)  % 给出的是数字滤波器的指标，WP,WS是数字角频率，Ws,Wp是模拟角频率
%butterworth低通滤波器设计
if WP<=0
   error('Passband edge must be larger than 0')
end
if WS<=WP
   error('Stopband edge must be larger than Passband edge')
end
if (AP<=0)|(AS<0)
   error('PB ripple and/or SB attenuation must be large than 0')
end
Ws=2*FS*tan(WS/2);
Wp=2*FS*tan(WP/2);
es=10^(AS/10)-1;
ep=10^(AP/10)-1;
N=ceil(0.5*log10(es/ep)/log10(Ws/Wp))

%求极点Sk确定模拟滤波器转移函数
k=1;
S=zeros(1,N);   % ZEROS(M,N) or ZEROS([M,N]) is an M-by-N matrix of zeros.
for k=1:N
S(k)=Wp*exp(j*pi*(0.5+0.5*(2*k-1)/N));
end
aa=poly(S)   %aa为所设计的模拟滤波器的分母多项式的系数  POLY Convert roots to polynomial.
ba=[Wp^N];      %ba为所设计的模拟滤波器的分子多项式的系数

%以下为双线性变换求出数字滤波器分子分母系数向量的过程 《数字信号处理教程－MATLAB释义与实现》P285
lba=length(ba);laa=length(aa);ld=laa-lba;
if ld>=0
    ba=[zeros(1,ld),ba];
else aa=[zeros(1,-ld),aa];
end
Nz=FS*[2,-2];Dz=[1,1];   %双线性变换分子分母系数向量,注意这里究竟是[2,-2]还是[-2,2]决定了最后结果的形式，若是[2,-2]则最后结果也以降幂排序，最高项系数为1
bd=0;ad=0;                 %bd,ad系数向量初始化
for k=0:N
    plb=[1];pla=[1];       %双线性变换分子分母系数乘积项初始化
    for l=0:k-1
        plb=conv(plb,Dz);  %求双线性变换分母系数k次幂Dz^k
    end
    for l=0:N-k-1
        pla=conv(pla,Nz);  %求双线性变换分子系数(N-k)次幂Nz^(N-k)
    end
    bd=bd+ba(k+1)*conv(pla,plb);   %分子系数多项式向量求和
    ad=ad+aa(k+1)*conv(pla,plb);   %分母系数多项式向量求和
end
adl=ad(1);a=ad/adl;b=bd/adl;  %用分母系数多项式的首项使分子分母系数归一化
%双线性变换结束