用matlab設計的IIR濾波器源程序-IIR一階低通/高通

（6）IIR高通濾波（Butterworth）
% IIR Hightpass Use Butterworth
% copyright by Etual
clear;
fs=20;fpass=5;fstop=4;
Ap=0.5;As=10;
wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;
omegap=cot(wp/2);omegas=cot(ws/2);
ep=sqrt(10^(Ap/10)-1);
es=sqrt(10^(As/10)-1);
N=ceil(log(es/ep)/log(omegas/omegap));
omega0=omegap/ep^(1/N);
K=floor(N/2);
for i=1:K
theta(i)=pi*(N-1+2*i)/(2*N);
end
for i=1:K
G(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a1(i)=-2*(omega0^2-1)/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a2(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2);
end
if K(N/2)
G0=omega0/(omega0+1);a0=-(omega0-1)/(omega0+1);
end
w=(0+eps):pi/300:pi;
Hw2=1./(1+(cot(w/2)/omega0).^(2*N));
plot(w/pi,Hw2);
grid;

（7）IIR帶通濾波（Butterworth）
% IIR Bandpass Use Butterworth
% copyright by Etual
clear;
fs=20;fpa=2;fpb=4;fsa=1.5;fsb=4.5;
Ap=0.0877;As=16.9897;
wpa=2*pi*fpa/fs;wpb=2*pi*fpb/fs;wsa=2*pi*fsa/fs;wsb=2*pi*fsb/fs;
c=sin(wpa+wpb)/(sin(wpa)+sin(wpb));
omegap=abs((c-cos(wpb))/sin(wpb));
omegasa=(c-cos(wsa))/sin(wsa);omegasb=(c-cos(wsb))/sin(wsb);
omegas=min(abs(omegasa),abs(omegasb));
ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);
N=ceil(log(es/ep)/log(omegas/omegap));
omega0=omegap/ep^(1/N);
K=floor(N/2);
for i=1:K
theta(i)=pi*(N-1+2*i)/(2*N);
end
for i=1:K
G(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a1(i)=4*c*(omega0*cos(theta(i))-1)/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a2(i)=2*(2*c^2+1-omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a3(i)=-(4*c*(omega0*cos(theta(i))+1))/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a4(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2);
end
if K(N/2)
G0=omega0/(1+omega0);a0(1)=-2*c/(1+omega0);a0(2)=(1-omega0)/(1+omega0);
end
w=(0+eps):pi/300:pi;
Hw2=1./(1+((c-cos(w))./(omega0*sin(w))).^(2*N));
plot(w/pi,Hw2);
grid;

（8）IIR帶阻濾波（Butterworth）
% IIR Bandstop Use Butterworth
% copyright by Etual
clear;
fs=20;fpa=1.5;fpb=4.5;fsa=2;fsb=4;
Ap=0.5;As=10;
wpa=2*pi*fpa/fs;wpb=2*pi*fpb/fs;wsa=2*pi*fsa/fs;wsb=2*pi*fsb/fs;
c=sin(wpa+wpb)/(sin(wpa)+sin(wpb));
omegap=abs(sin(wpb)/(c-cos(wpb)));
omegasa=sin(wsa)/(cos(wsa)-c);omegasb=sin(wsb)/(cos(wsb)-c);
omegas=min(abs(omegasa),abs(omegasb));
ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);
N=ceil(log(es/ep)/log(omegas/omegap));
omega0=omegap/ep^(1/N);
K=floor(N/2);
theta=zeros(1,K);
for i=1:K
theta(i)=pi*(N-1+2*i)/(2*N);
end
G=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);
for i=1:K
G(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a1(i)=2*(omega0^2-1)/(1-2*omega0*cos(theta(i))+omega0^2);
end
for i=1:K
a2(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2);
end
if K(N/2)
G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);
end
w=(0+eps):pi/300:pi;
Hw2=1./(1+(sin(w)./(omega0*(c-cos(w)))).^(2*N));
plot(w/pi,Hw2);
grid;