octave_packages/signal-1.1.3/fir1.m

   1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
   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 ## usage: b = fir1(n, w [, type] [, window] [, noscale])
  17 ##
  18 ## Produce an order n FIR filter with the given frequency cutoff,
  19 ## returning the n+1 filter coefficients in b.
  20 ##
  21 ## n: order of the filter (1 less than the length of the filter)
  22 ## w: band edges
  23 ##    strictly increasing vector in range [0, 1]
  24 ##    singleton for highpass or lowpass, vector pair for bandpass or
  25 ##    bandstop, or vector for alternating pass/stop filter.
  26 ## type: choose between pass and stop bands
  27 ##    'high' for highpass filter, cutoff at w
  28 ##    'stop' for bandstop filter, edges at w = [lo, hi]
  29 ##    'DC-0' for bandstop as first band of multiband filter
  30 ##    'DC-1' for bandpass as first band of multiband filter
  31 ## window: smoothing window
  32 ##    defaults to hamming(n+1) row vector
  33 ##    returned filter is the same shape as the smoothing window
  34 ## noscale: choose whether to normalize or not
  35 ##    'scale': set the magnitude of the center of the first passband to 1
  36 ##    'noscale': don't normalize
  37 ##
  38 ## To apply the filter, use the return vector b:
  39 ##       y=filter(b,1,x);
  40 ##
  41 ## Examples:
  42 ##   freqz(fir1(40,0.3));
  43 ##   freqz(fir1(15,[0.2, 0.5], 'stop'));  # note the zero-crossing at 0.1
  44 ##   freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));
  45
  46 ## TODO: Consider using exact expression (in terms of sinc) for the
  47 ## TODO:    impulse response rather than relying on fir2.
  48 ## TODO: Find reference to the requirement that order be even for
  49 ## TODO:    filters that end high.  Figure out what to do with the
  50 ## TODO:    window in these cases
  51 function b = fir1(n, w, varargin)
  52
  53   if nargin < 2 || nargin > 5
  54     print_usage;
  55   endif
  56
  57   ## Assign default window, filter type and scale.
  58   ## If single band edge, the first band defaults to a pass band to
  59   ## create a lowpass filter.  If multiple band edges, the first band
  60   ## defaults to a stop band so that the two band case defaults to a
  61   ## band pass filter.  Ick.
  62   window  = [];
  63   scale   = 1;
  64   ftype   = (length(w)==1);
  65
  66   ## sort arglist, normalize any string
  67   for i=1:length(varargin)
  68     arg = varargin{i};
  69     if ischar(arg), arg=lower(arg);end
  70     if isempty(arg) continue; end  # octave bug---can't switch on []
  71     switch arg
  72       case {'low','stop','dc-1'},             ftype  = 1;
  73       case {'high','pass','bandpass','dc-0'}, ftype  = 0;
  74       case {'scale'},                         scale  = 1;
  75       case {'noscale'},                       scale  = 0;
  76       otherwise                               window = arg;
  77     end
  78   endfor
  79
  80   ## build response function according to fir2 requirements
  81   bands = length(w)+1;
  82   f = zeros(1,2*bands);
  83   f(1) = 0; f(2*bands)=1;
  84   f(2:2:2*bands-1) = w;
  85   f(3:2:2*bands-1) = w;
  86   m = zeros(1,2*bands);
  87   m(1:2:2*bands) = rem([1:bands]-(1-ftype),2);
  88   m(2:2:2*bands) = m(1:2:2*bands);
  89
  90   ## Increment the order if the final band is a pass band.  Something
  91   ## about having a nyquist frequency of zero causing problems.
  92   if rem(n,2)==1 && m(2*bands)==1,
  93     warning("n must be even for highpass and bandstop filters. Incrementing.");
  94     n = n+1;
  95     if isvector(window) && isreal(window) && !ischar(window)
  96       ## Extend the window using interpolation
  97       M = length(window);
  98       if M == 1,
  99         window = [window; window];
 100       elseif M < 4
 101         window = interp1(linspace(0,1,M),window,linspace(0,1,M+1),'linear');
 102       else
 103         window = interp1(linspace(0,1,M),window,linspace(0,1,M+1),'spline');
 104       endif
 105     endif
 106   endif
 107
 108   ## compute the filter
 109   b = fir2(n, f, m, 512, 2, window);
 110
 111   ## normalize filter magnitude
 112   if scale == 1
 113     ## find the middle of the first band edge
 114     ## find the frequency of the normalizing gain
 115     if m(1) == 1
 116       ## if the first band is a passband, use DC gain
 117       w_o = 0;
 118     elseif f(4) == 1
 119       ## for a highpass filter,
 120       ## use the gain at half the sample frequency
 121       w_o = 1;
 122     else
 123       ## otherwise, use the gain at the center
 124       ## frequency of the first passband
 125       w_o = f(3) + (f(4)-f(3))/2;
 126     endif
 127
 128     ## compute |h(w_o)|^-1
 129     renorm = 1/abs(polyval(b, exp(-1i*pi*w_o)));
 130
 131     ## normalize the filter
 132     b = renorm*b;
 133   endif
 134 endfunction
 135
 136 %!demo
 137 %! freqz(fir1(40,0.3));
 138 %!demo
 139 %! freqz(fir1(15,[0.2, 0.5], 'stop'));  # note the zero-crossing at 0.1
 140 %!demo
 141 %! freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));
 142
 143 %!assert(fir1(2, .5, 'low', @hanning, 'scale'), [0 1 0]');
 144 %!assert(fir1(2, .5, 'low', "hanning", 'scale'), [0 1 0]');
 145 %!assert(fir1(2, .5, 'low', hanning(3), 'scale'), [0 1 0]');
 146
 147 %!assert(fir1(10,.5,'noscale'), fir1(10,.5,'low','hamming','noscale'));
 148 %!assert(fir1(10,.5,'high'), fir1(10,.5,'high','hamming','scale'));
 149 %!assert(fir1(10,.5,'boxcar'), fir1(10,.5,'low','boxcar','scale'));
 150 %!assert(fir1(10,.5,'hanning','scale'), fir1(10,.5,'scale','hanning','low'));
 151 %!assert(fir1(10,.5,'haNNing','NOscale'), fir1(10,.5,'noscale','Hanning','LOW'));
 152 %!assert(fir1(10,.5,'boxcar',[]), fir1(10,.5,'boxcar'));