## Copyright (C) 2018 Ionescu Vlad <vlad.inf@gmail.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/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{h} =} fir (@var{n}, @var{f})
## @deftypefnx {Function File} {@var{h} =} fir (@var{n}, @var{f}, @var{a})
## @deftypefnx {Function File} {@var{h} =} fir (@var{n}, @var{f}, @var{a}, @var{win})
## @deftypefnx {Function File} {@var{h} =} fir (@var{n}, @var{f}, @var{a}, @var{win}, @var{scale})
##
## Create a single, or multiband FIR filter of order @var{n} with the frequency
## band edges @var{f}, thus @var{h} will have a length of @var{n}+1. For types
## I & II, specifying either a band edge at DC, or one at Nyquist, but not both,
## will mean the difference between having a DC, or not. For types III & IV, the
## same difference will be by either specifying both DC and Nyquist, or none.
##
## @var{a} specifies the gain of each passband. If @var{f} starts at 0, the
## first passband is between 0 and f(1), the 2nd between f(2) and f(3), if they
## exist, and so on. Similar for @var{f} with 1 at the end.
##
## An optional shaping @var{win} can be given as a vector with length @var{n}+1.
## Defaults to a rectangular (boxcar) window of length @var{n}+1.
##
## @var{scale} is a flag with the value "scale", or "y", meaning the filter's
## coefficients will be normalized such that the magnitude response at DC, or
## the Nyquist, or the center of the first passband is 1, depending on the case.
##
## To apply the filter, use the return vector @var{h} with the @code{filter}
## function, for example @code{y = filter (h, 1, x)}.
##
## Examples:
## @example
## freqz (fir (41, [0 0.37]));                  # type II lowpass
## freqz (fir (16, [0.41 1], "scale"));         # type I highpass, normalized
## freqz (fir (24, [0 0.3 0.6], [0.5 1]);       # type I bandstop
## freqz (fir (33, [0.22 0.55 1], 0.8, "y"));   # type II bandstop
## freqz (fir (80, [0 0.12 0.34 0.56 0.78 0.9], [1 -2 0.1], bartlett (81))); # multiband with DC
## freqz (fir (80, [0.12 0.34 0.56 0.78 0.9 1], [-1 2 -0.1], hanning (81))); # multiband without DC
## freqz (fir (51, [0 1], 0.7, hamming (52)));       # type III Hilbert transformer
## Unorthodox examples. The purpose is not to be different, but to be able to
## exemplify clearly what happens if you *do* go out of the recommended way.
## There's a difference between saying "Don't do that!" and "If you do it, this
## is what will happen.".
## freqz (fir (21, [0.35 1], "y"));     # any type II that ends in Nyquist
## freqz (fir (38, [0 0.23 0.45 1]));   # any type III with DC and/or Nyquist
## freqz (fir (55, [0 0.28 1]));        # any type IV starting with DC
## freqz (fir (44, [0 0.4 0.6 1]));     # multiband Hilbert transformer
## freqz (fir (99, [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1], [0.2 -0.4 0.6 -0.8 1], kaiser(100, 5), "y"));
## @end example
## @seealso{filter, fir1, fir2}
## @end deftypefn

