bp_butt_design_bi.m

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

% Name: bp_butt_design_bi
% This program is used to design the IIR digital highpass 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
%======================================================================
% Input the specifications of the highpass digital filter
alphap=3;%input('Type in the data ap=:');
alphas=18;%input('Type in the data as=:');
wu=0.4*pi;%input('Type in the passband cutoff frequency wp=:');
wl=0.3*pi;%input('Type in the stopband cutoff frequency ws=:');
ws2=0.5*pi;
ws1=0.2*pi;
T=1;
%========================================================================
% specifications of the bandpass digital filter to specifications
% of bandpass analog filter
omegau=2*tan(wu/2)/T;
omegal=2*tan(wl/2)/T;
omegas2=2*tan(ws2/2)/T;
omegas1=2*tan(ws1/2)/T;
omega0=sqrt(omegau*omegal);
B=omegau-omegal;
nu=omegau/B;
ns2=omegas2/B;
n0=omega0/B;
%==========================================================================
%  Analog highpass filter to analog lowpass filter
lamds=(ns2^2-n0^2)/ns2;
lamdp=1;

%=====================================================
% Design the analog lowpass filter
ksp=sqrt((10^(0.1*alphap)-1)/(10^(0.1*alphas)-1));
lamdsp=lamds/lamdp;
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 lowpass filter are %2.0f\n');
b,a
%==========================================================================
% Analog lowpass filter to analog bandpass filter
[bb,ab]=lp2bp(a,b,omega0,B);
bb,ab
%==========================================================================
[Ha,wa]=freqs(bb,ab);Ha0=max(abs(Ha));
%=================================================================

fprintf('\n***The coefficient vectors of the digital filter are %2.0f\n');
[bz,az]=bilinear(bb,ab,1/T);
bz,az
[H,w]=freqz(bz,az);H0=max(abs(H));
subplot(211)
plot(w/pi,abs(H)/H0);
title('Amplitude response'),ylabel('|Ha|'),axis([0,1,0,1.1]),
grid on
subplot(212)
plot(w/pi,20*log(abs(H)/H0));
ylabel('dB'),axis([0,1,-40,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