design_chebyshev.m

经典《信号与系统》教程的matlab例程,对深入理解信号与系统相关概念有很大帮助

% Name: desi_chebyshev
% This program is used to design the type 1 Chebyshev lowpass analog filter
% The specifications are:
% alphap: The maximal passband attenuation in dB
% alphas: The minimum stopband attenuation in dB
% wp: The passband cutoff frequency in radians
% ws: The stopband cutoff frequency in radian
% N:  The order of the lowpass filter
clear
format long
alphap=0.1;%input('Type in the data ap=:');
alphas=60;%input('Type in the data as=:');
wp=2*pi*3000;%input('Type in the passband cutoff frequency wp=:');
ws=2*pi*12000;%input('Type in the stopband cutoff frequency ws=:');
lamdp=wp/wp;lamds=ws/wp;       % Compute the normalized frequencies

k1=sqrt((10.^(0.1*alphas)-1)/(10.^(0.1*alphap)-1)); % Compute the order N
N=ceil(acosh(k1)/acosh(lamds))
e=sqrt(10^(0.1*alphap)-1);
a=e*2^(N-1);
ksai=asinh(1/e)/N;

for k=1:N
    realp(k)=-sinh(ksai)*sin(pi*(2*k-1)/(2*N));  % Compute the normalized poles
    imagp(k)=cosh(ksai)*cos(pi*(2*k-1)/(2*N));
end
p=realp+j*imagp
s=wp*realp+j*wp*imagp;
for i=1:fix(N/2)              % Form the binomial with real coefficients
    coef2(i)=-2*(real(p(i)));
    coef3(i)=real(p(i)*p(N-i+1));
end
if N/2~=fix(N/2)
    coef2(fix(N/2)+1)=realp(fix(N/2)+1);
    coef3(fix(N/2)+1)=0;
    b=ones(fix(N/2)+1,1);
    coef2=coef2';
    coef3=coef3';
    b=[b,coef2,coef3]
    sos=zeros(fix(N/2)+1,2);
    I=ones(fix(N/2)+1,1);
    sos=[sos,I,b];
else
    b=ones((N/2),1);
    coef2=coef2';
    coef3=coef3';
    b=[b,coef2,coef3]
    sos=zeros((N/2),2);
    I=ones((N/2),1);
    sos=[sos,I,b];
end

K=1/a;
[num,den]=sos2tf(sos,K)
w=0:0.01:4;
subplot(221)
[h,w]=freqs(num,den,w);
plot(w,abs(h))
title('The normalized magnitude response')
xlabel('Normalized frequency')

for i=1:N+1
    den(i)=(wp^(i-1))*den(i);
end
num=(wp^(N))*num;
w=0:0.5:4*wp;
[h,w]=freqs(num,den,w);
subplot(222)
plot(w/(2*pi),abs(h));%axis([0,4,0,1.2])
title('The magnitude response')
xlabel('Frequency')

%coef2=coef2';
%coef3=coef3';
%b=[b,coef2,coef3]
%sos=zeros(fix(N/2)+1,2);
%I=ones(fix(N/2)+1,1);
%sos=[sos,I,b]
    
%((realp(1))^2)/((sinh(ksai))^2)+((imagp(1))^2)/((cosh(ksai))^2)