function h = fir (N, F, varargin)

  ## First check the required number of arguments
  if (nargin < 2 || nargin > 5)
    print_usage;
  endif

  ## Order must be a one one element vector...
  if (length (N) != 1)
    error ("The order (N) must be a vector of size 1.")
  endif
  ## ...and proper valued
  if ((N <= 0) || ischar (N))
    error ("The order (N) must be positive definite.")
  endif

  ## The frequencies must be positive values, ...
  if (min (F) < 0)
    error ("The frequencies must be positive.")
  endif
  ## ... strictly increasing values, ...
  if (min (diff (F)) <= 0)
    error ("The frequencies must be strictly increasing.")
  endif
  ## ... and their range between 0 and 1, closed interval at both ends.
  if ((F(1) < 0) || (F(end) > 1))
    error ("The frequencies must lie in the interval [0..1], with 1 being Nyquist.")
  endif
  ## Only real numbers
  if (iscomplex (F) || ischar (F))
    error ("The frequency vector must be made of real numbers, only.")
  endif

  ## Default values
  N = fix (N);  # silently consider the integer part of N
  lenF = length(F);
  fType = (F(1) == 0 && F(end) == 1) || (F(1) > 0 && F(end) < 1);
  scale = 0;
  win = ones (1, N+1);
  DC = (F(1) == 0);
  ## A has different lengths
  if (! fType)  # types I, II
    if (DC)
      A = ones (1, ceil (lenF/2));
    else
      A = ones (1, floor (lenF/2));
    endif
  else          # types III, IV
    if (DC)
      A = ones (1, floor (lenF/2));
    else
      A = ones (1, lenF/2);
    endif
  endif

  ## Out of the 3 possible variable arguments, 1 will be string
  for k=1:length (varargin)
    if (isempty (varargin))
      continue; # as per fir1.m, line 78, this may be to circumvent a bug
    endif
    switch (varargin{k})
      case {"scale" "y"}
        scale = 1;
      otherwise
        if (ischar (varargin{k}))
          error ("The flag must be \"scale\", or \"y\".")
        endif
        if (length (varargin{k}) == N+1 )
          win = varargin{k};
        else
          A = varargin{k};
        endif
    endswitch
  endfor

  ## Some helpers
  oddN = mod (N+1, 2);
  lenA = length (A);
  oddF = mod(lenF, 2);
  midN = ceil ((N + 1)/2);

  ## Handle errors from varargin
  ## A needs to have different lengths, depending on the type of FIR
  if (! fType) # types I & II
    if (DC && lenA != ceil (lenF/2))
      error ("For type I & II with DC, the length of the amplitude vector must be ceil(length(F)/2).")
    elseif (! DC && lenA != floor (lenF/2))
      error ("For type I & II without DC, the length of the amplitude vector must be floor(length(F)/2).")
    endif
  else  # types III & IV
    if (DC && lenA != floor (lenF/2))
      error ("For type III & IV with DC, the length of the amplitude vector must be floor(length(F)/2).")
    elseif (!DC && lenA != lenF/2)
      error ("For type III & IV without DC, the length of the amplitude vector must be length(F)/2.")
    endif
  endif
  ## Only real numbers
  if (iscomplex (A) || ischar (A))
    error ("The values of the amplitude vector must be real numbers, only.")
  endif

  ## window
  if (length (win) != N+1)
    error ("The length of the window must be the same as the filter's.")
  endif
  ## Only real numbers
  if (iscomplex (win) || ischar (win))
    error ("The values of the window vector must be real numbers, only.")
  endif

  ## Make sure the vectors are columns
  F = F(:)';
  ## F needs some adjusting since A must have one value per band, but F can
  ## be, for ex., [0 0.3] (types I,II) and [0 0.3 1] (types III,IV), both being
  ## lowpass, but different lengths. A needs the same treatment.
  if ( (DC && ! fType && oddF) || (! DC && fType && oddF) )
    F = [F 1];
    lenF += 1;
    oddF = 1 - oddF;
  endif
  win = win(:)';
  A = [A(:)' ; A(:)'](:)';
  ## The lengths of A and F are interdependent.
  if ( (! DC && ! fType && oddF) || (DC && fType && oddF) )
    A = [A 0];
  endif

  ## Some more helpers
  W = F*pi;
  endF = lenF - (! DC && oddN);

  ## Begin carnage
  n = [-N/2:1:N/2];
  h = zeros (1, N+1);
  ## Types I & II
  if (! fType)
    ## sin(0) or sin(pi) are useless, just account for the sign
    for k=1:endF
      h += (-1)^k*A(k)*sin(W(k)*n)./(pi*n).*win;
    endfor
    ## Even N will result in h(midN)=inf. Here, all the frequencies are needed.
    if (oddN)
      h(midN) = 0;
      for k=1:lenF
        h(midN) += (-1)^k*A(k)*F(k);
      endfor
    endif
  ## Types III & IV
  else
    ## For types III & IV, cos(0) and cos(pi) matter.
    for k=1:lenF
      h += (-1)^(k + 1)*A(k)*cos(W(k)*n)./(pi*n).*win;
    endfor
    if (oddN)
      h(midN) = 0; # midpoint is null for type III
    endif
  endif

  ## Normalizing for unity gain
  if (scale)
    ## Not valid for types III & IV, one, because there can be no DC, and two,
    ## their sum is 0, so they are considered, by all practical means, bandpass
    ## and/or bandstop variants, which means the normalizing is done at the
    ## center of the first bandwidth, whether F(1)=0 or not. The same is valid
    ## for type II with Nyquist (so to speak).
    if (DC && ! fType)
      h = A(1)*h/sum (h);
    elseif (! DC && ! fType && oddN)
      h = A(end)*h/abs (polyval (h, -1));
    else
      h = A(1)*h/abs (polyval (h, exp(-0.5i*(W(1) + W(2))) ));
    endif
  endif

