hi_butt_design_bi.m

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

% Name: hi_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=15;%input('Type in the data as=:');
wp=0.8*pi;%input('Type in the passband cutoff frequency wp=:');
ws=0.44*pi;%input('Type in the stopband cutoff frequency ws=:');
T=1;
%========================================================================
% specifications of the highpass digital filter to specifications
% of highpass analog filter
omegap=2*tan(wp/2)/T;
omegas=2*tan(ws/2)/T;
%==========================================================================
%  Analog highpass filter to analog lowpass filter
omegap=1/omegap;
omegas=1/omegas;
%=====================================================
% Design the analog lowpass filter
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=omegap;
for i=1:N
    b(i+1)=b(i+1)*omegac^(i);
end
a=a*omegac^N;a,b
%==========================================================================
% Analog lowpass filter to analog highpass filter
bh=fliplr(b);
ah=zeros(1,N+1);
ah(1)=a;

%==========================================================================
[Ha,wa]=freqs(ah,bh);Ha0=max(abs(Ha));
%=================================================================
subplot(221)
plot(wa/pi,abs(Ha)/Ha0);
title('Amplitude response'),ylabel('|Ha|'),axis([0,1,0,1.1]),
grid on
subplot(223)
plot(wa/pi,20*log(abs(Ha)/Ha0));
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]=bilinear(ah,bh,1/T);
bz,az
[H,w]=freqz(bz,az);H0=max(abs(H));
subplot(222)
plot(w/pi,abs(H)/H0);
title('Amplitude response'),ylabel('|Ha|'),axis([0,1,0,1.1]),
grid on
subplot(224)
plot(w/pi,20*log(abs(H)/H0));
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