octave_packages/signal-1.1.3/besselap.m

   1 ## Copyright (C) 2009 Thomas Sailer
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## Return bessel analog filter prototype.
  17 ##
  18 ## References:
  19 ##
  20 ## http://en.wikipedia.org/wiki/Bessel_polynomials
  21
  22 function [zero, pole, gain] = besselap (n)
  23
  24   if (nargin>1 || nargin<1)
  25     print_usage;
  26   end
  27
  28   ## interpret the input parameters
  29   if (!(length(n)==1 && n == round(n) && n > 0))
  30     error ("besselap: filter order n must be a positive integer");
  31   end
  32
  33   p0=1;
  34   p1=[1 1];
  35   for nn=2:n;
  36     px=(2*nn-1)*p1;
  37     py=[p0 0 0];
  38     px=prepad(px,max(length(px),length(py)),0);
  39     py=prepad(py,length(px));
  40     p0=p1;
  41     p1=px+py;
  42   endfor;
  43   % p1 now contains the reverse bessel polynomial for n
  44
  45   % scale it by replacing s->s/w0 so that the gain becomes 1
  46   p1=p1.*p1(length(p1)).^((length(p1)-1:-1:0)/(length(p1)-1));
  47
  48   zero=[];
  49   pole=roots(p1);
  50   gain=1;
  51
  52 endfunction