chebyshev.m

来自「是有关DSP中著名的契比雪夫的窗函数,还是不错的,欢迎使用」· M 代码 · 共 89 行

function chebyshev(wp,ws,ap,as,Fs,type)
e=sqrt(10.^(ap/10)-1);
wc=wp;
N=1;
if ws<=wp
   while 10*log10(e.^2*(cos(N*acos(ws/wp))).^2+1)<as
       N=N+1;
   end
else
   while 10*log10(e.^2*(cosh(N*acosh(ws/wp))).^2+1)<as
       N=N+1;
   end
end
f=0:0.001:0.5;
if ws<=wp
    A=1./(1+e.^2.*(cos(N*acos(2*pi*f/wc))).^2);
else 
    A=1./(e.^2.*(cosh(N*acosh(2*pi*f/wc))).^2+1);
end
subplot(2,1,1);
plot(f,sqrt(A));
title('chebyshev模拟低通滤波器');
c=1/e+sqrt(1/e.^2+1);
a=(sqrt(c)-1/sqrt(c))/2;
b=(sqrt(c)+1/sqrt(c))/2;
for k=1:N
    s(k)=(-a*sin((2*k-1)*pi/(2*N))+j*b*cos((2*k-1)*pi/(2*N)))*wc;
end
pp=poly(s);
k0=pp(length(pp));
fprintf('chebyshev模拟低通滤波器H(s)：\n');
fprintf('                 %.4g\n',k0);
fprintf('-----------------------------------------------------------------------------------------------------------------\n');
pstr1(pp,'s');
fprintf('\n');

f=0.01:0.0001:0.49;
H=1;
if type=='l'
    c=1;
    for k=1:N
        H=H./(((1-exp(-j*2*pi*f))./(1+exp(-j*2*pi*f)))-s(k));
        c=c/(1-s(k));
    end
    v=(1+s)./(1-s);
    pz=poly(v);
    fprintf('数字低通滤波器H(z)：\n');
    fprintf('             %.4g(z+1).^%.0f\n',c,N);
    fprintf('----------------------------------------------------------------------------------------------------------------\n');
    pstr1(pz,'z');
    fprintf('\n');
    
end
if type=='h'
    c=1;
    for k=1:N
        H=H./(((1+exp(-j*2*pi*f))./(1-exp(-j*2*pi*f)))-s(k));
        c=c/(s(k)-1);
    end
    v=(s+1)./(s-1);
    pz=poly(v);
    fprintf('数字高通滤波器H(z)：\n');
    fprintf('             %.4g(z-1).^%.0f\n',c,N);
    fprintf('----------------------------------------------------------------------------------------------------------------\n');
    pstr1(pz,'z');
    fprintf('\n');
    
end
if type=='bp'
    c=1;
    for k=1:N
        H=H./(a*(exp(j*2*2*pi*f)-2*b*exp(j*2*pi*f)+1)./(exp(j*2*2*pi*f)-1)-s(k));
        c=c/(a-s(k));
    end
    p=[a-s -2*a*b a+s];
    v=roots(p);
    pz=poly(v);
    fprintf('数字带通滤波器H(z)：\n');
    fprintf('             %.4g(z.^2-1).^%.0f\n',c,N);
    fprintf('----------------------------------------------------------------------------------------------------------------\n');
    pstr1(pz,'z');
    fprintf('\n');
end
subplot(2,1,2);
plot(f,20*log10(abs(H)));
title('数字滤波器');