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

（9）IIR低通濾波（chebyshev 1）
% IIR Lowpass Use Chebyshev Type 1
% copyright by Etual
clear;
fs=20;fpass=4;fstop=5;
Ap=0.5;As=10;
wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;
omegap=tan(wp/2);omegas=tan(ws/2);
ep=sqrt(10^(Ap/10)-1);
es=sqrt(10^(As/10)-1);
e=es/ep;w=omegas/omegap;
N=ceil(log(e+sqrt(e^2-1))/log(w+sqrt(w^2-1)));
a=log(1/ep+sqrt(1/ep^2+1))/N;
omega0=omegap*sinh(a);
K=floor(N/2);
theta=zeros(1,K);omega=zeros(1,K);
for i=1:K
theta(i)=pi*(N-1+2*i)/(2*N);
end
for i=1:K
omega(i)=omegap*sin(theta(i));
end
G=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);
for i=1:K
G(i)=(omega0^2+omega(i)^2)/(1-2*omega0*cos(theta(i))+omega0^2+omega(i)^2);
end
for i=1:K
a1(i)=2*(omega0^2+omega(i)^2-1)/(1-2*omega0*cos(theta(i))+omega0^2+omega(i)^2);
end
for i=1:K
a2(i)=(1+2*omega0*cos(theta(i))+omega0^2+omega(i)^2)/(1-2*omega0*cos(theta(i))+omega0^2+omega(i)^2);
end
if K(N/2)
G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);
else
H0=sqrt(1/(1+ep^2));
end
f=0:1/300:10;
Hf2=1./(1+ep^2*(cheby(N,tan(pi*f/fs)/omegap)).^2);
plot(f,abs(Hf2));
grid;
（9）IIR低通濾波（chebyshev 1）
% IIR Lowpass Use Chebyshev Type 2
% copyright by Etual
clear;
fs=20;fpass=4;fstop=5;
Ap=0.5;As=10;
wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;
omegap=tan(wp/2);omegas=tan(ws/2);
ep=sqrt(10^(Ap/10)-1);
es=sqrt(10^(As/10)-1);
e=es/ep;w=omegas/omegap;
N=ceil(log(e+sqrt(e^2-1))/log(w+sqrt(w^2-1)));
a=log(es+sqrt(es^2+1))/N;
omega0=omegas/sinh(a);
K=floor(N/2);
for i=1:K
theta(i)=pi*(N-1+2*i)/(2*N);
end
for i=1:K
omega(i)=omegas/sin(theta(i));
end
for i=1:K
G(i)=(1+omega(i)^-2)/(1-2*omega0^-1*cos(theta(i))+omega0^-2+omega(i)^-2);
end
for i=1:K
a1(i)=2*(1-omega0^-2+omega(i)^-2)/(1-2*omega0^-1*cos(theta(i))+omega0^-2+omega(i)^-2);
end
for i=1:K
a2(i)=(1+2*omega0^-1*cos(theta(i))+omega0^-2+omega(i)^-2)/(1-2*omega0^-1*cos(theta(i))+omega0^-2+omega(i)^-2);
end
for i=1:K
b1(i)=2*(1-omega(i))/(1+omega(i));
end
if K(N/2)
G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);
else
H0=sqrt(1/(1+ep^2));
end
f=(0+eps):1/100:10;
Hf2=(cheby(N,omegas./tan(pi*f/fs))).^2./((cheby(N,omegas./tan(pi*f/fs))).^2+es^2);
plot(f,abs(Hf2));
grid;

（10）chebyshev 中用到的函數 cheby.m
function CN=cheby(N,x)
if x=1
CN=cos(N*acos(x));
else
CN=cosh(N*log(x+sqrt(x.^2-1)));
end