experiment_3.m

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

% Name: experiment_3.m
% 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,close all
format short
alphap=0.5;%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=2*tan(wp/2)/T;
omegas=2*tan(ws/2)/T;
%==========================================================================
ksp=sqrt((10^(0.1*alphap)-1)/(10^(0.1*alphas)-1));
lamdsp=omegas/omegap;
N=ceil(-log(ksp)/log(lamdsp));
omegac=omegas*(10^(0.1*alphas)-1)^(-1/(2*N));
fprintf('\n***The order N of Butterworth is %2.0f\n',N);
%===========================================================================
wn=wp/pi;
[z,p,k]=buttap(N);
[b,a]=zp2tf(z,p,k);
fprintf('\n***The coefficient vectors of Butterworth are %2.0f\n');
b
a
%==========================================================================
ba=b;aa=a;
for m=1:N;
    aa(m+1)=a(m+1)*omegac^(m);
    ba(m+1)=b(m+1)*omegac^(m);
end
fprintf('\n***The non-normlized coefficient vectors of Butterworth are %2.0f\n');
ba,aa
%==========================================================================
fprintf('\n***The coefficient vector of IIR digital filter in direct form are %2.0f\n');
[Bz,Az]=bilinear(ba,aa,1/T)
fprintf('\n***The coefficient vector of IIR digital filter in parallel form are %2.0f\n');
[cz,bz,az] = dir2par(Bz,Az)
%==========================================================================
[H,w]=freqz(Bz,Az,400);
%==========================================================================
x=[-4 -2 0 -4 -6 -4 -2 -4 -6 -6 -4 -4 -6 -6 -2 6 12 8 0 -16 -38 -60 -84 -90 -66 -32 -4 -2 -4 8 12 12 10 6 6 6 4 0 0 0 0 0 -2 -4 0 0 0 -2 -2 0 0 -2 -2 -2 -2 0];
y=filter(Bz,Az,x);
L=length(x);
n=0:L-1;
X=fft(x,L);Y=fft(y,L);
subplot(221)
stem(n,x,'.');title('The input sequence x(n)')
axis([0,L-1,min(x),max(x)]),
subplot(223)
stem(n,y,'.');title('The output sequence y(n)')
axis([0,L-1,min(y),max(y)]),xlabel('Index n')
subplot(222)
stem(n*(2/L),abs(X),'.');title('The DFT of the sequence x(n)')
axis([0,1,0,max(abs(X))]),
subplot(224)
stem(n*(2/L),abs(Y),'.');title('The DFT of the sequence y(n)')
axis([0,1,0,max(abs(Y))]),xlabel('Frequency *pi')

figure(2)
subplot(211)
plot(w/pi,abs(H));
title('Amplitude response of the IIR digital filter'),ylabel('|H|'),axis([0,1,0,1.1]),
grid on
subplot(212)
plot(w/pi,20*log(abs(H)));
ylabel('dB'),axis([0,1,-80,3]),
xlabel('Frequency (pi)'),grid on
%==========================================================================