gr_remez.cc

这是用python语言写的一个数字广播的信号处理工具包。利用它

         {
            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 coeffiecients
 * 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
 ********************/

static bool
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 reponse 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 [2 * 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
 ********************/

static 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(int 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;
}

//////////////////////////////////////////////////////////////////////////////
//
//			    GNU Radio interface
//
//////////////////////////////////////////////////////////////////////////////


static void
punt (const std::string msg) 
{
  std::cerr << msg << '\n';
  throw std::runtime_error (msg);
}

std::vector<double>
gr_remez (int order,
	  const std::vector<double> &arg_bands,
	  const std::vector<double> &arg_response,
	  const std::vector<double> &arg_weight,
	  const std::string filter_type,
	  int grid_density
	  ) throw (std::runtime_error)
{
  int numtaps = order + 1;
  if (numtaps < 4)
    punt ("gr_remez: number of taps must be >= 3");

  int numbands = arg_bands.size () / 2;
  LOCAL_BUFFER (double, bands, numbands * 2);
  if (numbands < 1 || arg_bands.size () % 2 == 1)
    punt ("gr_remez: must have an even number of band edges");

  for (unsigned int i = 1; i < arg_bands.size (); i++){
    if (arg_bands[i] < arg_bands[i-1])
      punt ("gr_remez: band edges must be nondecreasing");
  }

  if (arg_bands[0] < 0 || arg_bands[arg_bands.size () - 1] > 1)
    punt ("gr_remez: band edges must be in the range [0,1]");

  for (int i = 0; i < 2 * numbands; i++)
    bands[i] = arg_bands[i] / 2;		// FIXME why / 2?

  LOCAL_BUFFER (double, response, numbands * 2);
  if (arg_response.size () != arg_bands.size ())
    punt ("gr_remez: must have one response magnitude for each band edge");

  for (int i = 0; i < 2 * numbands; i++)
    response[i] = arg_response[i];

  LOCAL_BUFFER (double, weight, numbands);
  for (int i = 0; i < numbands; i++)
    weight[i] = 1.0;

  if (arg_weight.size () != 0){
    if ((int) arg_weight.size () != numbands)
      punt ("gr_remez: need one weight for each band [=length(band)/2]");
    for (int i = 0; i < numbands; i++)
      weight[i] = arg_weight [i];
  }
  
  int itype = 0;
  if (filter_type == "bandpass") 
    itype = BANDPASS;
  else if (filter_type == "differentiator") 
    itype = DIFFERENTIATOR;
  else if (filter_type == "hilbert") 
    itype = HILBERT;
  else
    punt ("gr_remez: unknown ftype '" + filter_type + "'");

  if (grid_density < 16)
    punt ("gr_remez: grid_density is too low; must be >= 16");

  LOCAL_BUFFER (double, coeff, numtaps + 5);		// FIXME why + 5?
  int err = remez (coeff, numtaps, numbands,
		   bands, response, weight, itype, grid_density);

  if (err == -1)
    punt ("gr_remez: failed to converge");

  if (err == -2)
    punt ("gr_remez: insufficient extremals -- cannot continue");

  if (err == -3)
    punt ("gr_remez: too many extremals -- cannot continue");

  return std::vector<double> (&coeff[0], &coeff[numtaps]);
}



#if 0
/* == Octave interface starts here ====================================== */

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

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

  int numtaps = NINT (args(0).double_value()) + 1; // #coeff = filter order+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.length()/2;
  OCTAVE_LOCAL_BUFFER(double, bands, numbands*2);
  if (numbands < 1 || o_bands.length()%2 == 1) {
    error("remez: must have an even number of band edges");
    return retval;
  }
  for (i=1; i < o_bands.length(); 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.length() != o_bands.length()) {
    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_real_matrix()) {
      ColumnVector o_weight(args(3).vector_value());
      if (o_weight.length() != 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 if (args(3).is_string())
      stype = args(3).string_value();
    else if (args(3).is_real_scalar())
      density = NINT(args(3).double_value());
    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() && !args(3).is_real_scalar())
      density = 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() 
	&& !args(3).is_real_scalar())
      density = NINT(args(4).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);

 */

#endif