Filter Design

Example of the design of a two stage, second order infinite impulse response(IIR):

% ADC sample rate:
Fs = 22050;   % Hz
 
% Filter Response
%    'low'     : Low pass 
%    'high'    : High Pass 
%    'bandpass': Band Pass 
%    'stop'    : Band Stop 
response = 'high';
 
% Wp is the pass band edge and Ws is the stop band edge
Wp = 993;      % Hz   Pass band corner frequency
Ws = 11000;    % Hz   Stop band corner frequency
 
Rp = 3;       % dB   Pass band Ripple
Rs = 60;      % dB   Stop band Ripple
 
% Desired Filter Gain
G = 1;
 
% Filter type:
%    Butterworth
%    Chebyshev_I
%    Chebyshev_II
%    Elliptic
type = 'Butterworth';
 
% Frequencies are normalized to [0,1], corresponding to the range [0,Fnyq]
Fnyq = Fs/2;  % Hz   Nyquist Frequency
Wp = Wp/Fnyq;       
Ws = Ws/Fnyq;
 
switch (type)
    case 'Butterworth'
        [order,Wn] = buttord(Wp,Ws,Rp,Rs);
    case 'Chebyshev_I'
        [order,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
    case 'Chebyshev_II'
        [order,Wn]=cheb2ord(Wp,Ws,Rp,Rs);
    case 'Elliptic'
        [order,Wn]=ellipord(Wp,Ws,Rp,Rs);
    otherwise
        disp('ERROR: unsupported filter type')
        return
end
 
%    - Two stage second order => 4th Order
if (order > 4)
    disp('It is not possible to design a two stage second order')
    disp('IIR filter with the following design specifications:')
    disp('Wp: ',Wp)
    disp('Ws: ',Ws)
    disp('Rp: ',Rp)
    disp('Rs: ',Rs)
    return
else
    order = 4;
end
 
switch (type)
    case 'Butterworth'
        [Z,P,K] = butter(order,Wn,response);
    case 'Chebyshev_I'
        [Z,P,K]=cheby1(order,Rp,Wn,response);
    case 'Chebyshev_II'
        [Z,P,K]=cheby2(order,Rs,Wn,response);
    case 'Elliptic'
        [Z,P,K]=ellip(order,Rp,Rs,Wn,response);
end
 
% converts a discrete-time zero-pole-gain representation 
% to an equivalent second-order section representation
sos = zp2sos(Z,P,K);
 
% Number of biquads (no of rows)
nbiq = size(sos,1);
if (nbiq ~= 2)
   disp('ERROR: Wrong number of biquads for a two stage second order filter')
   disp(nbiq)
   return
end
 
disp(sprintf('\nMatlab second-order section representation'));
disp(sos)
 
% System analysis
[B, A] = sos2tf(sos);
[H, f] = freqz(B, A, 512, Fs);
sys = tf(B, A, 1/Fs);
 
disp(sprintf('\nPoles distances from the unit circle:'));
disp(abs(eig(sys)));
 
S = allmargin(sys);
if (S.Stable)
  disp('THE FILTER IS STABLE');
else
  disp('THE FILTER IS UNSTABLE');
end
 
figure('Name','Frequency Response')
  plot(f,abs(H));
  grid;
  title('Frequency Response')
  xlabel('Hertz');
  ylabel('Magnitude Response');
 
figure('Name','Stability Analysis')
  % Pole-Zero Map
  subplot(2,1,1);
  pzmap(sys);
 
  % Step Response
  subplot(2,1,2);
  step(sys);
  grid;

dsp/matlab/filter_design.txt · Last modified: 2010/08/10 (external edit)
CC Attribution-Noncommercial-Share Alike 3.0 Unported
Valid CSS Driven by DokuWiki Recent changes RSS feed Valid XHTML 1.0