/usr/share/octave/packages/signal-1.3.0/besself.m is in octave-signal 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 | ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
## Copyright (C) 2003 Doug Stewart <dastew@sympatico.ca>
## Copyright (C) 2009 Thomas Sailer <t.sailer@alumni.ethz.ch>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{a}, @var{b}] =} besself (@var{n}, @var{W})
## @deftypefnx {Function File} {[@var{a}, @var{b}] =} besself (@var{n}, @var{W}, "high")
## @deftypefnx {Function File} {[@var{a}, @var{b}] =} besself (@var{n}, [@var{Wl}, @var{Wh}])
## @deftypefnx {Function File} {[@var{a}, @var{b}] =} besself (@var{n}, [@var{Wl}, @var{Wh}], "stop")
## @deftypefnx {Function File} {[@var{z}, @var{p}, @var{g}] =} besself (@dots{})
## @deftypefnx {Function File} {[@var{a}, @var{b}, @var{c}, @var{d}] =} besself (@dots{})
## @deftypefnx {Function File} {[@dots{}] =} besself (@dots{}, "z")
## Generate a bessel filter.
## Default is a Laplace space (s) filter.
##
## [b,a] = besself(n, Wc)
## low pass filter with cutoff pi*Wc radians
##
## [b,a] = besself(n, Wc, 'high')
## high pass filter with cutoff pi*Wc radians
##
## [b,a] = besself(n, [Wl, Wh])
## band pass filter with edges pi*Wl and pi*Wh radians
##
## [b,a] = besself(n, [Wl, Wh], 'stop')
## band reject filter with edges pi*Wl and pi*Wh radians
##
## [z,p,g] = besself(...)
## return filter as zero-pole-gain rather than coefficients of the
## numerator and denominator polynomials.
##
## [...] = besself(...,'z')
## return a discrete space (Z) filter, W must be less than 1.
##
## [a,b,c,d] = besself(...)
## return state-space matrices
##
## References:
##
## Proakis & Manolakis (1992). Digital Signal Processing. New York:
## Macmillan Publishing Company.
## @end deftypefn
function [a, b, c, d] = besself (n, W, varargin)
if (nargin>4 || nargin<2) || (nargout>4 || nargout<2)
print_usage;
endif
## interpret the input parameters
if (!(length(n)==1 && n == round(n) && n > 0))
error ("besself: filter order n must be a positive integer");
endif
stop = 0;
digital = 0;
for i=1:length(varargin)
switch varargin{i}
case 's', digital = 0;
case 'z', digital = 1;
case { 'high', 'stop' }, stop = 1;
case { 'low', 'pass' }, stop = 0;
otherwise, error ("besself: expected [high|stop] or [s|z]");
endswitch
endfor
[r, c]=size(W);
if (!(length(W)<=2 && (r==1 || c==1)))
error ("besself: frequency must be given as w0 or [w0, w1]");
elseif (!(length(W)==1 || length(W) == 2))
error ("besself: only one filter band allowed");
elseif (length(W)==2 && !(W(1) < W(2)))
error ("besself: first band edge must be smaller than second");
endif
if ( digital && !all(W >= 0 & W <= 1))
error ("besself: critical frequencies must be in (0 1)");
elseif ( !digital && !all(W >= 0 ))
error ("besself: critical frequencies must be in (0 inf)");
endif
## Prewarp to the band edges to s plane
if digital
T = 2; # sampling frequency of 2 Hz
W = 2/T*tan(pi*W/T);
endif
## Generate splane poles for the prototype bessel filter
[zero, pole, gain] = besselap(n);
## splane frequency transform
[zero, pole, gain] = sftrans(zero, pole, gain, W, stop);
## Use bilinear transform to convert poles to the z plane
if digital
[zero, pole, gain] = bilinear(zero, pole, gain, T);
endif
## convert to the correct output form
if nargout==2,
a = real(gain*poly(zero));
b = real(poly(pole));
elseif nargout==3,
a = zero;
b = pole;
c = gain;
else
## output ss results
[a, b, c, d] = zp2ss (zero, pole, gain);
endif
endfunction
|