signal-1.4.1/src/remez.cc

/*

Copyright (C) 1995, 1998 Jake Janovetz <janovetz@uiuc.edu>
Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
Copyright (C) 2000 Kai Habel <kahacjde@linux.zrz.tu-berlin.de>

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; see the file COPYING.  If not, see
<https://www.gnu.org/licenses/>.

*/

/**************************************************************************
 *  There appear to be some problems with the routine Search. See comments
 *  therein [search for PAK:].  I haven't looked closely at the rest
 *  of the code---it may also have some problems.
 *************************************************************************/

#include <octave/oct.h>
#include <cmath>

#ifndef OCTAVE_LOCAL_BUFFER
#include <vector>
#define OCTAVE_LOCAL_BUFFER(T, buf, size) \
  std::vector<T> buf ## _vector (size); \
  T *buf = &(buf ## _vector[0])
#endif

#include "octave-compat.h"

#define CONST const
#define BANDPASS       1
#define DIFFERENTIATOR 2
#define HILBERT        3

#define NEGATIVE       0
#define POSITIVE       1

#define Pi             3.1415926535897932
#define Pi2            6.2831853071795865

#define GRIDDENSITY    16
#define MAXITERATIONS  40

/*******************
 * CreateDenseGrid
 *=================
 * Creates the dense grid of frequencies from the specified bands.
 * Also creates the Desired Frequency Response function (D[]) and
 * the Weight function (W[]) on that dense grid
 *
 *
 * INPUT:
 * ------
 * int      r        - 1/2 the number of filter coefficients
 * int      numtaps  - Number of taps in the resulting filter
 * int      numband  - Number of bands in user specification
 * double   bands[]  - User-specified band edges [2*numband]
 * double   des[]    - Desired response per band [2*numband]
 * double   weight[] - Weight per band [numband]
 * int      symmetry - Symmetry of filter - used for grid check
 * int      griddensity
 *
 * OUTPUT:
 * -------
 * int    gridsize   - Number of elements in the dense frequency grid
 * double Grid[]     - Frequencies (0 to 0.5) on the dense grid [gridsize]
 * double D[]        - Desired response on the dense grid [gridsize]
 * double W[]        - Weight function on the dense grid [gridsize]
 *******************/

void CreateDenseGrid(int r, int numtaps, int numband, const double bands[],
                     const double des[], const double weight[], int gridsize,
                     double Grid[], double D[], double W[],
                     int symmetry, int griddensity)
{
  int i, j, k, band;
  double delf, lowf, highf, grid0;

  delf = 0.5/(griddensity*r);

  /*
   * For differentiator, hilbert,
   *   symmetry is odd and Grid[0] = max(delf, bands[0])
   */
  grid0 = (symmetry == NEGATIVE) && (delf > bands[0]) ? delf : bands[0];

  j=0;
  for (band=0; band < numband; band++)
    {
      lowf = (band==0 ? grid0 : bands[2*band]);
      highf = bands[2*band + 1];
      k = (int)((highf - lowf)/delf + 0.5);   /* .5 for rounding */
      if (band == 0 && symmetry == NEGATIVE)
        k--;
      for (i=0; i<k; i++)
        {
          D[j] = des[2*band] + i*(des[2*band+1]-des[2*band])/(k-1);
          W[j] = weight[band];
          Grid[j] = lowf;
          lowf += delf;
          j++;
        }
      Grid[j-1] = highf;
    }

/*
 * Similar to above, if odd symmetry, last grid point can't be .5
 *  - but, if there are even taps, leave the last grid point at .5
 */
  if ((symmetry == NEGATIVE) &&
       (Grid[gridsize-1] > (0.5 - delf)) &&
       (numtaps % 2))
    {
      Grid[gridsize-1] = 0.5-delf;
    }
}


/********************
 * InitialGuess
 *==============
 * Places Extremal Frequencies evenly throughout the dense grid.
 *
 *
 * INPUT:
 * ------
 * int r        - 1/2 the number of filter coefficients
 * int gridsize - Number of elements in the dense frequency grid
 *
 * OUTPUT:
 * -------
 * int Ext[]    - Extremal indexes to dense frequency grid [r+1]
 ********************/

void InitialGuess(int r, int Ext[], int gridsize)
{
  int i;

  for (i=0; i<=r; i++)
    Ext[i] = i * (gridsize-1) / r;
}


/***********************
 * CalcParms
 *===========
 *
 *
 * INPUT:
 * ------
 * int    r      - 1/2 the number of filter coefficients
 * int    Ext[]  - Extremal indexes to dense frequency grid [r+1]
 * double Grid[] - Frequencies (0 to 0.5) on the dense grid [gridsize]
 * double D[]    - Desired response on the dense grid [gridsize]
 * double W[]    - Weight function on the dense grid [gridsize]
 *
 * OUTPUT:
 * -------
 * double ad[]   - 'b' in Oppenheim & Schafer [r+1]
 * double x[]    - [r+1]
 * double y[]    - 'C' in Oppenheim & Schafer [r+1]
 ***********************/

void CalcParms(int r, int Ext[], double Grid[], double D[], double W[],
               double ad[], double x[], double y[])
{
  int i, j, k, ld;
  double sign, xi, delta, denom, numer;

  /*
   * Find x[]
   */
  for (i=0; i<=r; i++)
    x[i] = cos(Pi2 * Grid[Ext[i]]);

  /*
   * Calculate ad[]  - Oppenheim & Schafer eq 7.132
   */
  ld = (r-1)/15 + 1;         /* Skips around to avoid round errors */
  for (i=0; i<=r; i++)
    {
      denom = 1.0;
      xi = x[i];
      for (j=0; j<ld; j++)
        {
          for (k=j; k<=r; k+=ld)
            if (k != i)
              denom *= 2.0*(xi - x[k]);
        }
      if (fabs(denom)<0.00001)
        denom = 0.00001;
      ad[i] = 1.0/denom;
    }

  /*
   * Calculate delta  - Oppenheim & Schafer eq 7.131
   */
  numer = denom = 0;
  sign = 1;
  for (i=0; i<=r; i++)
    {
      numer += ad[i] * D[Ext[i]];
      denom += sign * ad[i]/W[Ext[i]];
      sign = -sign;
    }
  delta = numer/denom;
  sign = 1;

  /*
   * Calculate y[]  - Oppenheim & Schafer eq 7.133b
   */
  for (i=0; i<=r; i++)
    {
      y[i] = D[Ext[i]] - sign * delta/W[Ext[i]];
      sign = -sign;
    }
}


/*********************
 * ComputeA
 *==========
 * Using values calculated in CalcParms, ComputeA calculates the
 * actual filter response at a given frequency (freq).  Uses
 * eq 7.133a from Oppenheim & Schafer.
 *
 *
 * INPUT:
 * ------
 * double freq - Frequency (0 to 0.5) at which to calculate A
 * int    r    - 1/2 the number of filter coefficients
 * double ad[] - 'b' in Oppenheim & Schafer [r+1]
 * double x[]  - [r+1]
 * double y[]  - 'C' in Oppenheim & Schafer [r+1]
 *
 * OUTPUT:
 * -------
 * Returns double value of A[freq]
 *********************/

double ComputeA(double freq, int r, double ad[], double x[], double y[])
{
  int i;
  double xc, c, denom, numer;

  denom = numer = 0;
  xc = cos(Pi2 * freq);
  for (i=0; i<=r; i++)
    {
      c = xc - x[i];
      if (fabs(c) < 1.0e-7)
        {
          numer = y[i];
          denom = 1;
          break;
        }
      c = ad[i]/c;
      denom += c;
      numer += c*y[i];
    }
  return numer/denom;
}


/************************
 * CalcError
 *===========
 * Calculates the Error function from the desired frequency response
 * on the dense grid (D[]), the weight function on the dense grid (W[]),
 * and the present response calculation (A[])
 *
 *
 * INPUT:
 * ------
 * int    r      - 1/2 the number of filter coefficients
 * double ad[]   - [r+1]
 * double x[]    - [r+1]
 * double y[]    - [r+1]
 * int gridsize  - Number of elements in the dense frequency grid
 * double Grid[] - Frequencies on the dense grid [gridsize]
 * double D[]    - Desired response on the dense grid [gridsize]
 * double W[]    - Weight function on the desnse grid [gridsize]
 *
 * OUTPUT:
 * -------
 * double E[]    - Error function on dense grid [gridsize]
 ************************/

void CalcError(int r, double ad[], double x[], double y[],
               int gridsize, double Grid[],
               double D[], double W[], double E[])
{
  int i;
  double A;

  for (i=0; i<gridsize; i++)
    {
      A = ComputeA(Grid[i], r, ad, x, y);
      E[i] = W[i] * (D[i] - A);
    }
}

/************************
 * Search
 *========
 * Searches for the maxima/minima of the error curve.  If more than
 * r+1 extrema are found, it uses the following heuristic (thanks
 * Chris Hanson):
 * 1) Adjacent non-alternating extrema deleted first.
 * 2) If there are more than one excess extrema, delete the
 *    one with the smallest error.  This will create a non-alternation
 *    condition that is fixed by 1).
 * 3) If there is exactly one excess extremum, delete the smaller
 *    of the first/last extremum
 *
 *
 * INPUT:
 * ------
 * int    r        - 1/2 the number of filter coefficients
 * int    Ext[]    - Indexes to Grid[] of extremal frequencies [r+1]
 * int    gridsize - Number of elements in the dense frequency grid
 * double E[]      - Array of error values.  [gridsize]
 * OUTPUT:
 * -------
 * int    Ext[]    - New indexes to extremal frequencies [r+1]
 ************************/
int Search(int r, int Ext[],
           int gridsize, double E[])
{
  int i, j, k, l, extra;     /* Counters */
  int up, alt;
  int *foundExt;             /* Array of found extremals */

  /*
   * Allocate enough space for found extremals.
   */
  foundExt = (int *)malloc((2*r) * sizeof(int));
  k = 0;

  /*
   * Check for extremum at 0.
   */
  if (((E[0]>0.0) && (E[0]>E[1])) ||
      ((E[0]<0.0) && (E[0]<E[1])))
    foundExt[k++] = 0;

  /*
   * Check for extrema inside dense grid
   */
  for (i=1; i<gridsize-1; i++)
    {
      if (((E[i]>=E[i-1]) && (E[i]>E[i+1]) && (E[i]>0.0)) ||
          ((E[i]<=E[i-1]) && (E[i]<E[i+1]) && (E[i]<0.0))) {
        // PAK: we sometimes get too many extremal frequencies
        if (k >= 2*r) return -3;
        foundExt[k++] = i;
      }
    }

  /*
   * Check for extremum at 0.5
   */
  j = gridsize-1;
  if (((E[j]>0.0) && (E[j]>E[j-1])) ||
       ((E[j]<0.0) && (E[j]<E[j-1]))) {
    if (k >= 2*r) return -3;
    foundExt[k++] = j;
  }

  // PAK: we sometimes get not enough extremal frequencies
  if (k < r+1) return -2;


  /*
   * Remove extra extremals
   */
  extra = k - (r+1);
  assert(extra >= 0);

  while (extra > 0)
    {
      if (E[foundExt[0]] > 0.0)
        up = 1;                /* first one is a maxima */
      else
        up = 0;                /* first one is a minima */

      l=0;
      alt = 1;
      for (j=1; j<k; j++)
        {
          if (fabs(E[foundExt[j]]) < fabs(E[foundExt[l]]))
            l = j;               /* new smallest error. */
          if ((up) && (E[foundExt[j]] < 0.0))
            up = 0;             /* switch to a minima */
          else if ((!up) && (E[foundExt[j]] > 0.0))
            up = 1;             /* switch to a maxima */
          else
            {
              alt = 0;
              // PAK: break now and you will delete the smallest overall
              // extremal.  If you want to delete the smallest of the
              // pair of non-alternating extremals, then you must do:
              //
              // if (fabs(E[foundExt[j]]) < fabs(E[foundExt[j-1]])) l=j;
              // else l=j-1;
              break;              /* Ooops, found two non-alternating */
            }                     /* extrema.  Delete smallest of them */
        }  /* if the loop finishes, all extrema are alternating */

      /*
       * If there's only one extremal and all are alternating,
       * delete the smallest of the first/last extremals.
       */
      if ((alt) && (extra == 1))
        {
          if (fabs(E[foundExt[k-1]]) < fabs(E[foundExt[0]]))
            /* Delete last extremal */
            l = k-1;
            // PAK: changed from l = foundExt[k-1];
          else
            /* Delete first extremal */
            l = 0;
            // PAK: changed from l = foundExt[0];
        }

      for (j=l; j<k-1; j++)        /* Loop that does the deletion */
        {
          foundExt[j] = foundExt[j+1];
          assert(foundExt[j]<gridsize);
        }
      k--;
      extra--;
    }

  for (i=0; i<=r; i++)
    {
      assert(foundExt[i]<gridsize);
      Ext[i] = foundExt[i];       /* Copy found extremals to Ext[] */
    }

  free(foundExt);
  return 0;
}


/*********************
 * FreqSample
 *============
 * Simple frequency sampling algorithm to determine the impulse
 * response h[] from A's found in ComputeA
 *
 *
 * INPUT:
 * ------
 * int      N        - Number of filter coefficients
 * double   A[]      - Sample points of desired response [N/2]
 * int      symmetry - Symmetry of desired filter
 *
 * OUTPUT:
 * -------
 * double h[] - Impulse Response of final filter [N]
 *********************/
void FreqSample(int N, double A[], double h[], int symm)
{
  int n, k;
  double x, val, M;

  M = (N-1.0)/2.0;
  if (symm == POSITIVE)
    {
      if (N%2)
        {
          for (n=0; n<N; n++)
            {
              val = A[0];
              x = Pi2 * (n - M)/N;
              for (k=1; k<=M; k++)
                val += 2.0 * A[k] * cos(x*k);
              h[n] = val/N;
            }
        }
      else
        {
          for (n=0; n<N; n++)
            {
              val = A[0];
              x = Pi2 * (n - M)/N;
              for (k=1; k<=(N/2-1); k++)
                val += 2.0 * A[k] * cos(x*k);
              h[n] = val/N;
            }
        }
    }
  else
    {
      if (N%2)
        {
          for (n=0; n<N; n++)
            {
              val = 0;
              x = Pi2 * (n - M)/N;
              for (k=1; k<=M; k++)
                val += 2.0 * A[k] * sin(x*k);
              h[n] = val/N;
            }
        }
      else
        {
          for (n=0; n<N; n++)
            {
              val = A[N/2] * sin(Pi * (n - M));
              x = Pi2 * (n - M)/N;
              for (k=1; k<=(N/2-1); k++)
                val += 2.0 * A[k] * sin(x*k);
              h[n] = val/N;
            }
        }
    }
}

/*******************
 * isDone
 *========
 * Checks to see if the error function is small enough to consider
 * the result to have converged.
 *
 * INPUT:
 * ------
 * int    r     - 1/2 the number of filter coefficients
 * int    Ext[] - Indexes to extremal frequencies [r+1]
 * double E[]   - Error function on the dense grid [gridsize]
 *
 * OUTPUT:
 * -------
 * Returns 1 if the result converged
 * Returns 0 if the result has not converged
 ********************/

int isDone(int r, int Ext[], double E[])
{
  int i;
  double min, max, current;

  min = max = fabs(E[Ext[0]]);
  for (i=1; i<=r; i++)
    {
      current = fabs(E[Ext[i]]);
      if (current < min)
        min = current;
      if (current > max)
        max = current;
    }
  return (((max-min)/max) < 0.0001);
}

/********************
 * remez
 *=======
 * Calculates the optimal (in the Chebyshev/minimax sense)
 * FIR filter impulse response given a set of band edges,
 * the desired response on those bands, and the weight given to
 * the error in those bands.
 *
 * INPUT:
 * ------
 * int     numtaps     - Number of filter coefficients
 * int     numband     - Number of bands in filter specification
 * double  bands[]     - User-specified band edges [2 * numband]
 * double  des[]       - User-specified band responses [numband]
 * double  weight[]    - User-specified error weights [numband]
 * int     type        - Type of filter
 *
 * OUTPUT:
 * -------
 * double h[]      - Impulse response of final filter [numtaps]
 * returns         - true on success, false on failure to converge
 ********************/

int remez(double h[], int numtaps,
          int numband, const double bands[],
          const double des[], const double weight[],
          int type, int griddensity)
{
  double *Grid, *W, *D, *E;
  int    i, iter, gridsize, r, *Ext;
  double *taps, c;
  double *x, *y, *ad;
  int    symmetry;

  if (type == BANDPASS)
    symmetry = POSITIVE;
  else
    symmetry = NEGATIVE;

  r = numtaps/2;                  /* number of extrema */
  if ((numtaps%2) && (symmetry == POSITIVE))
    r++;

  /*
   * Predict dense grid size in advance for memory allocation
   *   .5 is so we round up, not truncate
   */
  gridsize = 0;
  for (i=0; i<numband; i++)
    {
      gridsize += (int)(2*r*griddensity*(bands[2*i+1] - bands[2*i]) + .5);
    }
  if (symmetry == NEGATIVE)
    {
      gridsize--;
    }

  /*
   * Dynamically allocate memory for arrays with proper sizes
   */
  Grid = (double *)malloc(gridsize * sizeof(double));
  D = (double *)malloc(gridsize * sizeof(double));
  W = (double *)malloc(gridsize * sizeof(double));
  E = (double *)malloc(gridsize * sizeof(double));
  Ext = (int *)malloc((r+1) * sizeof(int));
  taps = (double *)malloc((r+1) * sizeof(double));
  x = (double *)malloc((r+1) * sizeof(double));
  y = (double *)malloc((r+1) * sizeof(double));
  ad = (double *)malloc((r+1) * sizeof(double));

  /*
   * Create dense frequency grid
   */
  CreateDenseGrid(r, numtaps, numband, bands, des, weight,
                  gridsize, Grid, D, W, symmetry, griddensity);
  InitialGuess(r, Ext, gridsize);

  /*
   * For Differentiator: (fix grid)
   */
  if (type == DIFFERENTIATOR)
    {
      for (i=0; i<gridsize; i++)
        {
          /* D[i] = D[i]*Grid[i]; */
          if (D[i] > 0.0001)
            W[i] = W[i]/Grid[i];
        }
    }

  /*
   * For odd or Negative symmetry filters, alter the
   * D[] and W[] according to Parks McClellan
   */
  if (symmetry == POSITIVE)
    {
      if (numtaps % 2 == 0)
        {
          for (i=0; i<gridsize; i++)
            {
              c = cos(Pi * Grid[i]);
              D[i] /= c;
              W[i] *= c;
            }
        }
    }
  else
    {
      if (numtaps % 2)
        {
          for (i=0; i<gridsize; i++)
            {
              c = sin(Pi2 * Grid[i]);
              D[i] /= c;
              W[i] *= c;
            }
        }
      else
        {
          for (i=0; i<gridsize; i++)
            {
              c = sin(Pi * Grid[i]);
              D[i] /= c;
              W[i] *= c;
            }
        }
    }

  /*
   * Perform the Remez Exchange algorithm
   */
  for (iter=0; iter<MAXITERATIONS; iter++)
    {
      CalcParms(r, Ext, Grid, D, W, ad, x, y);
      CalcError(r, ad, x, y, gridsize, Grid, D, W, E);
      int err = Search(r, Ext, gridsize, E);
      if (err) return err;
      for(i=0; i <= r; i++) assert(Ext[i]<gridsize);
      if (isDone(r, Ext, E))
        break;
    }

  CalcParms(r, Ext, Grid, D, W, ad, x, y);

  /*
   * Find the 'taps' of the filter for use with Frequency
   * Sampling.  If odd or Negative symmetry, fix the taps
   * according to Parks McClellan
   */
  for (i=0; i<=numtaps/2; i++)
    {
      if (symmetry == POSITIVE)
        {
          if (numtaps%2)
            c = 1;
          else
            c = cos(Pi * (double)i/numtaps);
        }
      else
        {
          if (numtaps%2)
            c = sin(Pi2 * (double)i/numtaps);
          else
            c = sin(Pi * (double)i/numtaps);
        }
      taps[i] = ComputeA((double)i/numtaps, r, ad, x, y)*c;
    }

  /*
   * Frequency sampling design with calculated taps
   */
  FreqSample(numtaps, taps, h, symmetry);

  /*
   * Delete allocated memory
   */
  free(Grid);
  free(W);
  free(D);
  free(E);
  free(Ext);
  free(x);
  free(y);
  free(ad);
  return iter<MAXITERATIONS?0:-1;
}


/* == Octave interface starts here ====================================== */

DEFUN_DLD (remez, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Loadable Function} {@var{b} =} remez (@var{n}, @var{f}, @var{a})\n\
@deftypefnx {Loadable Function} {@var{b} =} remez (@var{n}, @var{f}, @var{a}, @var{w})\n\
@deftypefnx {Loadable Function} {@var{b} =} remez (@var{n}, @var{f}, @var{a}, @var{w}, @var{ftype})\n\
@deftypefnx {Loadable Function} {@var{b} =} remez (@var{n}, @var{f}, @var{a}, @var{w}, @var{ftype}, @var{griddensity})\n\
Parks-McClellan optimal FIR filter design.\n\
@table @var\n\
@item n\n\
gives the number of taps in the returned filter\n\
@item f\n\
gives frequency at the band edges [b1 e1 b2 e2 b3 e3 @dots{}]\n\
@item a\n\
gives amplitude at the band edges [a(b1) a(e1) a(b2) a(e2) @dots{}]\n\
@item w\n\
gives weighting applied to each band\n\
@item ftype\n\
is \"bandpass\", \"hilbert\" or \"differentiator\"\n\
@item griddensity\n\
determines how accurately the filter will be\n\
constructed. The minimum value is 16, but higher numbers are\n\
slower to compute.\n\
@end table\n\
\n\
Frequency is in the range (0, 1), with 1 being the Nyquist frequency.\n\
@end deftypefn")
{
  octave_value_list retval;
  int i;

  int nargin = args.length();
  if (nargin < 3 || nargin > 6) {
    print_usage ();
    return retval;
  }

  int numtaps = octave::math::nint (args(0).double_value ()) + 1;
  if (numtaps < 4) {
    error("remez: number of taps must be an integer greater than 3");
    return retval;
  }

  ColumnVector o_bands(args(1).vector_value());
  int numbands = o_bands.numel()/2;
  OCTAVE_LOCAL_BUFFER(double, bands, numbands*2);
  if (numbands < 1 || o_bands.numel()%2 == 1) {
    error("remez: must have an even number of band edges");
    return retval;
  }
  for (i=1; i < o_bands.numel(); i++) {
    if (o_bands(i)<o_bands(i-1)) {
      error("band edges must be nondecreasing");
      return retval;
    }
  }
  if (o_bands(0) < 0 || o_bands(1) > 1) {
    error("band edges must be in the range [0,1]");
    return retval;
  }
  for(i=0; i < 2*numbands; i++) bands[i] = o_bands(i)/2.0;

  ColumnVector o_response(args(2).vector_value());
  OCTAVE_LOCAL_BUFFER (double, response, numbands*2);
  if (o_response.numel() != o_bands.numel()) {
    error("remez: must have one response magnitude for each band edge");
    return retval;
  }
  for(i=0; i < 2*numbands; i++) response[i] = o_response(i);

  std::string stype = std::string("bandpass");
  int density = 16;
  OCTAVE_LOCAL_BUFFER (double, weight, numbands);
  for (i=0; i < numbands; i++) weight[i] = 1.0;
  if (nargin > 3) {
    if (args(3).is_string())
      stype = args(3).string_value();
    else if (args(3).is_real_matrix() || args(3).is_real_scalar()) {
      ColumnVector o_weight(args(3).vector_value());
      if (o_weight.numel() != numbands) {
        error("remez: need one weight for each band [=length(band)/2]");
        return retval;
      }
      for (i=0; i < numbands; i++) weight[i] = o_weight(i);
    }
    else {
      error("remez: incorrect argument list");
      return retval;
    }
  }
  if (nargin > 4) {
    if (args(4).is_string() && !args(3).is_string())
      stype = args(4).string_value();
    else if (args(4).is_real_scalar())
      density = octave::math::nint (args(4).double_value ());
    else {
      error("remez: incorrect argument list");
      return retval;
    }
  }
  if (nargin > 5) {
    if (args(5).is_real_scalar()
        && !args(4).is_real_scalar())
      density = octave::math::nint (args(5).double_value ());
    else {
      error("remez: incorrect argument list");
      return retval;
    }
  }

  int itype;
  if (stype == "bandpass")
    itype = BANDPASS;
  else if (stype == "differentiator")
    itype = DIFFERENTIATOR;
  else if (stype == "hilbert")
    itype = HILBERT;
  else {
    error("remez: unknown ftype '%s'", stype.data());
    return retval;
  }

  if (density < 16) {
    error("remez: griddensity is too low; must be greater than 16");
    return retval;
  }

  OCTAVE_LOCAL_BUFFER (double, coeff, numtaps+5);
  int err = remez(coeff,numtaps,numbands,bands,response,weight,itype,density);

  if (err == -1)
    warning("remez: -- failed to converge -- returned filter may be bad.");
  else if (err == -2) {
    error("remez: insufficient extremals--cannot continue");
    return retval;
  }
  else if (err == -3) {
    error("remez: too many extremals--cannot continue");
    return retval;
  }

  ColumnVector h(numtaps);
  while(numtaps--) h(numtaps) = coeff[numtaps];

  return octave_value(h);
}

/*
%!test
%! b = [
%!    0.0415131831103279
%!    0.0581639884202646
%!   -0.0281579212691008
%!   -0.0535575358002337
%!   -0.0617245915143180
%!    0.0507753178978075
%!    0.2079018331396460
%!    0.3327160895375440
%!    0.3327160895375440
%!    0.2079018331396460
%!    0.0507753178978075
%!   -0.0617245915143180
%!   -0.0535575358002337
%!   -0.0281579212691008
%!    0.0581639884202646
%!    0.0415131831103279];
%! assert(remez(15,[0,0.3,0.4,1],[1,1,0,0]),b,1e-14);

 */