octave_packages/signal-1.1.3/ellip.m

   1 ## Copyright (C) 2001 Paulo Neis <p_neis@yahoo.com.br>
   2 ## Copyright (C) 2003 Doug Stewart <dastew@sympatico.ca>
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## N-ellip 0.2.1
  18 ##usage: [Zz, Zp, Zg] = ellip(n, Rp, Rs, Wp, stype,'s')
  19 ##
  20 ## Generate an Elliptic or Cauer filter (discrete and contnuious).
  21 ##
  22 ## [b,a] = ellip(n, Rp, Rs, Wp)
  23 ##  low pass filter with order n, cutoff pi*Wp radians, Rp decibels
  24 ##  of ripple in the passband and a stopband Rs decibels down.
  25 ##
  26 ## [b,a] = ellip(n, Rp, Rs, Wp, 'high')
  27 ##  high pass filter with cutoff pi*Wp...
  28 ##
  29 ## [b,a] = ellip(n, Rp, Rs, [Wl, Wh])
  30 ##  band pass filter with band pass edges pi*Wl and pi*Wh ...
  31 ##
  32 ## [b,a] = ellip(n, Rp, Rs, [Wl, Wh], 'stop')
  33 ##  band reject filter with edges pi*Wl and pi*Wh, ...
  34 ##
  35 ## [z,p,g] = ellip(...)
  36 ##  return filter as zero-pole-gain.
  37 ##
  38 ## [...] = ellip(...,'s')
  39 ##     return a Laplace space filter, W can be larger than 1.
  40 ##
  41 ## [a,b,c,d] = ellip(...)
  42 ##  return  state-space matrices
  43 ##
  44 ## References:
  45 ##
  46 ## - Oppenheim, Alan V., Discrete Time Signal Processing, Hardcover, 1999.
  47 ## - Parente Ribeiro, E., Notas de aula da disciplina TE498 -  Processamento
  48 ##   Digital de Sinais, UFPR, 2001/2002.
  49 ## - Kienzle, Paul, functions from Octave-Forge, 1999 (http://octave.sf.net).
  50
  51
  52 function [a,b,c,d] = ellip(n, Rp, Rs, W, varargin)
  53
  54   if (nargin>6 || nargin<4) || (nargout>4 || nargout<2)
  55     print_usage;
  56   endif
  57
  58   ## interpret the input parameters
  59   if (!(length(n)==1 && n == round(n) && n > 0))
  60     error ("ellip: filter order n must be a positive integer");
  61   endif
  62
  63
  64   stop = 0;
  65   digital = 1;
  66   for i=1:length(varargin)
  67     switch varargin{i}
  68     case 's', digital = 0;
  69     case 'z', digital = 1;
  70     case { 'high', 'stop' }, stop = 1;
  71     case { 'low',  'pass' }, stop = 0;
  72     otherwise,  error ("ellip: expected [high|stop] or [s|z]");
  73     endswitch
  74   endfor
  75
  76   [r, c]=size(W);
  77   if (!(length(W)<=2 && (r==1 || c==1)))
  78     error ("ellip: frequency must be given as w0 or [w0, w1]");
  79   elseif (!(length(W)==1 || length(W) == 2))
  80     error ("ellip: only one filter band allowed");
  81   elseif (length(W)==2 && !(W(1) < W(2)))
  82     error ("ellip: first band edge must be smaller than second");
  83   endif
  84
  85   if ( digital && !all(W >= 0 & W <= 1))
  86     error ("ellip: critical frequencies must be in (0 1)");
  87   elseif ( !digital && !all(W >= 0 ))
  88     error ("ellip: critical frequencies must be in (0 inf)");
  89   endif
  90
  91   if (Rp < 0)
  92     error("ellip: passband ripple must be positive decibels");
  93   endif
  94
  95   if (Rs < 0)
  96     error("ellip: stopband ripple must be positive decibels");
  97   end
  98
  99
 100   ##Prewarp the digital frequencies
 101   if digital
 102     T = 2;       # sampling frequency of 2 Hz
 103     W = tan(pi*W/T);
 104   endif
 105
 106   ##Generate s-plane poles, zeros and gain
 107   [zero, pole, gain] = ncauer(Rp, Rs, n);
 108
 109   ## splane frequency transform
 110   [zero, pole, gain] = sftrans(zero, pole, gain, W, stop);
 111
 112   ## Use bilinear transform to convert poles to the z plane
 113   if digital
 114     [zero, pole, gain] = bilinear(zero, pole, gain, T);
 115   endif
 116
 117
 118   ## convert to the correct output form
 119   if nargout==2,
 120     a = real(gain*poly(zero));
 121     b = real(poly(pole));
 122   elseif nargout==3,
 123     a = zero;
 124     b = pole;
 125     c = gain;
 126   else
 127     ## output ss results
 128     [a, b, c, d] = zp2ss (zero, pole, gain);
 129   endif
 130
 131 endfunction
 132
 133 %!demo
 134 %! clc
 135 %! disp('---------------------------> NELLIP 0.2 EXAMPLE <-------------------------')
 136 %! x=input("Let's calculate the filter order: [ENTER]");
 137 %! disp("")
 138 %! x=input("[n, Ws] = ellipord([.1 .2],.4,1,90); [ENTER]");
 139 %! [n, Ws] = ellipord([.1 .2],.4,1,90)
 140 %! disp("")
 141 %! x=input("Let's calculate the filter: [ENTER]");
 142 %! disp("")
 143 %! x=input("[b,a]=ellip(5,1,90,[.1,.2]);  [ENTER]");
 144 %! [b,a]=ellip(5,1,90,[.1,.2])
 145 %! disp("")
 146 %! x=input("Let's calculate the frequency response: [ENTER]");
 147 %! disp("")
 148 %! x=input("[h,w]=freqz(b,a);  [ENTER]");
 149 %! [h,w]=freqz(b,a);
 150 %!
 151 %! xlabel("Frequency");
 152 %! ylabel("abs(H[w])[dB]");
 153 %! axis([0,1,-100,0]);
 154 %! plot(w./pi, 20*log10(abs(h)), ';;')
 155 %!
 156 %! hold('on');
 157 %! x=ones(1,length(h));
 158 %! plot(w./pi, x.*-1, ';-1 dB;')
 159 %! plot(w./pi, x.*-90, ';-90 dB;')
 160 %! hold('off');
 161 %!
 162 %! xlabel("")
 163 %! ylabel("")
 164 %! clc