myiir.m

matlb编写的IIR滤波器程序

%通带波纹AP(dB)
%阻带衰减AS(dB)
%通带截止频率为WP(rad/sec)
%阻带截止频率为WS(rad/sec)
%滤波器的采样周期T
%输入滤波器参数
AP=input('AP=');
if AP<0
    error
end
AS=input('AS=');
if AS<0
    error
end
WP=input('WP=');
WS=input('WS=');
if WP<0
    error
end
if WS<WP
    error
end
T=input('T=');
%%%%%%%%%%%%%%%%%%
%转换为模拟原型指标
Fs=1/T;                %得到采样频率
OmegaP=(2/T)*tan(WP/2);%带通截止频率预修正
OmegaS=(2/T)*tan(WS/2);%阻带截止频率预修正
%%%%%%%%%%%%%%%%%%
%计算N和OmegaC
N=ceil((log10((10^(AP/10)-1)/(10^(AS/10)-1)))/(2*log10(OmegaP/OmegaS)));
OmegaC=OmegaS/((10^(AS/10)-1)^(1/(2*N)));
%%%%%%%%%%%%%%%%%%
%计算Ha(s)的系数
[z0,p0,k0]=buttap(N);  %归一化巴特沃思原型设计函数
     p=p0*OmegaC;      %得非归一化零极点
     z=z0*OmegaC;
     k=k0*OmegaC^N;
     b0=k*real(poly(z));
     a0=real(poly(p));
 %%%%%%%%%%%%%%%%%
%双线形变换
[b,a]=bilinear(b0,a0,Fs);%b为分子系数，a为分母系数
%绘图
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';w=(w(1:1:501))'
mag=abs(H);
pha=angle(H);
figure;
subplot(2,1,1),plot(w/pi,mag);%绘制幅度响应曲线
subplot(2,1,2),plot(w/pi,pha/pi);%绘制相位响应曲线
figure;
freqz(b,a,512,Fs);%绘制频率响应曲线