% Program 8.1, p511. m-file for the analog filter design
% (Example 8.19) % program name: prog81.m
%
FN=1000/2; fc=300; % cut offf frequency
N=5; % filter order
[z, p, k]=buttap(N); % create an analog filter
w=linspace(0, FN/fc, 1000); % plot the response of filter
h=freqs(k*poly(z), poly(p), w);
f=fc*w;
plot(f, 20*log10(abs(h))), grid
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
=============================================
% % Program 8.2, p511. m-file for the design of the
% impulse invariant filter (Example 8.19) % program name: prog82.m
%
Fs=1000; % sampling frequency
fc=300; % cutoff frequency
WC=2*pi*fc; % cutoff frequency in radian
N=5; % filter order
[b,a]=butter(N,WC,'s'); % create an analog filter
[z, p, k]=butter(N, WC, 's');
[bz, az]=impinvar(b,a,Fs); % determine coeffs of IIR filter
[h, f]=freqz(bz, az, 512,Fs);
plot(f, 20*log10(abs(h))), grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
==============================================
% % Program 8.3, p512. m-file for the design of the BZT filter
% (Example 8.19) % program name: prog83.m %
Fs=1000; % sampling frequency
FN=Fs/2;
fc=300; % cutoff frequency
N=5;
[z, p, k]=butter(N, fc/FN);
[h, f]=freqz(k*poly(z), poly(p), 512, Fs);
plot(f, 20*log10(abs(h))), grid
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
=====================================================
% % m-file for Example 8B.1 (Program 8B.1, p564).
% Program name: prog8b1.m %
N=2; % Filter order
Fs=1280; % Sampling frequency
fc=150; % Cutoff frequency
WC=2*pi*fc; % Cutoff frequency in radian
%
%Create an analogue filter
%
[b,a]=butter(N,WC,'s');
[z,p,k]=butter(N, WC, 's');
%
%Convert analogue filter into Discrete IIR filter
%
[bz, az]=impinvar(b, a, Fs); %Determine coeffs of IIR filter
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(bz, az, 512, Fs);
plot(f, 20*log10(abs(H)))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response (dB)')
subplot(2,1,2) % Plot pole-zero diagram
zplane(bz, az)
zz=roots(bz); % Determine poles and zeros
pz=roots(az);
============================================================
% % m-file for Example 8B.2 (Program 8B.2, p566).
%Program name: prog8b2.m %
N=2; % Filter order
Fs=1280; % Sampling frequency
FN=Fs/2;
fc=150;% Cutoff frequency
Fc=fc/FN; % Normalized Cutoff frequency
[b,a]=butter(N,Fc); % Create and digitize analogue filter.
[z,p,k]=butter(N, Fc);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
b
a
==========================================================
% % m-file for Example 8B.3 (Program 8B.3, p567)
% Program name: prog8b3.m
%
Fs=2000; % Sampling frequency
FN=Fs/2;
fc1=200/FN;
fc2=300/FN;
[b,a]=butter(4,[fc1, fc2]); % Create and digitize analogue filter.
[z,p,k]=butter(4, [fc1, fc2]);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
=========================================================
% % m-file for Example 8B.4 (Program 8B.4, p569)
% Program name: prog8b4.m %
Fs=8000; % Sampling frequency
Ap=3;
As=20;
wp=500/4000;
ws=2000/4000;
[N, wc]=buttord(wp, ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,500/4000); % Digitise filter
[b, a]=butter(N, 500/4000);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
=======================================================
% % m-file for Example 8B.5 (Program 8B.5, p571)
% Program name: prog8b5.m
%
Fs=8000; % Sampling frequency
Ap=3;
As=20;
wp=2000/4000;
ws=500/4000;
[N, wc]=buttord(wp, ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,2000/4000, 'high'); % Digitise filter
[b, a]=butter(N, 2000/4000, 'high');
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
====================================================
% % m-file for Example 8B.6 (Program 8B.6, p572)
%Program name: prog8b6.m
%Band pass filter
%
Fs=1000; % Sampling frequency
Ap=3;
As=20;
Wp=[200/500, 300/500]; % Band edge frequencies
Ws=[50/500, 450/500];
[N, Wc]=buttord(Wp, Ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,Wp); % Digitise filter
[b, a]=butter(N, Wp);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
=======================================================
% % m-file for Example 8B.7 (Program 8B.7, p574)
%Program name: prog8b7.m % Band stop filter
%
Fs=1000; %Sampling frequency
Ap=3;
As=20;
Wp=[50/500, 450/500]; % Band edge frequencies
Ws=[200/500, 300/500];
[N, Wc]=buttord(Wp, Ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,Ws, 'stop'); % Digitise filter
[b, a]=butter(N, Ws, 'stop');
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
============================================================
% % m-file for the design Example 8B.8 (Program 8B.8, p576)
%Program Name: prog8b8.m %
Ap=0.25;
As=45;
Fs=100000;
Wp=[20500/50000, 23500/50000]; % Band edge frequencies
Ws=[19000/50000, 25000/50000];
[N,Wc]=ellipord(Wp, Ws, Ap, As); % Determine filter order
[b, a]=ellip(N, Ap, As, Wc); % Determine filter coeffs
[z, p, k]=ellip(N, Ap, As,Wc); % Determine poles and zeros
sos=zp2sos(z, p,k); % Convert to second order sections
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, 20*log10(abs(H)))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response (dB)')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
============================================================
function IIRBZT
%Design program for IIR filters with Butterworth, Chebyshev, Inverse
%Chebyshev or Elliptic Characteristic using bilinear transformation
%output arguments:
%n - estimated minimun filter order.
%Wn - cutoff frequencies.
%b,a - coefficients of digital filter
%ftype - 'low' for lowpass filter
% 'high' for highpass filter
% 'bandpass' for bandpass filter
% 'stop' for bandstop filter
%ctype - 1 for Butterworth filter
% 2 for Chebyshev filter
% 3 for Inverse Chebyshev filter
% 4 for Elliptic filter
%steps of designing a digital filter implied in the called functions:
% step 1: estimate the minimum order of the filter from specifications
% step 2: get analog, pre-warped frequencies
% step 3: convert to low-pass prototype estimate
% step 4: Get n-th order analog lowpass prototype with desired filter characters
% step 5: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
% step 6: Use Bilinear transformation to find discrete equivalent:
clear all; format;
Fs = 100000; %sampling frequency(Hz).
Wp = [20500 23500]/(Fs/2); %passband edge frequency normalised by Fs/2.
Ws = [19000 25000]/(Fs/2); %stopband edge frewquency normalised by Fs/2.
Rp = 0.25; %passband attenuation in dB.
Rs = 45; %stopband attenuation in dB
ctype = 4; %character of filter
ftype = 'bandpass'; %type of filter
if ctype==1
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn,ftype);
elseif ctype==2
[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
[b,a]=cheby1(n,Rp,Wn,ftype);
elseif ctype==3
[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs);
[b,a]=cheby2(n,Rs,Wn,ftype);
elseif ctype==4
[n,Wn]=ellipord(Wp,Ws,Rp,Rs);
[b,a]=ellip(n,Rp,Rs,Wn,ftype);
end
%Output the result
disp('Numerator coefficients (in descending powers of z):'); disp(b);
disp('Denominator coefficients (in descending powers of z):'); disp(a);
freqz(b,a,1024,Fs);
=======================================================
function Impinv
%Impulse invariance method of anolog-to-digital filter conversion
%a,b -- s-plane coefficients
%az,bz -- digital filter coefficients
clear all;
b = 1;
a = [1.84496 1.920675 1];
[bz,az]=impinvar(b,a) %get z-plane coefficients using impulse Inv.
freqz(bz,az,1024);
=======================================================
% % Program name - prob825.m (problem 8.25(1), p553)
% Comparison of magnitude and phase responses of an analogue
% and equivalent BZT and impulse invariant discrete-time filters
% and pole-zero diagrams (elliptic low pass filter)
%
Fs=10000; FN=Fs/2;
fp=1000; fs=3000;
wp=fp*2*pi; ws=fs*2*pi;
Ap=1; As=60;
%
% Calculate filter coefficients and frequency responses
[N, wc]=ellipord(wp, ws, Ap, As,'s'); % analog filter
[B, A]=ellip(N, Ap, As, wc, 's');
[bBZT, aBZT]=bilinear(B,A,Fs); % BZT filter
[bIIT,aIIT]=impinvar(B,A,Fs); % Impulse invariance filter
%
% Compute frequency response
[Ha, wa]=freqs(B,A);
[HBZT, fBZT]=freqz(bBZT, aBZT, 512, Fs);
[HIIT, fIIT]=freqz(bIIT, aIIT, 512, Fs);
%
% Plot magnitude frequency responses
%
figure(1); % Plot analogue magnitude response
plot(wa/(2*pi), 20*log10(abs(Ha)))
hold on
figure (1);
plot(fBZT, 20*log10(abs(HBZT)), 'r:') % Plot BZT magnitude response
hold on
figure (1);
plot(fIIT, 20*log10(abs(HIIT)), 'g:') % Plot Impinv magnitude response
legend('Analog', 'BZT', 'Imp Invar');
axis([0 10000 -120 0])
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
title('Filter magnitude responses')
hold off;
%
% Plot phase responses
%
figure(2);
plot(wa/(2*pi), angle(Ha)*180/pi) % Plot analogue phase response
hold on
figure (2);
plot(fBZT, angle(HBZT)*180/pi, 'r:') % Plot BZT phase response
hold on
figure(2);
plot(fIIT, angle(HIIT)*180/pi, 'g.') %ImpInvar phase response
legend('Analog', 'BZT','Imp Invar');
axis([0 10000 -360 360])
ylabel('Phase (Degrees)')
xlabel('Frequency (Hz)')
title('Filter Phase Responses')
hold off
%
% Plot pole-zero diagrams
%
figure (3);
zplane(bBZT, aBZT)
title('Pole-zero diagram - BZT filter')
figure (4);
zplane(bIIT, aIIT)
xmin=-1; xmax=1; ymin=-1; ymax=1; % Scale the z-plane for the Impulse Inva.
axis([xmin xmax ymin ymax])
title('Pole-zero diagram - impulse invariance filter')
=====================================================
% % Program name - prob8261a.m (modified problem 8.26(1), p553)
% Comparison of characteristic features of discrete-time Butterworth,
% Chebyshev 1,Chebyshev 2 and elliptic lowpass filters.
%
Fs=8000;
fp=500/4000;
fs=2000/5000;
Ap=0.1; As=60;
%
% Calculate filter coefficients
[N1, wc1]=buttord(fp, fs, Ap, As); % filter order
[N2, wc2]=cheb1ord(fp, fs, Ap, As);
[N3, wc3]=cheb2ord(fp, fs, Ap, As);
[N4, wc4]=ellipord(fp, fs, Ap, As);
[b1, a1]=butter(4, fp); % Butterworth filter
[b2, a2]=cheby1(4, Ap, fp);
[b3, a3]=cheby2(4, As, fs);
[b4, a4]=ellip(4, Ap, As, fp);
%
% Calculate frequency responses
%
[H1, f1]=freqz(b1, a1, 512, Fs);
[H2, f2]=freqz(b2, a2, 512, Fs);
[H3, f3]=freqz(b3, a3, 512, Fs);
[H4, f4]=freqz(b4, a4, 512, Fs);
%
% Plot magnitude frequency responses
%
figure (1);
plot(f1, 20*log10(abs(H1)), 'r:') % Plot Butterworth magnitude response
hold on
figure (1);
plot(f2, 20*log10(abs(H2)), 'b--') % Plot Cheby1 magnitude response
hold on
figure (1);
plot(f3, 20*log10(abs(H3)), 'g-.') % Plot Cheby2 magnitude response
hold on
figure (1);
plot(f4, 20*log10(abs(H4)), 'k-') % Plot elliptic magnitude response
hold on
legend('Butworth', 'Cheby1', 'Cheby2', 'Ellip');
axis([0 4000 -100 10])
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
title('Filter magnitude responses')
hold off;
%
% Plot phase responses
%
figure (2);
plot(f1, angle(H1)*180/pi, 'r:') % Plot BZT phase response
hold on
figure (2);
plot(f2, angle(H2)*180/pi, 'b--') % Plot Cheby 1 phase response
hold on
figure (2);
plot(f3, angle(H3)*180/pi, 'g-.') % Plot Cheby 2 phase response
hold on
figure (2);
plot(f4, angle(H4)*180/pi, 'k-') % Plot elliptic phase response
hold on
figure(2);
legend('Butter', 'Cheby1','Cheby2', 'Ellip');
axis([0 4000 -360 360])
ylabel('Phase (Degrees)')
xlabel('Frequency (Hz)')
title('Filter Phase Responses')
hold off
%
% Plot pole-zero diagrams
%
figure (3);
zplane(b1, a1)
title('Pole-zero diagram - Butterworth filter')
figure (4);
zplane(b2, a2)
title('Pole-zero diagram - Chebyshev 1 filter')
figure (5);
zplane(b3, a3)
title('Pole-zero diagram - Chebyshev 2 filter')
figure (6);
zplane(b4, a4)
title('Pole-zero diagram - Elliptic filter')
N1
N2
N3
N4
% (Example 8.19) % program name: prog81.m
%
FN=1000/2; fc=300; % cut offf frequency
N=5; % filter order
[z, p, k]=buttap(N); % create an analog filter
w=linspace(0, FN/fc, 1000); % plot the response of filter
h=freqs(k*poly(z), poly(p), w);
f=fc*w;
plot(f, 20*log10(abs(h))), grid
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
=============================================
% % Program 8.2, p511. m-file for the design of the
% impulse invariant filter (Example 8.19) % program name: prog82.m
%
Fs=1000; % sampling frequency
fc=300; % cutoff frequency
WC=2*pi*fc; % cutoff frequency in radian
N=5; % filter order
[b,a]=butter(N,WC,'s'); % create an analog filter
[z, p, k]=butter(N, WC, 's');
[bz, az]=impinvar(b,a,Fs); % determine coeffs of IIR filter
[h, f]=freqz(bz, az, 512,Fs);
plot(f, 20*log10(abs(h))), grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
==============================================
% % Program 8.3, p512. m-file for the design of the BZT filter
% (Example 8.19) % program name: prog83.m %
Fs=1000; % sampling frequency
FN=Fs/2;
fc=300; % cutoff frequency
N=5;
[z, p, k]=butter(N, fc/FN);
[h, f]=freqz(k*poly(z), poly(p), 512, Fs);
plot(f, 20*log10(abs(h))), grid
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
=====================================================
% % m-file for Example 8B.1 (Program 8B.1, p564).
% Program name: prog8b1.m %
N=2; % Filter order
Fs=1280; % Sampling frequency
fc=150; % Cutoff frequency
WC=2*pi*fc; % Cutoff frequency in radian
%
%Create an analogue filter
%
[b,a]=butter(N,WC,'s');
[z,p,k]=butter(N, WC, 's');
%
%Convert analogue filter into Discrete IIR filter
%
[bz, az]=impinvar(b, a, Fs); %Determine coeffs of IIR filter
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(bz, az, 512, Fs);
plot(f, 20*log10(abs(H)))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response (dB)')
subplot(2,1,2) % Plot pole-zero diagram
zplane(bz, az)
zz=roots(bz); % Determine poles and zeros
pz=roots(az);
============================================================
% % m-file for Example 8B.2 (Program 8B.2, p566).
%Program name: prog8b2.m %
N=2; % Filter order
Fs=1280; % Sampling frequency
FN=Fs/2;
fc=150;% Cutoff frequency
Fc=fc/FN; % Normalized Cutoff frequency
[b,a]=butter(N,Fc); % Create and digitize analogue filter.
[z,p,k]=butter(N, Fc);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
b
a
==========================================================
% % m-file for Example 8B.3 (Program 8B.3, p567)
% Program name: prog8b3.m
%
Fs=2000; % Sampling frequency
FN=Fs/2;
fc1=200/FN;
fc2=300/FN;
[b,a]=butter(4,[fc1, fc2]); % Create and digitize analogue filter.
[z,p,k]=butter(4, [fc1, fc2]);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
=========================================================
% % m-file for Example 8B.4 (Program 8B.4, p569)
% Program name: prog8b4.m %
Fs=8000; % Sampling frequency
Ap=3;
As=20;
wp=500/4000;
ws=2000/4000;
[N, wc]=buttord(wp, ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,500/4000); % Digitise filter
[b, a]=butter(N, 500/4000);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
=======================================================
% % m-file for Example 8B.5 (Program 8B.5, p571)
% Program name: prog8b5.m
%
Fs=8000; % Sampling frequency
Ap=3;
As=20;
wp=2000/4000;
ws=500/4000;
[N, wc]=buttord(wp, ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,2000/4000, 'high'); % Digitise filter
[b, a]=butter(N, 2000/4000, 'high');
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
====================================================
% % m-file for Example 8B.6 (Program 8B.6, p572)
%Program name: prog8b6.m
%Band pass filter
%
Fs=1000; % Sampling frequency
Ap=3;
As=20;
Wp=[200/500, 300/500]; % Band edge frequencies
Ws=[50/500, 450/500];
[N, Wc]=buttord(Wp, Ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,Wp); % Digitise filter
[b, a]=butter(N, Wp);
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
=======================================================
% % m-file for Example 8B.7 (Program 8B.7, p574)
%Program name: prog8b7.m % Band stop filter
%
Fs=1000; %Sampling frequency
Ap=3;
As=20;
Wp=[50/500, 450/500]; % Band edge frequencies
Ws=[200/500, 300/500];
[N, Wc]=buttord(Wp, Ws, Ap, As); % Determine filter order
[zz, pz, kz]=butter(N,Ws, 'stop'); % Digitise filter
[b, a]=butter(N, Ws, 'stop');
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, abs(H))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response ')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
============================================================
% % m-file for the design Example 8B.8 (Program 8B.8, p576)
%Program Name: prog8b8.m %
Ap=0.25;
As=45;
Fs=100000;
Wp=[20500/50000, 23500/50000]; % Band edge frequencies
Ws=[19000/50000, 25000/50000];
[N,Wc]=ellipord(Wp, Ws, Ap, As); % Determine filter order
[b, a]=ellip(N, Ap, As, Wc); % Determine filter coeffs
[z, p, k]=ellip(N, Ap, As,Wc); % Determine poles and zeros
sos=zp2sos(z, p,k); % Convert to second order sections
subplot(2,1,1) % Plot magnitude freq. response
[H, f]=freqz(b, a, 512, Fs);
plot(f, 20*log10(abs(H)))
xlabel('Frequency (Hz)')
ylabel('Magnitude Response (dB)')
subplot(2,1,2) % Plot pole-zero diagram
zplane(b, a)
============================================================
function IIRBZT
%Design program for IIR filters with Butterworth, Chebyshev, Inverse
%Chebyshev or Elliptic Characteristic using bilinear transformation
%output arguments:
%n - estimated minimun filter order.
%Wn - cutoff frequencies.
%b,a - coefficients of digital filter
%ftype - 'low' for lowpass filter
% 'high' for highpass filter
% 'bandpass' for bandpass filter
% 'stop' for bandstop filter
%ctype - 1 for Butterworth filter
% 2 for Chebyshev filter
% 3 for Inverse Chebyshev filter
% 4 for Elliptic filter
%steps of designing a digital filter implied in the called functions:
% step 1: estimate the minimum order of the filter from specifications
% step 2: get analog, pre-warped frequencies
% step 3: convert to low-pass prototype estimate
% step 4: Get n-th order analog lowpass prototype with desired filter characters
% step 5: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
% step 6: Use Bilinear transformation to find discrete equivalent:
clear all; format;
Fs = 100000; %sampling frequency(Hz).
Wp = [20500 23500]/(Fs/2); %passband edge frequency normalised by Fs/2.
Ws = [19000 25000]/(Fs/2); %stopband edge frewquency normalised by Fs/2.
Rp = 0.25; %passband attenuation in dB.
Rs = 45; %stopband attenuation in dB
ctype = 4; %character of filter
ftype = 'bandpass'; %type of filter
if ctype==1
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn,ftype);
elseif ctype==2
[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
[b,a]=cheby1(n,Rp,Wn,ftype);
elseif ctype==3
[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs);
[b,a]=cheby2(n,Rs,Wn,ftype);
elseif ctype==4
[n,Wn]=ellipord(Wp,Ws,Rp,Rs);
[b,a]=ellip(n,Rp,Rs,Wn,ftype);
end
%Output the result
disp('Numerator coefficients (in descending powers of z):'); disp(b);
disp('Denominator coefficients (in descending powers of z):'); disp(a);
freqz(b,a,1024,Fs);
=======================================================
function Impinv
%Impulse invariance method of anolog-to-digital filter conversion
%a,b -- s-plane coefficients
%az,bz -- digital filter coefficients
clear all;
b = 1;
a = [1.84496 1.920675 1];
[bz,az]=impinvar(b,a) %get z-plane coefficients using impulse Inv.
freqz(bz,az,1024);
=======================================================
% % Program name - prob825.m (problem 8.25(1), p553)
% Comparison of magnitude and phase responses of an analogue
% and equivalent BZT and impulse invariant discrete-time filters
% and pole-zero diagrams (elliptic low pass filter)
%
Fs=10000; FN=Fs/2;
fp=1000; fs=3000;
wp=fp*2*pi; ws=fs*2*pi;
Ap=1; As=60;
%
% Calculate filter coefficients and frequency responses
[N, wc]=ellipord(wp, ws, Ap, As,'s'); % analog filter
[B, A]=ellip(N, Ap, As, wc, 's');
[bBZT, aBZT]=bilinear(B,A,Fs); % BZT filter
[bIIT,aIIT]=impinvar(B,A,Fs); % Impulse invariance filter
%
% Compute frequency response
[Ha, wa]=freqs(B,A);
[HBZT, fBZT]=freqz(bBZT, aBZT, 512, Fs);
[HIIT, fIIT]=freqz(bIIT, aIIT, 512, Fs);
%
% Plot magnitude frequency responses
%
figure(1); % Plot analogue magnitude response
plot(wa/(2*pi), 20*log10(abs(Ha)))
hold on
figure (1);
plot(fBZT, 20*log10(abs(HBZT)), 'r:') % Plot BZT magnitude response
hold on
figure (1);
plot(fIIT, 20*log10(abs(HIIT)), 'g:') % Plot Impinv magnitude response
legend('Analog', 'BZT', 'Imp Invar');
axis([0 10000 -120 0])
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
title('Filter magnitude responses')
hold off;
%
% Plot phase responses
%
figure(2);
plot(wa/(2*pi), angle(Ha)*180/pi) % Plot analogue phase response
hold on
figure (2);
plot(fBZT, angle(HBZT)*180/pi, 'r:') % Plot BZT phase response
hold on
figure(2);
plot(fIIT, angle(HIIT)*180/pi, 'g.') %ImpInvar phase response
legend('Analog', 'BZT','Imp Invar');
axis([0 10000 -360 360])
ylabel('Phase (Degrees)')
xlabel('Frequency (Hz)')
title('Filter Phase Responses')
hold off
%
% Plot pole-zero diagrams
%
figure (3);
zplane(bBZT, aBZT)
title('Pole-zero diagram - BZT filter')
figure (4);
zplane(bIIT, aIIT)
xmin=-1; xmax=1; ymin=-1; ymax=1; % Scale the z-plane for the Impulse Inva.
axis([xmin xmax ymin ymax])
title('Pole-zero diagram - impulse invariance filter')
=====================================================
% % Program name - prob8261a.m (modified problem 8.26(1), p553)
% Comparison of characteristic features of discrete-time Butterworth,
% Chebyshev 1,Chebyshev 2 and elliptic lowpass filters.
%
Fs=8000;
fp=500/4000;
fs=2000/5000;
Ap=0.1; As=60;
%
% Calculate filter coefficients
[N1, wc1]=buttord(fp, fs, Ap, As); % filter order
[N2, wc2]=cheb1ord(fp, fs, Ap, As);
[N3, wc3]=cheb2ord(fp, fs, Ap, As);
[N4, wc4]=ellipord(fp, fs, Ap, As);
[b1, a1]=butter(4, fp); % Butterworth filter
[b2, a2]=cheby1(4, Ap, fp);
[b3, a3]=cheby2(4, As, fs);
[b4, a4]=ellip(4, Ap, As, fp);
%
% Calculate frequency responses
%
[H1, f1]=freqz(b1, a1, 512, Fs);
[H2, f2]=freqz(b2, a2, 512, Fs);
[H3, f3]=freqz(b3, a3, 512, Fs);
[H4, f4]=freqz(b4, a4, 512, Fs);
%
% Plot magnitude frequency responses
%
figure (1);
plot(f1, 20*log10(abs(H1)), 'r:') % Plot Butterworth magnitude response
hold on
figure (1);
plot(f2, 20*log10(abs(H2)), 'b--') % Plot Cheby1 magnitude response
hold on
figure (1);
plot(f3, 20*log10(abs(H3)), 'g-.') % Plot Cheby2 magnitude response
hold on
figure (1);
plot(f4, 20*log10(abs(H4)), 'k-') % Plot elliptic magnitude response
hold on
legend('Butworth', 'Cheby1', 'Cheby2', 'Ellip');
axis([0 4000 -100 10])
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
title('Filter magnitude responses')
hold off;
%
% Plot phase responses
%
figure (2);
plot(f1, angle(H1)*180/pi, 'r:') % Plot BZT phase response
hold on
figure (2);
plot(f2, angle(H2)*180/pi, 'b--') % Plot Cheby 1 phase response
hold on
figure (2);
plot(f3, angle(H3)*180/pi, 'g-.') % Plot Cheby 2 phase response
hold on
figure (2);
plot(f4, angle(H4)*180/pi, 'k-') % Plot elliptic phase response
hold on
figure(2);
legend('Butter', 'Cheby1','Cheby2', 'Ellip');
axis([0 4000 -360 360])
ylabel('Phase (Degrees)')
xlabel('Frequency (Hz)')
title('Filter Phase Responses')
hold off
%
% Plot pole-zero diagrams
%
figure (3);
zplane(b1, a1)
title('Pole-zero diagram - Butterworth filter')
figure (4);
zplane(b2, a2)
title('Pole-zero diagram - Chebyshev 1 filter')
figure (5);
zplane(b3, a3)
title('Pole-zero diagram - Chebyshev 2 filter')
figure (6);
zplane(b4, a4)
title('Pole-zero diagram - Elliptic filter')
N1
N2
N3
N4
No comments:
Post a Comment