octave-signal 1.2.2-1build1 / usr / share / octave / packages / signal-1.2.2 / fir1.m

## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
##
## 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/>.

## usage: b = fir1(n, w [, type] [, window] [, noscale])
##
## Produce an order n FIR filter with the given frequency cutoff,
## returning the n+1 filter coefficients in b.  
##
## n: order of the filter (1 less than the length of the filter)
## w: band edges
##    strictly increasing vector in range [0, 1]
##    singleton for highpass or lowpass, vector pair for bandpass or
##    bandstop, or vector for alternating pass/stop filter.
## type: choose between pass and stop bands
##    'high' for highpass filter, cutoff at w
##    'stop' for bandstop filter, edges at w = [lo, hi]
##    'DC-0' for bandstop as first band of multiband filter
##    'DC-1' for bandpass as first band of multiband filter
## window: smoothing window
##    defaults to hamming(n+1) row vector
##    returned filter is the same shape as the smoothing window
## noscale: choose whether to normalize or not
##    'scale': set the magnitude of the center of the first passband to 1
##    'noscale': don't normalize
##
## To apply the filter, use the return vector b:
##       y=filter(b,1,x);
##
## Examples:
##   freqz(fir1(40,0.3));
##   freqz(fir1(15,[0.2, 0.5], 'stop'));  # note the zero-crossing at 0.1
##   freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));

## TODO: Consider using exact expression (in terms of sinc) for the
## TODO:    impulse response rather than relying on fir2.
## TODO: Find reference to the requirement that order be even for
## TODO:    filters that end high.  Figure out what to do with the
## TODO:    window in these cases
function b = fir1(n, w, varargin)

  if nargin < 2 || nargin > 5
    print_usage;
  endif
  
  ## Assign default window, filter type and scale.
  ## If single band edge, the first band defaults to a pass band to 
  ## create a lowpass filter.  If multiple band edges, the first band 
  ## defaults to a stop band so that the two band case defaults to a 
  ## band pass filter.  Ick.
  window  = [];
  scale   = 1;
  ftype   = (length(w)==1);

  ## sort arglist, normalize any string
  for i=1:length(varargin)
    arg = varargin{i}; 
    if ischar(arg), arg=lower(arg);end
    if isempty(arg) continue; end  # octave bug---can't switch on []
    switch arg
      case {'low','stop','dc-1'},             ftype  = 1;
      case {'high','pass','bandpass','dc-0'}, ftype  = 0;
      case {'scale'},                         scale  = 1;
      case {'noscale'},                       scale  = 0;
      otherwise                               window = arg;
    end
  endfor

  ## build response function according to fir2 requirements
  bands = length(w)+1;
  f = zeros(1,2*bands);
  f(1) = 0; f(2*bands)=1;
  f(2:2:2*bands-1) = w;
  f(3:2:2*bands-1) = w;
  m = zeros(1,2*bands);
  m(1:2:2*bands) = rem([1:bands]-(1-ftype),2);
  m(2:2:2*bands) = m(1:2:2*bands);

  ## Increment the order if the final band is a pass band.  Something
  ## about having a nyquist frequency of zero causing problems.
  if rem(n,2)==1 && m(2*bands)==1, 
    warning("n must be even for highpass and bandstop filters. Incrementing.");
    n = n+1;
    if isvector(window) && isreal(window) && !ischar(window)
      ## Extend the window using interpolation
      M = length(window);
      if M == 1,
        window = [window; window];
      elseif M < 4
        window = interp1(linspace(0,1,M),window,linspace(0,1,M+1),'linear');
      else
        window = interp1(linspace(0,1,M),window,linspace(0,1,M+1),'spline');
      endif
    endif
  endif

  ## compute the filter
  b = fir2(n, f, m, 512, 2, window);

  ## normalize filter magnitude
  if scale == 1
    ## find the middle of the first band edge
    ## find the frequency of the normalizing gain
    if m(1) == 1
      ## if the first band is a passband, use DC gain
      w_o = 0;
    elseif f(4) == 1
      ## for a highpass filter,
      ## use the gain at half the sample frequency
      w_o = 1;
    else
      ## otherwise, use the gain at the center
      ## frequency of the first passband
      w_o = f(3) + (f(4)-f(3))/2;
    endif

    ## compute |h(w_o)|^-1
    renorm = 1/abs(polyval(b, exp(-1i*pi*w_o)));

    ## normalize the filter
    b = renorm*b;
  endif
endfunction

%!demo
%! freqz(fir1(40,0.3));
%!demo
%! freqz(fir1(15,[0.2, 0.5], 'stop'));  # note the zero-crossing at 0.1
%!demo
%! freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));

%!assert(fir1(2, .5, 'low', @hanning, 'scale'), [0 1 0]);
%!assert(fir1(2, .5, 'low', "hanning", 'scale'), [0 1 0]);
%!assert(fir1(2, .5, 'low', hanning(3), 'scale'), [0 1 0]);

%!assert(fir1(10,.5,'noscale'), fir1(10,.5,'low','hamming','noscale'));
%!assert(fir1(10,.5,'high'), fir1(10,.5,'high','hamming','scale'));
%!assert(fir1(10,.5,'boxcar'), fir1(10,.5,'low','boxcar','scale'));
%!assert(fir1(10,.5,'hanning','scale'), fir1(10,.5,'scale','hanning','low'));
%!assert(fir1(10,.5,'haNNing','NOscale'), fir1(10,.5,'noscale','Hanning','LOW'));
%!assert(fir1(10,.5,'boxcar',[]), fir1(10,.5,'boxcar'));