endfunction

%% expected output
%!test
%! x = [-0.01497698351644062; ...
%!      -0.01020620096284455; ...
%!      0.008051172211853279; ...
%!      0.01813306531629943; ...
%!      0.006248855546642409; ...
%!      -0.0151891541716538; ...
%!      -0.02000749470712845; ...
%!      3.703551396872747*10^-4; ...
%!      0.02352639910729612; ...
%!      0.01987710435709456; ...
%!      -0.0107156702183056; ...
%!      -0.03346910762203006; ...
%!      -0.01647493965701502; ...
%!      0.02756343124616879; ...
%!      0.0468059024239567; ...
%!      0.006350821104335874; ...
%!      -0.06144311077699433; ...
%!      -0.07272786646455694; ...
%!      0.02972318687964276; ...
%!      0.2090466916579446; ...
%!      0.3495187814185787; ...
%!      0.3495187814185787; ...
%!      0.2090466916579446; ...
%!      0.02972318687964276; ...
%!      -0.07272786646455694; ...
%!      -0.06144311077699433; ...
%!      0.006350821104335874; ...
%!      0.0468059024239567; ...
%!      0.02756343124616879; ...
%!      -0.01647493965701502; ...
%!      -0.03346910762203006; ...
%!      -0.0107156702183056; ...
%!      0.01987710435709456; ...
%!      0.02352639910729612; ...
%!      3.703551396872747*10^-4; ...
%!      -0.02000749470712845; ...
%!      -0.0151891541716538; ...
%!      0.006248855546642409; ...
%!      0.01813306531629943; ...
%!      0.008051172211853279; ...
%!      -0.01020620096284455; ...
%!      -0.01497698351644062];
%! N = 41;
%! f = [0 0.37];
%! h = fir (N, f)';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [0.02929434394789889; ...
%!      -0.01725630770401515; ...
%!      -0.05029262346635695; ...
%!      -0.009516036365963014; ...
%!      0.06880171142241152; ...
%!      0.06704690212228634; ...
%!      -0.08148695470602121; ...
%!      -0.2920772488516647; ...
%!      0.5637616640054226; ...
%!      -0.2920772488516647; ...
%!      -0.08148695470602121; ...
%!      0.06704690212228634; ...
%!      0.06880171142241152; ...
%!      -0.009516036365963014; ...
%!      -0.05029262346635695; ...
%!      -0.01725630770401515; ...
%!      0.02929434394789889];
%! N = 16;
%! f = [0.41 1];
%! h = fir (N, f, "scale")';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [0.002977709252466263; ...
%!      -0.03922634046876217; ...
%!      2.923628874389531*10^-17; ...
%!      0.04794330501737597; ...
%!      -0.004466563878699571; ...
%!      -0.01970232494336491; ...
%!      0.03486362717956649; ...
%!      -0.03183098861837909; ...
%!      -0.09906990495878168; ...
%!      0.07875981297063382; ...
%!      0.1692316012429296; ...
%!      -0.1739716377561419; ...
%!      0.55; ...
%!      -0.1739716377561419; ...
%!      0.1692316012429296; ...
%!      0.07875981297063382; ...
%!      -0.09906990495878168; ...
%!      -0.03183098861837909; ...
%!      0.03486362717956649; ...
%!      -0.01970232494336491; ...
%!      -0.004466563878699571; ...
%!      0.04794330501737597; ...
%!      2.923628874389531*10^-17; ...
%!      -0.03922634046876217; ...
%!      0.002977709252466263];
%! N = 24;
%! f = [0 0.3 0.6];
%! A = [0.5 1];
%! h = fir (N, f, A)';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [0.01003380588433263; ...
%!      0.03054943447655625; ...
%!      0.008069327033012825; ...
%!      -0.01911240979492978; ...
%!      -0.006279127634478077; ...
%!      -0.003006753731602169; ...
%!      -0.03402121285782708; ...
%!      -0.02364433277622251; ...
%!      0.03557252267223805; ...
%!      0.04108655060846322; ...
%!      1.32015720630261*10^-4; ...
%!      0.0235093138332214; ...
%!      0.05190867492026852; ...
%!      -0.06184924359196657; ...
%!      -0.1849894389661927; ...
%!      -0.05455510212580623; ...
%!      0.2040318218295396; ...
%!      0.2040318218295396; ...
%!      -0.05455510212580623; ...
%!      -0.1849894389661927; ...
%!      -0.06184924359196657; ...
%!      0.05190867492026852; ...
%!      0.0235093138332214; ...
%!      1.32015720630261*10^-4; ...
%!      0.04108655060846322; ...
%!      0.03557252267223805; ...
%!      -0.02364433277622251; ...
%!      -0.03402121285782708; ...
%!      -0.003006753731602169; ...
%!      -0.006279127634478077; ...
%!      -0.01911240979492978; ...
%!      0.008069327033012825; ...
%!      0.03054943447655625; ...
%!      0.01003380588433263];
%! N = 33;
%! f = [0.22 0.55 1];
%! A = 0.8;
%! h = fir (N, f, A, "y")';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [0.0; ...
%!      4.532610419962665*10^-5; ...
%!      0.00132766833268789; ...
%!      8.097206407907734*10^-4; ...
%!      0.001167270162370672; ...
%!      0.00210626333276858; ...
%!      -0.002064214102702842; ...
%!      -0.005746380022157374; ...
%!      0.001587810100392391; ...
%!      0.006789748295554875; ...
%!      -0.00227048018592197; ...
%!      -0.00990281544117809; ...
%!      -0.004170409103326456; ...
%!      -0.003526181767557801; ...
%!      -0.006604587816182529; ...
%!      0.01050422624406504; ...
%!      0.01657359695603568; ...
%!      -0.006021304773640971; ...
%!      -0.0181685269284046; ...
%!      0.01046490805824807; ...
%!      0.02703487969081995; ...
%!      0.007569361930096066; ...
%!      0.00745866402394203; ...
%!      0.00902796074637719; ...
%!      -0.02946707469549715; ...
%!      -0.04147353109022093; ...
%!      0.02303271450092412; ...
%!      0.03921057487778715; ...
%!      -0.03484620320766304; ...
%!      -0.0711946676779947; ...
%!      -0.01262910533114663; ...
%!      -0.009198947246326734; ...
%!      -0.01454443729490718; ...
%!      0.113250549558954; ...
%!      0.1195572497537248; ...
%!      -0.095345938066831; ...
%!      -0.1603576217899679; ...
%!      0.2415531631760952; ...
%!      0.4761051923170519; ...
%!      0.03827346739107588; ...
%!      -0.3080000000000001; ...
%!      0.03827346739107588; ...
%!      0.4761051923170519; ...
%!      0.2415531631760952; ...
%!      -0.1603576217899679; ...
%!      -0.095345938066831; ...
%!      0.1195572497537248; ...
%!      0.113250549558954; ...
%!      -0.01454443729490718; ...
%!      -0.009198947246326734; ...
%!      -0.01262910533114663; ...
%!      -0.0711946676779947; ...
%!      -0.03484620320766304; ...
%!      0.03921057487778715; ...
%!      0.02303271450092412; ...
%!      -0.04147353109022093; ...
%!      -0.02946707469549715; ...
%!      0.00902796074637719; ...
%!      0.00745866402394203; ...
%!      0.007569361930096066; ...
%!      0.02703487969081995; ...
%!      0.01046490805824807; ...
%!      -0.0181685269284046; ...
%!      -0.006021304773640971; ...
%!      0.01657359695603568; ...
%!      0.01050422624406504; ...
%!      -0.006604587816182529; ...
%!      -0.003526181767557801; ...
%!      -0.004170409103326456; ...
%!      -0.00990281544117809; ...
%!      -0.00227048018592197; ...
%!      0.006789748295554875; ...
%!      0.001587810100392391; ...
%!      -0.005746380022157374; ...
%!      -0.002064214102702842; ...
%!      0.00210626333276858; ...
%!      0.001167270162370672; ...
%!      8.097206407907734*10^-4; ...
%!      0.00132766833268789; ...
%!      4.532610419962665*10^-5; ...
%!      0.0];
%! N = 80;
%! f = [0 0.12 0.34 0.56 0.78 0.9];
%! A = [1 -2 0.1];
%! h = fir (N, f, A, bartlett (81))';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [0.0; ...
%!      5.589398138030331*10^-5; ...
%!      2.700696169574606*10^-5; ...
%!      -8.992997725840461*10^-5; ...
%!      -4.570786490232379*10^-5; ...
%!      3.741377529790033*10^-4; ...
%!      0.00182284924480958; ...
%!      -0.002148300213245028; ...
%!      -4.089525659678325*10^-5; ...
%!      0.001338057355777996; ...
%!      -0.007173373007655817; ...
%!      -9.789965904564333*10^-5; ...
%!      0.00177972635233203; ...
%!      8.265849130635226*10^-4; ...
%!      -0.004574060517063696; ...
%!      -0.01139688205022777; ...
%!      0.01510697742131124; ...
%!      0.001009233660211206; ...
%!      -0.004931453084429298; ...
%!      0.03270939591637084; ...
%!      -0.002890820867463352; ...
%!      -0.008948006659384581; ...
%!      -7.597847147421062*10^-4; ...
%!      0.0140877258884251; ...
%!      0.03351731045599608; ...
%!      -0.05719992113622888; ...
%!      4.482364821912782*10^-4; ...
%!      0.006725364027416759; ...
%!      -0.09524575706917161; ...
%!      0.01974606409253418; ...
%!      0.02960963943625397; ...
%!      0.002947420328344872; ...
%!      -0.04789873216001209; ...
%!      -0.07809370770080949; ...
%!      0.2030690054755283; ...
%!      0.004070495321897486; ...
%!      -0.02310233592934926; ...
%!      0.4698001262597655; ...
%!      -0.2288613011593811; ...
%!      -0.3708972906055919; ...
%!      0.2099999999999997; ...
%!      -0.3708972906055919; ...
%!      -0.2288613011593811; ...
%!      0.4698001262597655; ...
%!      -0.02310233592934926; ...
%!      0.004070495321897486; ...
%!      0.2030690054755283; ...
%!      -0.07809370770080949; ...
%!      -0.04789873216001209; ...
%!      0.002947420328344872; ...
%!      0.02960963943625397; ...
%!      0.01974606409253418; ...
%!      -0.09524575706917161; ...
%!      0.006725364027416759; ...
%!      4.482364821912782*10^-4; ...
%!      -0.05719992113622888; ...
%!      0.03351731045599608; ...
%!      0.0140877258884251; ...
%!      -7.597847147421062*10^-4; ...
%!      -0.008948006659384581; ...
%!      -0.002890820867463352; ...
%!      0.03270939591637084; ...
%!      -0.004931453084429298; ...
%!      0.001009233660211206; ...
%!      0.01510697742131124; ...
%!      -0.01139688205022777; ...
%!      -0.004574060517063696; ...
%!      8.265849130635226*10^-4; ...
%!      0.00177972635233203; ...
%!      -9.789965904564333*10^-5; ...
%!      -0.007173373007655817; ...
%!      0.001338057355777996; ...
%!      -4.089525659678325*10^-5; ...
%!      -0.002148300213245028; ...
%!      0.00182284924480958; ...
%!      3.741377529790033*10^-4; ...
%!      -4.570786490232379*10^-5; ...
%!      -8.992997725840461*10^-5; ...
%!      2.700696169574606*10^-5; ...
%!      5.589398138030331*10^-5; ...
%!      0.0];
%! N = 80;
%! f = [0.12 0.34 0.56 0.78 0.9 1];
%! A = [-1 2 -0.1];
%! h = fir (N, f, A, hanning (81))';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [-6.990334755408738*10^-4; ...
%!      -7.592742622395885*10^-4; ...
%!      -8.902570351910555*10^-4; ...
%!      -0.001099851794286741; ...
%!      -0.001396330324753902; ...
%!      -0.001788507076912539; ...
%!      -0.002285940466842752; ...
%!      -0.002899216466779759; ...
%!      -0.003640346078750408; ...
%!      -0.004523323335302385; ...
%!      -0.005564914287372944; ...
%!      -0.006785786070999286; ...
%!      -0.008212149570055808; ...
%!      -0.009878200050202643; ...
%!      -0.01182983780902589; ...
%!      -0.01413051804944165; ...
%!      -0.01687079425261369; ...
%!      -0.02018459160624273; ...
%!      -0.02427848614769582; ...
%!      -0.02948798939789099; ...
%!      -0.03639520195584896; ...
%!      -0.0461032907383726; ...
%!      -0.06098142761713091; ...
%!      -0.08719746951979813; ...
%!      -0.1473811565618083; ...
%!      -0.4452450390435943; ...
%!      0.4452450390435943; ...
%!      0.1473811565618083; ...
%!      0.08719746951979813; ...
%!      0.06098142761713091; ...
%!      0.0461032907383726; ...
%!      0.03639520195584896; ...
%!      0.02948798939789099; ...
%!      0.02427848614769582; ...
%!      0.02018459160624273; ...
%!      0.01687079425261369; ...
%!      0.01413051804944165; ...
%!      0.01182983780902589; ...
%!      0.009878200050202643; ...
%!      0.008212149570055808; ...
%!      0.006785786070999286; ...
%!      0.005564914287372944; ...
%!      0.004523323335302385; ...
%!      0.003640346078750408; ...
%!      0.002899216466779759; ...
%!      0.002285940466842752; ...
%!      0.001788507076912539; ...
%!      0.001396330324753902; ...
%!      0.001099851794286741; ...
%!      8.902570351910555*10^-4; ...
%!      7.592742622395885*10^-4; ...
%!      6.990334755408738*10^-4];
%! N = 51;
%! f = [0 1];
%! A = 0.7;
%! h = fir (N, f, A, hamming (52))';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [0.0574916731625867; ...
%!      -0.005054272988817951; ...
%!      0.03532635778447862; ...
%!      -0.0835833517925742; ...
%!      0.01201063075910171; ...
%!      -0.04541335600445939; ...
%!      0.1428167115246368; ...
%!      -0.03263529305854349; ...
%!      0.08045834518442294; ...
%!      -0.4337824723054122; ...
%!      0.3111772575077162; ...
%!      0.3111772575077162; ...
%!      -0.4337824723054122; ...
%!      0.08045834518442294; ...
%!      -0.03263529305854349; ...
%!      0.1428167115246368; ...
%!      -0.04541335600445939; ...
%!      0.01201063075910171; ...
%!      -0.0835833517925742; ...
%!      0.03532635778447862; ...
%!      -0.005054272988817951; ...
%!      0.0574916731625867];
%! N = 21;
%! f = [0.35 1];
%! h = fir (N, f, "y")';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [-0.02423205451259543; ...
%!      -8.175164094327593*10^-4; ...
%!      -0.02796814487513872; ...
%!      0.02675481709526449; ...
%!      -0.03075568867269888; ...
%!      -0.03088284566919951; ...
%!      -0.09526067516468953; ...
%!      -0.01113956304314307; ...
%!      -0.03201676505355007; ...
%!      0.01870978567577275; ...
%!      -0.07115186469221117; ...
%!      0.02257173940981391; ...
%!      -0.03502570831782201; ...
%!      0.01165336204853484; ...
%!      -0.2290630073420544; ...
%!      -0.1414569257946582; ...
%!      -0.2236756013656254; ...
%!      0.1713127493835815; ...
%!      -0.4476466399304864; ...
%!      0.0; ...
%!      0.4476466399304864; ...
%!      -0.1713127493835815; ...
%!      0.2236756013656254; ...
%!      0.1414569257946582; ...
%!      0.2290630073420544; ...
%!      -0.01165336204853484; ...
%!      0.03502570831782201; ...
%!      -0.02257173940981391; ...
%!      0.07115186469221117; ...
%!      -0.01870978567577275; ...
%!      0.03201676505355007; ...
%!      0.01113956304314307; ...
%!      0.09526067516468953; ...
%!      0.03088284566919951; ...
%!      0.03075568867269888; ...
%!      -0.02675481709526449; ...
%!      0.02796814487513872; ...
%!      8.175164094327593*10^-4; ...
%!      0.02423205451259543];
%! N = 38;
%! f = [0 0.23 0.45 1];
%! h = fir (N, f)';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [-0.004771346524584097; ...
%!      -0.01499888059857716; ...
%!      -0.02377746204981428; ...
%!      -0.02474797070490875; ...
%!      -0.01691363131328921; ...
%!      -0.005831645752269472; ...
%!      -2.921450462356479*10^-5; ...
%!      -0.004898135356496898; ...
%!      -0.01836947146126617; ...
%!      -0.03173339524606126; ...
%!      -0.03548803300800307; ...
%!      -0.02639318610001874; ...
%!      -0.01064276961366783; ...
%!      -3.888374623829752*10^-4; ...
%!      -0.005410955820201908; ...
%!      -0.02546479089470322; ...
%!      -0.04900624946042221; ...
%!      -0.06009348849171697; ...
%!      -0.04964808877428745; ...
%!      -0.02366261196673792; ...
%!      -0.002077225963670382; ...
%!      -0.007623378510068939; ...
%!      -0.05062092343143465; ...
%!      -0.1191573327909766; ...
%!      -0.1817119030051928; ...
%!      -0.2021630971766074; ...
%!      -0.1594329576701907; ...
%!      -0.06058898019463443; ...
%!      0.06058898019463443; ...
%!      0.1594329576701907; ...
%!      0.2021630971766074; ...
%!      0.1817119030051928; ...
%!      0.1191573327909766; ...
%!      0.05062092343143465; ...
%!      0.007623378510068939; ...
%!      0.002077225963670382; ...
%!      0.02366261196673792; ...
%!      0.04964808877428745; ...
%!      0.06009348849171697; ...
%!      0.04900624946042221; ...
%!      0.02546479089470322; ...
%!      0.005410955820201908; ...
%!      3.888374623829752*10^-4; ...
%!      0.01064276961366783; ...
%!      0.02639318610001874; ...
%!      0.03548803300800307; ...
%!      0.03173339524606126; ...
%!      0.01836947146126617; ...
%!      0.004898135356496898; ...
%!      2.921450462356479*10^-5; ...
%!      0.005831645752269472; ...
%!      0.01691363131328921; ...
%!      0.02474797070490875; ...
%!      0.02377746204981428; ...
%!      0.01499888059857716; ...
%!      0.004771346524584097];
%! N = 55;
%! f = [0 0.28 1];
%! h = fir (N, f)';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [-6.938893903907228*10^-18; ...
%!      -0.02094730684528045; ...
%!      0.0; ...
%!      -0.02315228651320469; ...
%!      -2.42861286636753*10^-17; ...
%!      -0.06774446983341569; ...
%!      3.816391647148976*10^-17; ...
%!      0.0; ...
%!      -3.816391647148976*10^-17; ...
%!      -0.08858892208985117; ...
%!      4.85722573273506*10^-17; ...
%!      -0.03999031306826262; ...
%!      0.0; ...
%!      -0.04887704930565423; ...
%!      -2.081668171172169*10^-17; ...
%!      -0.1645222838811522; ...
%!      7.632783294297951*10^-17; ...
%!      0.0; ...
%!      -4.163336342344337*10^-17; ...
%!      -0.3838853290560216; ...
%!      0.0; ...
%!      -0.4398934437508882; ...
%!      0.0; ...
%!      0.4398934437508882; ...
%!      0.0; ...
%!      0.3838853290560216; ...
%!      4.163336342344337*10^-17; ...
%!      0.0; ...
%!      -7.632783294297951*10^-17; ...
%!      0.1645222838811522; ...
%!      2.081668171172169*10^-17; ...
%!      0.04887704930565423; ...
%!      0.0; ...
%!      0.03999031306826262; ...
%!      -4.85722573273506*10^-17; ...
%!      0.08858892208985117; ...
%!      3.816391647148976*10^-17; ...
%!      0.0; ...
%!      -3.816391647148976*10^-17; ...
%!      0.06774446983341569; ...
%!      2.42861286636753*10^-17; ...
%!      0.02315228651320469; ...
%!      0.0; ...
%!      0.02094730684528045; ...
%!      6.938893903907228*10^-18];
%! N = 44;
%! f = [0 0.4 0.6 1];
%! h = fir (N, f)';
%! assert (x, h, 1e-15);

