## Copyright (C) 2014 James-Adam Renquinha Henri ## ## 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 . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{b}, @var{a}] =} legenp (@var{n}, @var{Wc}) ## low pass filter with cutoff pi*@var{Wc} radians ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} legenp (@var{n}, @var{Wc}, "high") ## high pass filter with cutoff pi*@var{Wc} radians ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} legenp (@var{n}, [@var{Wl}, @var{Wh}]) ## band pass filter with edges pi*@var{Wl} and pi*@var{Wh} radians ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} legenp (@var{n}, [@var{Wl}, @var{Wh}], "stop") ## band reject filter with edges pi*@var{Wl} and pi*@var{Wh} radians ## @deftypefnx {Function File} {[@var{z}, @var{p}, @var{g}] =} legenp (@dots{}) ## return filter as zero-pole-gain rather than coefficients of the numerator and denominator polynomials. ## @deftypefnx {Function File} {[@var{a}, @var{b}, @var{c}, @var{d}] =} legenp (@dots{}) ## return state-space matrices ## @deftypefnx {Function File} {[@dots{}] =} legenp (@dots{}, "s") ## return a Laplace space filter, @var{W} can be larger than 1. ## ## Generate a Legendre-Papoulis filter (also known as Optimal "L" filter). ## Default is a discrete space (Z) filter. ## ## Compared to Butterworth filters, they have a much sharper rolloff at the ## expense of a somewhat less flat passband; but it has monotonic frequency ## response unlike the Chebyshev type I. ## ## ## References: ## ## Notes on "L" (Optimal) Filters, by C. Bond, 2011 ## ## @seealso{butter, cheby1} ## @end deftypefn function [a, b, c, d] = legenp (n, W, varargin) # Please note that this code is entirely based on butter. # 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 ("legenp: filter order n must be a positive integer"); endif stop = 0; digital = 1; 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 ("legenp: expected [high|stop] or [s|z]"); endswitch endfor [r, c]=size (W); if (!(length (W) <= 2 && (r == 1 || c == 1))) error ("legenp: frequency must be given as w0 or [w0, w1]"); elseif (!(length (W) == 1 || length (W) == 2)) error ("legenp: only one filter band allowed"); elseif (length (W) == 2 && !(W(1) < W(2))) error ("legenp: first band edge must be smaller than second"); endif if (digital && !all (W >= 0 & W <= 1)) error ("legenp: critical frequencies must be in (0 1)"); elseif ( !digital && !all (W >= 0 )) error ("legenp: 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 Legendre-Papoulis filter C = 1; # default cutoff frequency # HERE: [zero, pole, gain] = legenpap (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 %!shared sf, sf2, off_db %! off_db = 0.5; %! ##Sampling frequency must be that high to make the low pass filters pass. %! sf = 6000; sf2 = sf/2; %! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)]; %!test %! ##Test low pass order 1 with 3dB @ 50Hz %! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)]; %! [b, a] = legenp ( 1, 50 / sf2 ); %! filtered = filter ( b, a, data ); %! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) ); %! assert ( [ damp_db( 4 ) - damp_db( 5 ), damp_db( 1 : 3 ) ], [ 6 0 0 -3 ], off_db ) %! ## Frequency response is monotonic %! H = freqz(b,a,2^16); %! assert(max(diff(abs(H))) < 1e-6); %!test %! ##Test low pass order 4 with 3dB @ 20Hz %! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(20,sf,10,1),sinetone(100,sf,10,1),sinetone(200,sf,10,1)]; %! [b, a] = legenp ( 4, 20 / sf2 ); %! filtered = filter ( b, a, data ); %! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) ); %! assert ( [ damp_db( 4 ) - damp_db( 5 ), damp_db( 1 : 3 ) ], [ 24 0 0 -3 ], off_db ) %! ## Frequency response is monotonic %! H = freqz(b,a,2^16); %! assert(max(diff(abs(H))) < 1e-6); %!test %! ##Test high pass order 1 with 3dB @ 50Hz %! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)]; %! [b, a] = legenp ( 1, 50 / sf2, "high" ); %! filtered = filter ( b, a, data ); %! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) ); %! assert ( [ damp_db( 2 ) - damp_db( 1 ), damp_db( 3 : end ) ], [ 6 -3 0 0 ], off_db ) %! ## Frequency response is monotonic (will always increase for an hi-pass) %! H = freqz(b,a,2^16); %! assert(min(diff(abs(H))) > -1e-6); %!test %! ##Test high pass order 4 with 3dB @ 50Hz %! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)]; %! [b, a] = legenp ( 4, 50 / sf2, "high" ); %! filtered = filter ( b, a, data ); %! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) ); %! assert ( [ damp_db( 2 ) - damp_db( 1 ), damp_db( 3 : end ) ], [ 24 -3 0 0 ], off_db ) %! ## Frequency response is monotonic (will always increase for an hi-pass) %! H = freqz(b,a,2^16); %! assert(min(diff(abs(H))) > -1e-6); %!demo %! sf = 800; sf2 = sf/2; %! data=[[1;zeros(sf-1,1)],sinetone(25,sf,1,1),sinetone(50,sf,1,1),sinetone(100,sf,1,1)]; %! [b,a]=legenp ( 1, 50 / sf2 ); %! filtered = filter(b,a,data); %! %! clf %! subplot ( columns ( filtered ), 1, 1) %! plot(filtered(:,1),";Impulse response;") %! subplot ( columns ( filtered ), 1, 2 ) %! plot(filtered(:,2),";25Hz response;") %! subplot ( columns ( filtered ), 1, 3 ) %! plot(filtered(:,3),";50Hz response;") %! subplot ( columns ( filtered ), 1, 4 ) %! plot(filtered(:,4),";100Hz response;")