octave_packages/signal-1.1.3/cheby1.m

   1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
   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 ## Generate an Chebyshev type I filter with Rp dB of pass band ripple.
  18 ##
  19 ## [b, a] = cheby1(n, Rp, Wc)
  20 ##    low pass filter with cutoff pi*Wc radians
  21 ##
  22 ## [b, a] = cheby1(n, Rp, Wc, 'high')
  23 ##    high pass filter with cutoff pi*Wc radians
  24 ##
  25 ## [b, a] = cheby1(n, Rp, [Wl, Wh])
  26 ##    band pass filter with edges pi*Wl and pi*Wh radians
  27 ##
  28 ## [b, a] = cheby1(n, Rp, [Wl, Wh], 'stop')
  29 ##    band reject filter with edges pi*Wl and pi*Wh radians
  30 ##
  31 ## [z, p, g] = cheby1(...)
  32 ##    return filter as zero-pole-gain rather than coefficients of the
  33 ##    numerator and denominator polynomials.
  34 ##
  35 ## [...] = cheby1(...,'s')
  36 ##     return a Laplace space filter, W can be larger than 1.
  37 ##
  38 ## [a,b,c,d] = cheby1(...)
  39 ##  return  state-space matrices
  40 ##
  41 ## References:
  42 ##
  43 ## Parks & Burrus (1987). Digital Filter Design. New York:
  44 ## John Wiley & Sons, Inc.
  45
  46 function [a,b,c,d] = cheby1(n, Rp, W, varargin)
  47
  48   if (nargin>5 || nargin<3) || (nargout>4 || nargout<2)
  49     print_usage;
  50   endif
  51
  52   ## interpret the input parameters
  53   if (!(length(n)==1 && n == round(n) && n > 0))
  54     error ("cheby1: filter order n must be a positive integer");
  55   endif
  56
  57   stop = 0;
  58   digital = 1;
  59   for i=1:length(varargin)
  60     switch varargin{i}
  61     case 's', digital = 0;
  62     case 'z', digital = 1;
  63     case { 'high', 'stop' }, stop = 1;
  64     case { 'low',  'pass' }, stop = 0;
  65     otherwise,  error ("cheby1: expected [high|stop] or [s|z]");
  66     endswitch
  67   endfor
  68
  69   [r, c]=size(W);
  70   if (!(length(W)<=2 && (r==1 || c==1)))
  71     error ("cheby1: frequency must be given as w0 or [w0, w1]");
  72   elseif (!(length(W)==1 || length(W) == 2))
  73     error ("cheby1: only one filter band allowed");
  74   elseif (length(W)==2 && !(W(1) < W(2)))
  75     error ("cheby1: first band edge must be smaller than second");
  76   endif
  77
  78   if ( digital && !all(W >= 0 & W <= 1))
  79     error ("cheby1: critical frequencies must be in (0 1)");
  80   elseif ( !digital && !all(W >= 0 ))
  81     error ("cheby1: critical frequencies must be in (0 inf)");
  82   endif
  83
  84   if (Rp < 0)
  85     error("cheby1: passband ripple must be positive decibels");
  86   end
  87
  88   ## Prewarp to the band edges to s plane
  89   if digital
  90     T = 2;       # sampling frequency of 2 Hz
  91     W = 2/T*tan(pi*W/T);
  92   endif
  93
  94   ## Generate splane poles and zeros for the chebyshev type 1 filter
  95   C = 1; # default cutoff frequency
  96   epsilon = sqrt(10^(Rp/10) - 1);
  97   v0 = asinh(1/epsilon)/n;
  98   pole = exp(1i*pi*[-(n-1):2:(n-1)]/(2*n));
  99   pole = -sinh(v0)*real(pole) + 1i*cosh(v0)*imag(pole);
 100   zero = [];
 101
 102   ## compensate for amplitude at s=0
 103   gain = prod(-pole);
 104   ## if n is even, the ripple starts low, but if n is odd the ripple
 105   ## starts high. We must adjust the s=0 amplitude to compensate.
 106   if (rem(n,2)==0)
 107     gain = gain/10^(Rp/20);
 108   endif
 109
 110   ## splane frequency transform
 111   [zero, pole, gain] = sftrans(zero, pole, gain, W, stop);
 112
 113   ## Use bilinear transform to convert poles to the z plane
 114   if digital
 115     [zero, pole, gain] = bilinear(zero, pole, gain, T);
 116   endif
 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