% This script illustrates basic processing and plotting functions
% for examining filter properties, specifying filter properties for
% the IIR filters in Matlab.  This script uses the IIR filters
% BUTTER (Butterworth maximally flat), CHEBY1 (Chebyshev Type 1 filter,
% ripple in passband), CHEBY2 (Chebyshev Type 2 filter, ripple in 
% stopband) and ELLIP (Elliptical filter, ripple in both pass and stop
% bands).  In general for a fixed order, the more ripple allowed for,
% the sharper the transition band.  FREQZ is used to compute the filter
% transfer functions, FILTER is used to compute the impulse response,
% and ROOTS to compute zero and pole locations of the filters.
% All examples will show the filter design and analysis for a
% 6th order band pass filter with a sampling rate of 8000Hz and
% a passband between 1000 and 2500 Hz, or in the case where the
% stopband edge frequencies are required, they will be taken at
% 850 and 2900 Hz.
%  A statement is used after the transfer function plots
%   >> set(gca,'Ylim', [-60 3])
%  This sets the limits on the Y-axis (for the current plot, gca stands
%  for "get current axis") so plot details can better be seen.  There
%  are regions in the stopband that go as low as -500 dB. If these are
%  plotted the small ripple in the passbands cannot be observed.
%
%  Written by Kevin D. Donohue ( donohue@engr.uky.edu ) Oct. 4, 2007 
 
 
 
clear
% General filter specification
n = 6;  %  Filter order
fs = 8000;  %  Sampling frequency
impdur = .01;  %  Impulse response duration in seconds
 
