iir.m

用Matlab实现的IIR滤波器

%用Butterworth滤波器逼近理想模拟低通滤波器
%用双线性变换法对模拟滤波器进行数字化
%Ap为通带最大衰减，单位为dB
%As为阻带最小衰减，单位为dB
%fp为截止频率，单位为Hz
%fs为阻带起始频率，单位为Hz
%Fs为抽样频率，单位为Hz
clear
fp=input('fp=');
fs=input('fs=');
Ap=input('Ap=');
As=input('As=');
Fs=input('Fs=');

%以下是用Butterworth滤波器逼近理想模拟低通滤波器
%先把给定的数字滤波器指标转换为模拟低通滤波器的技术指标
wp=2*pi*fp/Fs;                          %求通带截止频率wp,2*pi对应Fs
ws=2*pi*fs/Fs;                          %求阻带下限截止频率ws,
Wp=2*Fs*tan(wp/2);                      %将数字滤波器的技术要求转化为模拟滤波器的要求
                                        %Wp=(2/Ts)*tan(wp/2)  
Ws=2*Fs*tan(ws/2);                      %Ws=(2/Ts)*tan(ws/2)
%Butterworth滤波器的设计
N=ceil(0.5*(log10(10.^(As/10)-1)-log10(10.^(Ap/10)-1))/(log10(Ws)-log10(Wp)));
                                        %N为滤波器阶次
Wc=Ws/((10.^(As/10)-1).^(1/(2*N)));     %截止频率
sk=[];                                  %求极点,把左半面的值赋给sk
for k=1:2*N
    sk(k)=Wc*exp(j*pi*(2*k+N-1)/(2*N));
end

%用双线性变换将模拟滤波器数字化
%s=2*Fs*(z-1)/(1+z)
%可由Gs得Hz
c=2*Fs;
b=[1];
a=[1]; 
for j=1:N                                       % Ha(s)转化为H(Z)时
    b=conv(b,[1,1]);                            % 分子为 (1+1/z).^N
    a=conv(a,[-c-sk(j),c-sk(j)]);               % 分母为 (1+1/z).^N* ∏(s-sk) =∏ [c(1-1/z)-(1+1/z)*sk]=∏ [c-sk-(c+sk)1/z]
end
b=Wc.^N*b;
a=a(end:-1:1);                                  %使系数按(1/z)的升序排列
b=b(end:-1:1);
b=b/a(1);                                       %使a(1)为1
a=a/a(1);

%输出幅频响应
[H,w]=freqz(b,a);
figure(1)
f=w/max(w)*Fs/2;
H=20*log10(abs(H)/abs(H(1)));
plot(f,H);
grid;
xlabel('w  /rad');
ylabel('幅度 /db');
axis([0,Fs/2,-30,0]);