butt_design_im1.m

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

% Name: iir_butt_design
% This program is used to design the IIR digital lowpass filter
% using the bilinear transformation
% The analog lowpass filter is the Butterworth
% The specifications of the digital filter are:
% alphap: The maximal passband attenuation in dB
% alphas: The minimum stopband attenuation in dB
% wp: The passband cutoff frequency 
% ws: The stopband cutoff frequency
% N:  The order of the lowpass filter
clear
format short
alphap=1;%input('Type in the data ap=:');
alphas=15;%input('Type in the data as=:');
wp=0.2*pi;%input('Type in the passband cutoff frequency wp=:');
ws=0.3*pi;%input('Type in the stopband cutoff frequency ws=:');
T=1;
omegap=wp/T;
omegas=ws/T;
%==========================================================================
ksp=sqrt((10^(0.1*alphap)-1)/(10^(0.1*alphas)-1));
lamdsp=omegas/omegap;
N=ceil(-log(ksp)/log(lamdsp));
[b]=butterworth(N);a=1;
omegac=omegap*(10^(0.1*alphap)-1)^(-1/(2*N));
fprintf('\n***The order N of Butterworth is %2.0f\n',N);
%===========================================================================
fprintf('\n***The coefficient vectors of the analog filter are %2.0f\n');
b,a
fprintf('\n***The the 3-dB frequency the analog filter is %2.0f\n');
omegac
for i=1:N
    b(i+1)=b(i+1)*omegac^(i);
end
a=a*omegac^N;

%==========================================================================
[Ha,wa]=freqs(a,b);
%=================================================================
subplot(221)
plot(wa/pi,abs(Ha));
title('Amplitude response'),ylabel('|Ha|'),axis([0,1,0,1.1]),
grid on
subplot(223)
plot(wa/pi,20*log(abs(Ha)));
ylabel('dB'),axis([0,1,-80,0]),
xlabel('Frequency (pi)'),grid on
%==========================================================================
fprintf('\n***The coefficient vectors of the digital filter are %2.0f\n');
[bz,az]=impinvar(a,b,1/T);
bz,az
[H,w]=freqz(bz,az);
subplot(222)
plot(w/pi,abs(H));
title('Amplitude response'),ylabel('|Ha|'),axis([0,1,0,1.1]),
grid on
subplot(224)
plot(w/pi,20*log(abs(H)));
ylabel('dB'),axis([0,1,-80,0]),
xlabel('Frequency (pi)'),grid on
fprintf('\n***The parallel form of the digital filter %2.0f\n');
[c,b,a]=dir2par(bz,az);
c,b,a