% filter edge freqencies
wpband = [1000 2500];  %  Lower and upper cutoff frequencies for passband edges
wsband = [850 2900];  %  Lower and upper stopband edge frequencies (used in cheby2).
%  Note if a band pass or band stop filter is requested, the cut-off
%  frequency vector will have 2 elements to denote the upper and lower
%  cut-off values.  If 2 element are present, the filter order will double
%  the requested value (so you will see twice as many coefficients as
%  expected.
 
 
%  Maximum ripple values 
rp = 3;  %  passband ripple in dB
rs = 40; %  stopband ripple in dB
 
 
%  Butterworth filter is maximually flat.  There is no ripple in either
%  stopband or passband.  If a ripple value is given, it is equal to the 
%  difference between the maximum and minimum value in the passband
%  or the maximum peak in the stop band.  The butterworth filter
%  funtion in Matlab only requires the filter order and passband edge
%  frequencies.
 
%  A bandpass filter is assume if wpband is a vector of 2 elements
[bbut, abut] = butter(n,wpband/(fs/2)); 
 
%  Compute impulse response
imps = [1, zeros(1,round(impdur*fs)-1)];  %  Impulse signal input
ybut = filter(bbut,abut,imps);   %  Impulse response (windowing method)
 
%  Plot impuse response
figure(1)
plot([0:length(ybut)-1]/fs,ybut);  %  Normalize amplitude for plot
xlabel('seconds')
ylabel('Amplitude')
title('Butterworth impulse response')
 
%  Plot Transfer Functions
figure(2)
[h1,f1] = freqz(bbut,abut,1024,fs);
plot(f1,20*log10(abs(h1)))
ylabel('dB')
xlabel('Hz')
title('Butterwork TF magnitudes')
set(gca,'Ylim', [-60 3]);  %  Zoom in on Y-axis, set limit to -50 and 3
% Plot pole-zero diagram
figure(3)
zz = roots(bbut);  % Zeros
zp = roots(abut);  % Poles
%  Generate unit circle for reference
w = [-1:.001:1]*pi;
cpts = exp(j*w);
%  Plot unit circle with plots and zeros
plot(real(cpts),imag(cpts),'k-',real(zz),imag(zz),'bo',real(zp),imag(zp),'rx')
xlabel('Real')
ylabel('Imaginary')
title('Butterworth, Pole-Zeros Diagram')
 
%  Chebyshev type 1 filter is flat in the stopband, but has ripple in the passband.
%  The Chebyshev type 1 filter funtion in Matlab requires the 
%  filter order and passband edge frequencies, and passband ripple value
%  (difference between max and min TF magnitude value in the passband)
 
%  A bandpass filter is assume if wpband is a vector of 2 elements
[bche1, ache1] = cheby1(n, rp, wpband/(fs/2)); 
 
%  Compute impulse response
imps = [1, zeros(1,round(impdur*fs)-1)];  %  Impulse signal input
yche1 = filter(bche1,ache1,imps);   %  Impulse response (windowing method)
 
%  Plot impuse response
figure(4)
plot([0:length(yche1)-1]/fs,yche1);  %  Normalize amplitude for plot
xlabel('seconds')
ylabel('Amplitude')
title('Chebyshev Type 1 impulse response')
 
%  Plot Transfer Functions
figure(5)
[h1,f1] = freqz(bche1,ache1,1024,fs);
plot(f1,20*log10(abs(h1)))
ylabel('dB')
xlabel('Hz')
title('Chebyshev Type 1 TF magnitudes')
set(gca,'Ylim', [-60 3]);  %  Zoom in on Y-axis, set limit to -50 and 3
 
% Plot pole-zero diagram
figure(6)
zz = roots(bche1);  % Zeros
zp = roots(ache1);  % Poles
%  Generate unit circle for reference
w = [-1:.001:1]*pi;
cpts = exp(j*w);
%  Plot unit circle with plots and zeros
plot(real(cpts),imag(cpts),'k-',real(zz),imag(zz),'bo',real(zp),imag(zp),'rx')
xlabel('Real')
ylabel('Imaginary')
title('Chebyshev Type 1, Pole-Zeros Diagram')
 
 
%  Chebyshev type 2 filter is flat in the passband, but has ripple in the stopband.
%  The Chebyshev type 1 filter funtion in Matlab requires the 
%  filter order and passband edge frequencies, and stopband ripple value
%  (peak value in the stopband)
 
%  A bandpass filter is assume if wpband is a vector of 2 elements
[bche2, ache2] = cheby2(n, rs, wsband/(fs/2)); 
 
%  Compute impulse response
imps = [1, zeros(1,round(impdur*fs)-1)];  %  Impulse signal input
yche2 = filter(bche2,ache2,imps);   %  Impulse response (windowing method)
 
%  Plot impuse response
figure(7)
plot([0:length(yche2)-1]/fs,yche2);  %  Normalize amplitude for plot
xlabel('seconds')
ylabel('Amplitude')
title('Chebyshev Type 2 impulse response')
 
%  Plot Transfer Functions
figure(8)
[h1,f1] = freqz(bche2,ache2,1024,fs);
plot(f1,20*log10(abs(h1)))
ylabel('dB')
xlabel('Hz')
title('Chebyshev Type 2 TF magnitudes')
set(gca,'Ylim', [-60 3]);  %  Zoom in on Y-axis, set limit to -50 and 3
 
% Plot pole-zero diagram
figure(9)
zz = roots(bche2);  % Zeros
zp = roots(ache2);  % Poles
%  Generate unit circle for reference
w = [-1:.001:1]*pi;
cpts = exp(j*w);
%  Plot unit circle with plots and zeros
plot(real(cpts),imag(cpts),'k-',real(zz),imag(zz),'bo',real(zp),imag(zp),'rx')
xlabel('Real')
ylabel('Imaginary')
title('Chebyshev Type 1, Pole-Zeros Diagram')
 
 
%  Elliptical filter has ripple in both the pass and stopband.
%  The elliptical filter funtion in Matlab requires the 
%  filter order and passband edge frequencies, and stopband ripple value
%  (peak value in the stopband)
 
%  A bandpass filter is assume if wpband is a vector of 2 elements
[bell, aell] = ellip(n, rp, rs, wsband/(fs/2)); 
 
%  Compute impulse response
imps = [1, zeros(1,round(impdur*fs)-1)];  %  Impulse signal input
yell = filter(bell,aell,imps);   %  Impulse response (windowing method)
 
%  Plot impuse response
figure(10)
plot([0:length(yell)-1]/fs,yell);  %  Normalize amplitude for plot
xlabel('seconds')
ylabel('Amplitude')
title('Elliptical impulse response')
 
%  Plot Transfer Functions
figure(11)
[h1,f1] = freqz(bell,aell,1024,fs);
plot(f1,20*log10(abs(h1)))
ylabel('dB')
xlabel('Hz')
title('Elliptical TF magnitudes')
set(gca,'Ylim', [-60 3]);  %  Zoom in on Y-axis, set limit to -50 and 3
 
% Plot pole-zero diagram
figure(12)
zz = roots(bell);  % Zeros
zp = roots(aell);  % Poles
%  Generate unit circle for reference
w = [-1:.001:1]*pi;
cpts = exp(j*w);
%  Plot unit circle with plots and zeros
plot(real(cpts),imag(cpts),'k-',real(zz),imag(zz),'bo',real(zp),imag(zp),'rx')
xlabel('Real')
ylabel('Imaginary')
title('Elliptical filter, Pole-Zeros Diagram')