Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[6026] trunk/octave-forge/main/signal/inst
|
From: <ha...@us...> - 2009-07-18 21:06:20
|
Revision: 6026
http://octave.svn.sourceforge.net/octave/?rev=6026&view=rev
Author: hauberg
Date: 2009-07-18 21:06:08 +0000 (Sat, 18 Jul 2009)
Log Message:
-----------
New functions: besself and besselap (from Thomas Sailer)
Added Paths:
-----------
trunk/octave-forge/main/signal/inst/besselap.m
trunk/octave-forge/main/signal/inst/besself.m
Added: trunk/octave-forge/main/signal/inst/besselap.m
===================================================================
--- trunk/octave-forge/main/signal/inst/besselap.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/besselap.m 2009-07-18 21:06:08 UTC (rev 6026)
@@ -0,0 +1,52 @@
+## Copyright (C) 2009 Thomas Sailer
+##
+## 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 2 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/>.
+
+## Return bessel analog filter prototype.
+##
+## References:
+##
+## http://en.wikipedia.org/wiki/Bessel_polynomials
+
+function [zero, pole, gain] = besselap (n)
+
+ if (nargin>1 || nargin<1)
+ usage ("[z, p, g] = besselap (n)");
+ end
+
+ ## interpret the input parameters
+ if (!(length(n)==1 && n == round(n) && n > 0))
+ error ("besselap: filter order n must be a positive integer");
+ end
+
+ p0=1;
+ p1=[1 1];
+ for nn=2:n;
+ px=(2*nn-1)*p1;
+ py=[p0 0 0];
+ px=prepad(px,max(length(px),length(py)),0);
+ py=prepad(py,length(px));
+ p0=p1;
+ p1=px+py;
+ endfor;
+ % p1 now contains the reverse bessel polynomial for n
+
+ % scale it by replacing s->s/w0 so that the gain becomes 1
+ p1=p1.*p1(length(p1)).^((length(p1)-1:-1:0)/(length(p1)-1));
+
+ zero=[];
+ pole=roots(p1);
+ gain=1;
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/besself.m
===================================================================
--- trunk/octave-forge/main/signal/inst/besself.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/besself.m 2009-07-18 21:06:08 UTC (rev 6026)
@@ -0,0 +1,119 @@
+## Copyright (C) 2009 Thomas Sailer
+## Copyright (C) 1999 Paul Kienzle
+##
+## 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 2 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/>.
+
+## 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.
+
+## Author: Paul Kienzle <pki...@us...>
+## Modified by: Doug Stewart <da...@sy...> Feb, 2003
+
+function [a, b, c, d] = besself (n, W, varargin)
+
+ if (nargin>4 || nargin<2) || (nargout>4 || nargout<2)
+ usage ("[b, a] or [z, p, g] or [a,b,c,d] = besself (n, W [, 'ftype'][,'z'])");
+ end
+
+ ## interpret the input parameters
+ if (!(length(n)==1 && n == round(n) && n > 0))
+ error ("besself: filter order n must be a positive integer");
+ end
+
+ 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
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