%% expected output
%!test
%! x = [-2.596982003780794*10^-4; ...
%!      3.687425121492475*10^-4; ...
%!      -5.014856683353488*10^-4; ...
%!      0.001025890224635668; ...
%!      -0.002315442483754842; ...
%!      7.163914730616564*10^-4; ...
%!      5.55274557745488*10^-4; ...
%!      -1.017156221966768*10^-4; ...
%!      9.147404918569031*10^-5; ...
%!      3.076563909696012*10^-6; ...
%!      -3.527206528143375*10^-6; ...
%!      -1.379419236884723*10^-4; ...
%!      2.021831084524683*10^-4; ...
%!      -0.001460217755938636; ...
%!      -0.002505878517562597; ...
%!      0.01085376379283025; ...
%!      -0.006509135150844266; ...
%!      0.004364459613084872; ...
%!      -0.004481128825392679; ...
%!      0.004515902634334341; ...
%!      -0.004989519282046221; ...
%!      0.006045470501789038; ...
%!      -0.007194557768823587; ...
%!      0.0131261469494572; ...
%!      -0.02682003942293659; ...
%!      0.007604249646274445; ...
%!      0.005456669361565423; ...
%!      -9.335355913526482*10^-4; ...
%!      7.901597126014116*10^-4; ...
%!      2.518658395698518*10^-5; ...
%!      -2.754155888056701*10^-5; ...
%!      -0.001033490796384705; ...
%!      0.001461842444097427; ...
%!      -0.01024620657077106; ...
%!      -0.01716116565741117; ...
%!      0.07296656815564413; ...
%!      -0.04321809837188419; ...
%!      0.02880801137929008; ...
%!      -0.0296170703853042; ...
%!      0.03013070722156086; ...
%!      -0.03392636485319551; ...
%!      0.04236597433282478; ...
%!      -0.0526889967689672; ...
%!      0.1022473973980562; ...
%!      -0.2274936863299576; ...
%!      0.0725880402231166; ...
%!      0.06160871602082477; ...
%!      -0.01355974522888309; ...
%!      0.01755403624385729; ...
%!      0.001537976526426166; ...
%!      0.001537976526426166; ...
%!      0.01755403624385729; ...
%!      -0.01355974522888309; ...
%!      0.06160871602082477; ...
%!      0.0725880402231166; ...
%!      -0.2274936863299576; ...
%!      0.1022473973980562; ...
%!      -0.0526889967689672; ...
%!      0.04236597433282478; ...
%!      -0.03392636485319551; ...
%!      0.03013070722156086; ...
%!      -0.0296170703853042; ...
%!      0.02880801137929008; ...
%!      -0.04321809837188419; ...
%!      0.07296656815564413; ...
%!      -0.01716116565741117; ...
%!      -0.01024620657077106; ...
%!      0.001461842444097427; ...
%!      -0.001033490796384705; ...
%!      -2.754155888056701*10^-5; ...
%!      2.518658395698518*10^-5; ...
%!      7.901597126014116*10^-4; ...
%!      -9.335355913526482*10^-4; ...
%!      0.005456669361565423; ...
%!      0.007604249646274445; ...
%!      -0.02682003942293659; ...
%!      0.0131261469494572; ...
%!      -0.007194557768823587; ...
%!      0.006045470501789038; ...
%!      -0.004989519282046221; ...
%!      0.004515902634334341; ...
%!      -0.004481128825392679; ...
%!      0.004364459613084872; ...
%!      -0.006509135150844266; ...
%!      0.01085376379283025; ...
%!      -0.002505878517562597; ...
%!      -0.001460217755938636; ...
%!      2.021831084524683*10^-4; ...
%!      -1.379419236884723*10^-4; ...
%!      -3.527206528143375*10^-6; ...
%!      3.076563909696012*10^-6; ...
%!      9.147404918569031*10^-5; ...
%!      -1.017156221966768*10^-4; ...
%!      5.55274557745488*10^-4; ...
%!      7.163914730616564*10^-4; ...
%!      -0.002315442483754842; ...
%!      0.001025890224635668; ...
%!      -5.014856683353488*10^-4; ...
%!      3.687425121492475*10^-4; ...
%!      -2.596982003780794*10^-4];
%! N = 99;
%! f = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
%! A = [0.2 -0.4 0.6 -0.8 1];
%! h = fir (N, f, A, kaiser (100, 5))';
%! assert (x, h, 1e-15);

%% error tests
%!error h = fir ()
%!error h = fir (10)
%!error h = fir ([0 0.3])
%!error h = fir (9, [-0.1 1])
%!error h = fir (9, hamming (10))
%!error h = fir (9, [0.3 0.3])
%!error h = fir (9, [0 0.3 0.6], 1)
%!error h = fir (9, [0 0.3], hamming (11))
