octave_packages/signal-1.1.3/pei_tseng_notch.m

   1 ## Copyright (C) 2011 Alexander Klein <alexander.klein@math.uni-giessen.de>
   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 ## -*- texinfo -*-
  17 ## @deftypefn{Function File} { [ @var{b}, @var{a} ] } = pei_tseng_notch ( @var{frequencies}, @var{bandwidths}
  18 ## Return coefficients for an IIR notch-filter with one or more filter frequencies and according (very narrow) bandwidths
  19 ## to be used with @code{filter} or @code{filtfilt}.
  20 ## The filter construction is based on an allpass which performs a reversal of phase at the filter frequencies.
  21 ## Thus, the mean of the phase-distorted and the original signal has the respective frequencies removed.
  22 ## See the demo for an illustration.
  23 ##
  24 ## Original source:
  25 ## Pei, Soo-Chang, and Chien-Cheng Tseng
  26 ## "IIR Multiple Notch Filter Design Based on Allpass Filter"
  27 ## 1996 IEEE Tencon
  28 ## doi: 10.1109/TENCON.1996.608814)
  29 ## @end deftypefn
  30
  31 ## TODO: Implement Laplace-space frequencies and bandwidths, and perhaps better range checking for bandwidths?
  32
  33 function [ b, a ] = pei_tseng_notch ( frequencies, bandwidths )
  34
  35   err = nargchk ( 2, 2, nargin, "string" );
  36
  37   if ( err )
  38     error ( err );
  39   elseif ( !isvector ( frequencies ) || !isvector ( bandwidths ) )
  40     error ( "All arguments must be vectors!" )
  41   elseif ( length ( frequencies ) != length ( bandwidths ) )
  42     error ( "All arguments must be of equal length!" )
  43   elseif ( !all ( frequencies > 0 && frequencies < 1 ) )
  44     error ( "All frequencies must be in (0, 1)!" )
  45   elseif ( !all ( bandwidths > 0 && bandwidths < 1 ) )
  46     error ( "All bandwidths must be in (0, 1)!" )
  47   endif
  48
  49   ##Ensure row vectors
  50   frequencies = frequencies (:)';
  51   bandwidths  = bandwidths (:)';
  52
  53
  54   ##Normalise appropriately
  55   frequencies *= pi;
  56   bandwidths  *= pi;
  57   M2           = 2 * length ( frequencies );
  58
  59
  60   ##Splice centre and offset frequencies ( Equation 11 )
  61   omega = vec ( [ frequencies - bandwidths / 2; frequencies ] );
  62
  63   ##Splice centre and offset phases ( Equations 12 )
  64   factors = ( 1 : 2 : M2 );
  65   phi     = vec ( [ -pi * factors + pi / 2; -pi * factors ] );
  66
  67   ##Create linear equation
  68   t_beta = tan ( ( phi + M2 * omega ) / 2 );
  69
  70   Q = zeros ( M2 );
  71
  72   for k = 1 : M2
  73     Q ( : ,k ) = sin ( k .* omega ) - t_beta .* cos ( k .* omega );
  74   endfor
  75
  76   ##Compute coefficients of system function ( Equations 19, 20 ) ...
  77   h_a   = ( Q \ t_beta )' ;
  78   denom = [ 1, h_a ];
  79   num   = [ fliplr( h_a ), 1 ];
  80
  81   ##... and transform them to coefficients for difference equations
  82   a = denom;
  83   b = ( num + denom ) / 2;
  84
  85 endfunction
  86
  87 %!test
  88 %! ##2Hz bandwidth
  89 %! sf = 800; sf2 = sf/2;
  90 %! data=[sinetone(49,sf,10,1),sinetone(50,sf,10,1),sinetone(51,sf,10,1)];
  91 %! [b, a] = pei_tseng_notch ( 50 / sf2, 2 / sf2 );
  92 %! filtered = filter ( b, a, data );
  93 %! damp_db = 20 * log10 ( max ( filtered ( end - 1000 : end, : ) ) );
  94 %! assert ( damp_db, [ -3 -251.9 -3 ], 0.1 )
  95
  96 %!test
  97 %! ##1Hz bandwidth
  98 %! sf = 800; sf2 = sf/2;
  99 %! data=[sinetone(49.5,sf,10,1),sinetone(50,sf,10,1),sinetone(50.5,sf,10,1)];
 100 %! [b, a] = pei_tseng_notch ( 50 / sf2, 1 / sf2 );
 101 %! filtered = filter ( b, a, data );
 102 %! damp_db = 20 * log10 ( max ( filtered ( end - 1000 : end, : ) ) );
 103 %! assert ( damp_db, [ -3 -240.4 -3 ], 0.1 )
 104
 105 %!demo
 106 %! sf = 800; sf2 = sf/2;
 107 %! data=[[1;zeros(sf-1,1)],sinetone(49,sf,1,1),sinetone(50,sf,1,1),sinetone(51,sf,1,1)];
 108 %! [b,a]=pei_tseng_notch ( 50 / sf2, 2/sf2 );
 109 %! filtered = filter(b,a,data);
 110 %!
 111 %! clf
 112 %! subplot ( columns ( filtered ), 1, 1)
 113 %! plot(filtered(:,1),";Impulse response;")
 114 %! subplot ( columns ( filtered ), 1, 2 )
 115 %! plot(filtered(:,2),";49Hz response;")
 116 %! subplot ( columns ( filtered ), 1, 3 )
 117 %! plot(filtered(:,3),";50Hz response;")
 118 %! subplot ( columns ( filtered ), 1, 4 )
 119 %! plot(filtered(:,4),";51Hz response;")