git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/general/gr_remez.cc

   1 /**************************************************************************
   2  * Parks-McClellan algorithm for FIR filter design (C version)
   3  *-------------------------------------------------
   4  *  Copyright (c) 1995,1998  Jake Janovetz (janovetz@uiuc.edu)
   5  *  Copyright (c) 2004  Free Software Foundation, Inc.
   6  *
   7  *  This library is free software; you can redistribute it and/or
   8  *  modify it under the terms of the GNU Library General Public
   9  *  License as published by the Free Software Foundation; either
  10  *  version 2 of the License, or (at your option) any later version.
  11  *
  12  *  This library is distributed in the hope that it will be useful,
  13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  *  Library General Public License for more details.
  16  *
  17  *  You should have received a copy of the GNU Library General Public
  18  *  License along with this library; if not, write to the Free
  19  *  Foundation, Inc., 51 Franklin Street, Boston, MA  02110-1301  USA
  20  *
  21  *
  22  *  Sep 1999 - Paul Kienzle (pkienzle@cs.indiana.edu)
  23  *      Modified for use in octave as a replacement for the matlab function
  24  *      remez.mex.  In particular, magnitude responses are required for all
  25  *      band edges rather than one per band, griddensity is a parameter,
  26  *      and errors are returned rather than printed directly.
  27  *  Mar 2000 - Kai Habel (kahacjde@linux.zrz.tu-berlin.de)
  28  *      Change: ColumnVector x=arg(i).vector_value();
  29  *      to: ColumnVector x(arg(i).vector_value());
  30  *  There appear to be some problems with the routine Search. See comments
  31  *  therein [search for PAK:].  I haven't looked closely at the rest
  32  *  of the code---it may also have some problems.
  33  *************************************************************************/
  34
  35 /*
  36  * This code was extracted from octave.sf.net, and wrapped with
  37  * GNU Radio glue.
  38  */
  39
  40 #ifdef HAVE_CONFIG_H
  41 #include "config.h"
  42 #endif
  43 #include <gr_remez.h>
  44 #include <cmath>
  45 #include <assert.h>
  46 #include <iostream>
  47
  48
  49 #ifndef LOCAL_BUFFER
  50 #include <vector>
  51 #define LOCAL_BUFFER(T, buf, size) \
  52   std::vector<T> buf ## _vector (size); \
  53   T *buf = &(buf ## _vector[0])
  54 #endif
  55
  56
  57 #define CONST const
  58 #define BANDPASS       1
  59 #define DIFFERENTIATOR 2
  60 #define HILBERT        3
  61
  62 #define NEGATIVE       0
  63 #define POSITIVE       1
  64
  65 #define Pi             3.14159265358979323846
  66 #define Pi2            (2*Pi)
  67
  68 #define GRIDDENSITY    16
  69 #define MAXITERATIONS  40
  70
  71 /*******************
  72  * CreateDenseGrid
  73  *=================
  74  * Creates the dense grid of frequencies from the specified bands.
  75  * Also creates the Desired Frequency Response function (D[]) and
  76  * the Weight function (W[]) on that dense grid
  77  *
  78  *
  79  * INPUT:
  80  * ------
  81  * int      r        - 1/2 the number of filter coefficients
  82  * int      numtaps  - Number of taps in the resulting filter
  83  * int      numband  - Number of bands in user specification
  84  * double   bands[]  - User-specified band edges [2*numband]
  85  * double   des[]    - Desired response per band [2*numband]
  86  * double   weight[] - Weight per band [numband]
  87  * int      symmetry - Symmetry of filter - used for grid check
  88  * int      griddensity
  89  *
  90  * OUTPUT:
  91  * -------
  92  * int    gridsize   - Number of elements in the dense frequency grid
  93  * double Grid[]     - Frequencies (0 to 0.5) on the dense grid [gridsize]
  94  * double D[]        - Desired response on the dense grid [gridsize]
  95  * double W[]        - Weight function on the dense grid [gridsize]
  96  *******************/
  97
  98 static void
  99 CreateDenseGrid (int r, int numtaps, int numband, const double bands[],
 100                  const double des[], const double weight[], int gridsize,
 101                  double Grid[], double D[], double W[],
 102                  int symmetry, int griddensity)
 103 {
 104    int i, j, k, band;
 105    double delf, lowf, highf, grid0;
 106
 107    delf = 0.5/(griddensity*r);
 108
 109 /*
 110  * For differentiator, hilbert,
 111  *   symmetry is odd and Grid[0] = max(delf, bands[0])
 112  */
 113    grid0 = (symmetry == NEGATIVE) && (delf > bands[0]) ? delf : bands[0];
 114
 115    j=0;
 116    for (band=0; band < numband; band++)
 117    {
 118       lowf = (band==0 ? grid0 : bands[2*band]);
 119       highf = bands[2*band + 1];
 120       k = (int)((highf - lowf)/delf + 0.5);   /* .5 for rounding */
 121       for (i=0; i<k; i++)
 122       {
 123          D[j] = des[2*band] + i*(des[2*band+1]-des[2*band])/(k-1);
 124          W[j] = weight[band];
 125          Grid[j] = lowf;
 126          lowf += delf;
 127          j++;
 128       }
 129       Grid[j-1] = highf;
 130    }
 131
 132 /*
 133  * Similar to above, if odd symmetry, last grid point can't be .5
 134  *  - but, if there are even taps, leave the last grid point at .5
 135  */
 136    if ((symmetry == NEGATIVE) &&
 137        (Grid[gridsize-1] > (0.5 - delf)) &&
 138        (numtaps % 2))
 139    {
 140       Grid[gridsize-1] = 0.5-delf;
 141    }
 142 }
 143
 144
 145 /********************
 146  * InitialGuess
 147  *==============
 148  * Places Extremal Frequencies evenly throughout the dense grid.
 149  *
 150  *
 151  * INPUT:
 152  * ------
 153  * int r        - 1/2 the number of filter coefficients
 154  * int gridsize - Number of elements in the dense frequency grid
 155  *
 156  * OUTPUT:
 157  * -------
 158  * int Ext[]    - Extremal indexes to dense frequency grid [r+1]
 159  ********************/
 160
 161 static void
 162 InitialGuess (int r, int Ext[], int gridsize)
 163 {
 164    int i;
 165
 166    for (i=0; i<=r; i++)
 167       Ext[i] = i * (gridsize-1) / r;
 168 }
 169
 170
 171 /***********************
 172  * CalcParms
 173  *===========
 174  *
 175  *
 176  * INPUT:
 177  * ------
 178  * int    r      - 1/2 the number of filter coefficients
 179  * int    Ext[]  - Extremal indexes to dense frequency grid [r+1]
 180  * double Grid[] - Frequencies (0 to 0.5) on the dense grid [gridsize]
 181  * double D[]    - Desired response on the dense grid [gridsize]
 182  * double W[]    - Weight function on the dense grid [gridsize]
 183  *
 184  * OUTPUT:
 185  * -------
 186  * double ad[]   - 'b' in Oppenheim & Schafer [r+1]
 187  * double x[]    - [r+1]
 188  * double y[]    - 'C' in Oppenheim & Schafer [r+1]
 189  ***********************/
 190
 191 static void
 192 CalcParms (int r, int Ext[], double Grid[], double D[], double W[],
 193            double ad[], double x[], double y[])
 194 {
 195    int i, j, k, ld;
 196    double sign, xi, delta, denom, numer;
 197
 198 /*
 199  * Find x[]
 200  */
 201    for (i=0; i<=r; i++)
 202       x[i] = cos(Pi2 * Grid[Ext[i]]);
 203
 204 /*
 205  * Calculate ad[]  - Oppenheim & Schafer eq 7.132
 206  */
 207    ld = (r-1)/15 + 1;         /* Skips around to avoid round errors */
 208    for (i=0; i<=r; i++)
 209    {
 210        denom = 1.0;
 211        xi = x[i];
 212        for (j=0; j<ld; j++)
 213        {
 214           for (k=j; k<=r; k+=ld)
 215              if (k != i)
 216                 denom *= 2.0*(xi - x[k]);
 217        }
 218        if (fabs(denom)<0.00001)
 219           denom = 0.00001;
 220        ad[i] = 1.0/denom;
 221    }
 222
 223 /*
 224  * Calculate delta  - Oppenheim & Schafer eq 7.131
 225  */
 226    numer = denom = 0;
 227    sign = 1;
 228    for (i=0; i<=r; i++)
 229    {
 230       numer += ad[i] * D[Ext[i]];
 231       denom += sign * ad[i]/W[Ext[i]];
 232       sign = -sign;
 233    }
 234    delta = numer/denom;
 235    sign = 1;
 236
 237 /*
 238  * Calculate y[]  - Oppenheim & Schafer eq 7.133b
 239  */
 240    for (i=0; i<=r; i++)
 241    {
 242       y[i] = D[Ext[i]] - sign * delta/W[Ext[i]];
 243       sign = -sign;
 244    }
 245 }
 246
 247
 248 /*********************
 249  * ComputeA
 250  *==========
 251  * Using values calculated in CalcParms, ComputeA calculates the
 252  * actual filter response at a given frequency (freq).  Uses
 253  * eq 7.133a from Oppenheim & Schafer.
 254  *
 255  *
 256  * INPUT:
 257  * ------
 258  * double freq - Frequency (0 to 0.5) at which to calculate A
 259  * int    r    - 1/2 the number of filter coefficients
 260  * double ad[] - 'b' in Oppenheim & Schafer [r+1]
 261  * double x[]  - [r+1]
 262  * double y[]  - 'C' in Oppenheim & Schafer [r+1]
 263  *
 264  * OUTPUT:
 265  * -------
 266  * Returns double value of A[freq]
 267  *********************/
 268
 269 static double
 270 ComputeA (double freq, int r, double ad[], double x[], double y[])
 271 {
 272    int i;
 273    double xc, c, denom, numer;
 274
 275    denom = numer = 0;
 276    xc = cos(Pi2 * freq);
 277    for (i=0; i<=r; i++)
 278    {
 279       c = xc - x[i];
 280       if (fabs(c) < 1.0e-7)
 281       {
 282          numer = y[i];
 283          denom = 1;
 284          break;
 285       }
 286       c = ad[i]/c;
 287       denom += c;
 288       numer += c*y[i];
 289    }
 290    return numer/denom;
 291 }
 292
 293
 294 /************************
 295  * CalcError
 296  *===========
 297  * Calculates the Error function from the desired frequency response
 298  * on the dense grid (D[]), the weight function on the dense grid (W[]),
 299  * and the present response calculation (A[])
 300  *
 301  *
 302  * INPUT:
 303  * ------
 304  * int    r      - 1/2 the number of filter coefficients
 305  * double ad[]   - [r+1]
 306  * double x[]    - [r+1]
 307  * double y[]    - [r+1]
 308  * int gridsize  - Number of elements in the dense frequency grid
 309  * double Grid[] - Frequencies on the dense grid [gridsize]
 310  * double D[]    - Desired response on the dense grid [gridsize]
 311  * double W[]    - Weight function on the desnse grid [gridsize]
 312  *
 313  * OUTPUT:
 314  * -------
 315  * double E[]    - Error function on dense grid [gridsize]
 316  ************************/
 317
 318 static void
 319 CalcError (int r, double ad[], double x[], double y[],
 320            int gridsize, double Grid[],
 321            double D[], double W[], double E[])
 322 {
 323    int i;
 324    double A;
 325
 326    for (i=0; i<gridsize; i++)
 327    {
 328       A = ComputeA(Grid[i], r, ad, x, y);
 329       E[i] = W[i] * (D[i] - A);
 330    }
 331 }
 332
 333 /************************
 334  * Search
 335  *========
 336  * Searches for the maxima/minima of the error curve.  If more than
 337  * r+1 extrema are found, it uses the following heuristic (thanks
 338  * Chris Hanson):
 339  * 1) Adjacent non-alternating extrema deleted first.
 340  * 2) If there are more than one excess extrema, delete the
 341  *    one with the smallest error.  This will create a non-alternation
 342  *    condition that is fixed by 1).
 343  * 3) If there is exactly one excess extremum, delete the smaller
 344  *    of the first/last extremum
 345  *
 346  *
 347  * INPUT:
 348  * ------
 349  * int    r        - 1/2 the number of filter coefficients
 350  * int    Ext[]    - Indexes to Grid[] of extremal frequencies [r+1]
 351  * int    gridsize - Number of elements in the dense frequency grid
 352  * double E[]      - Array of error values.  [gridsize]
 353  * OUTPUT:
 354  * -------
 355  * int    Ext[]    - New indexes to extremal frequencies [r+1]
 356  ************************/
 357 static int
 358 Search (int r, int Ext[],
 359         int gridsize, double E[])
 360 {
 361    int i, j, k, l, extra;     /* Counters */
 362    int up, alt;
 363    int *foundExt;             /* Array of found extremals */
 364
 365 /*
 366  * Allocate enough space for found extremals.
 367  */
 368    foundExt = (int *)malloc((2*r) * sizeof(int));
 369    k = 0;
 370
 371 /*
 372  * Check for extremum at 0.
 373  */
 374    if (((E[0]>0.0) && (E[0]>E[1])) ||
 375        ((E[0]<0.0) && (E[0]<E[1])))
 376       foundExt[k++] = 0;
 377
 378 /*
 379  * Check for extrema inside dense grid
 380  */
 381    for (i=1; i<gridsize-1; i++)
 382    {
 383       if (((E[i]>=E[i-1]) && (E[i]>E[i+1]) && (E[i]>0.0)) ||
 384           ((E[i]<=E[i-1]) && (E[i]<E[i+1]) && (E[i]<0.0))) {
 385         // PAK: we sometimes get too many extremal frequencies
 386         if (k >= 2*r) return -3;
 387         foundExt[k++] = i;
 388       }
 389    }
 390
 391 /*
 392  * Check for extremum at 0.5
 393  */
 394    j = gridsize-1;
 395    if (((E[j]>0.0) && (E[j]>E[j-1])) ||
 396        ((E[j]<0.0) && (E[j]<E[j-1]))) {
 397      if (k >= 2*r) return -3;
 398      foundExt[k++] = j;
 399    }
 400
 401    // PAK: we sometimes get not enough extremal frequencies
 402    if (k < r+1) return -2;
 403
 404
 405 /*
 406  * Remove extra extremals
 407  */
 408    extra = k - (r+1);
 409    assert(extra >= 0);
 410
 411    while (extra > 0)
 412    {
 413       if (E[foundExt[0]] > 0.0)
 414          up = 1;                /* first one is a maxima */
 415       else
 416          up = 0;                /* first one is a minima */
 417
 418       l=0;
 419       alt = 1;
 420       for (j=1; j<k; j++)
 421       {
 422          if (fabs(E[foundExt[j]]) < fabs(E[foundExt[l]]))
 423             l = j;               /* new smallest error. */
 424          if ((up) && (E[foundExt[j]] < 0.0))
 425             up = 0;             /* switch to a minima */
 426          else if ((!up) && (E[foundExt[j]] > 0.0))
 427             up = 1;             /* switch to a maxima */
 428          else
 429          {
 430             alt = 0;
 431             // PAK: break now and you will delete the smallest overall
 432             // extremal.  If you want to delete the smallest of the
 433             // pair of non-alternating extremals, then you must do:
 434             //
 435             // if (fabs(E[foundExt[j]]) < fabs(E[foundExt[j-1]])) l=j;
 436             // else l=j-1;
 437             break;              /* Ooops, found two non-alternating */
 438          }                      /* extrema.  Delete smallest of them */
 439       }  /* if the loop finishes, all extrema are alternating */
 440
 441 /*
 442  * If there's only one extremal and all are alternating,
 443  * delete the smallest of the first/last extremals.
 444  */
 445       if ((alt) && (extra == 1))
 446       {
 447          if (fabs(E[foundExt[k-1]]) < fabs(E[foundExt[0]]))
 448            /* Delete last extremal */
 449            l = k-1;
 450            // PAK: changed from l = foundExt[k-1];
 451          else
 452            /* Delete first extremal */
 453            l = 0;
 454            // PAK: changed from l = foundExt[0];
 455       }
 456
 457       for (j=l; j<k-1; j++)        /* Loop that does the deletion */
 458       {
 459          foundExt[j] = foundExt[j+1];
 460          assert(foundExt[j]<gridsize);
 461       }
 462       k--;
 463       extra--;
 464    }
 465
 466    for (i=0; i<=r; i++)
 467    {
 468       assert(foundExt[i]<gridsize);
 469       Ext[i] = foundExt[i];       /* Copy found extremals to Ext[] */
 470    }
 471
 472    free(foundExt);
 473    return 0;
 474 }
 475
 476
 477 /*********************
 478  * FreqSample
 479  *============
 480  * Simple frequency sampling algorithm to determine the impulse
 481  * response h[] from A's found in ComputeA
 482  *
 483  *
 484  * INPUT:
 485  * ------
 486  * int      N        - Number of filter coefficients
 487  * double   A[]      - Sample points of desired response [N/2]
 488  * int      symmetry - Symmetry of desired filter
 489  *
 490  * OUTPUT:
 491  * -------
 492  * double h[] - Impulse Response of final filter [N]
 493  *********************/
 494 static void
 495 FreqSample (int N, double A[], double h[], int symm)
 496 {
 497    int n, k;
 498    double x, val, M;
 499
 500    M = (N-1.0)/2.0;
 501    if (symm == POSITIVE)
 502    {
 503       if (N%2)
 504       {
 505          for (n=0; n<N; n++)
 506          {
 507             val = A[0];
 508             x = Pi2 * (n - M)/N;
 509             for (k=1; k<=M; k++)
 510                val += 2.0 * A[k] * cos(x*k);
 511             h[n] = val/N;
 512          }
 513       }
 514       else
 515       {
 516          for (n=0; n<N; n++)
 517          {
 518             val = A[0];
 519             x = Pi2 * (n - M)/N;
 520             for (k=1; k<=(N/2-1); k++)
 521                val += 2.0 * A[k] * cos(x*k);
 522             h[n] = val/N;
 523          }
 524       }
 525    }
 526    else
 527    {
 528       if (N%2)
 529       {
 530          for (n=0; n<N; n++)
 531          {
 532             val = 0;
 533             x = Pi2 * (n - M)/N;
 534             for (k=1; k<=M; k++)
 535                val += 2.0 * A[k] * sin(x*k);
 536             h[n] = val/N;
 537          }
 538       }
 539       else
 540       {
 541           for (n=0; n<N; n++)
 542           {
 543              val = A[N/2] * sin(Pi * (n - M));
 544              x = Pi2 * (n - M)/N;
 545              for (k=1; k<=(N/2-1); k++)
 546                 val += 2.0 * A[k] * sin(x*k);
 547              h[n] = val/N;
 548           }
 549       }
 550    }
 551 }
 552
 553 /*******************
 554  * isDone
 555  *========
 556  * Checks to see if the error function is small enough to consider
 557  * the result to have converged.
 558  *
 559  * INPUT:
 560  * ------
 561  * int    r     - 1/2 the number of filter coeffiecients
 562  * int    Ext[] - Indexes to extremal frequencies [r+1]
 563  * double E[]   - Error function on the dense grid [gridsize]
 564  *
 565  * OUTPUT:
 566  * -------
 567  * Returns 1 if the result converged
 568  * Returns 0 if the result has not converged
 569  ********************/
 570
 571 static bool
 572 isDone (int r, int Ext[], double E[])
 573 {
 574    int i;
 575    double min, max, current;
 576
 577    min = max = fabs(E[Ext[0]]);
 578    for (i=1; i<=r; i++)
 579    {
 580       current = fabs(E[Ext[i]]);
 581       if (current < min)
 582          min = current;
 583       if (current > max)
 584          max = current;
 585    }
 586    return (((max-min)/max) < 0.0001);
 587 }
 588
 589 /********************
 590  * remez
 591  *=======
 592  * Calculates the optimal (in the Chebyshev/minimax sense)
 593  * FIR filter impulse response given a set of band edges,
 594  * the desired reponse on those bands, and the weight given to
 595  * the error in those bands.
 596  *
 597  * INPUT:
 598  * ------
 599  * int     numtaps     - Number of filter coefficients
 600  * int     numband     - Number of bands in filter specification
 601  * double  bands[]     - User-specified band edges [2 * numband]
 602  * double  des[]       - User-specified band responses [2 * numband]
 603  * double  weight[]    - User-specified error weights [numband]
 604  * int     type        - Type of filter
 605  *
 606  * OUTPUT:
 607  * -------
 608  * double h[]      - Impulse response of final filter [numtaps]
 609  * returns         - true on success, false on failure to converge
 610  ********************/
 611
 612 static int
 613 remez (double h[], int numtaps,
 614        int numband, const double bands[],
 615        const double des[], const double weight[],
 616        int type, int griddensity)
 617 {
 618    double *Grid, *W, *D, *E;
 619    int    i, iter, gridsize, r, *Ext;
 620    double *taps, c;
 621    double *x, *y, *ad;
 622    int    symmetry;
 623
 624    if (type == BANDPASS)
 625       symmetry = POSITIVE;
 626    else
 627       symmetry = NEGATIVE;
 628
 629    r = numtaps/2;                  /* number of extrema */
 630    if ((numtaps%2) && (symmetry == POSITIVE))
 631       r++;
 632
 633 /*
 634  * Predict dense grid size in advance for memory allocation
 635  *   .5 is so we round up, not truncate
 636  */
 637    gridsize = 0;
 638    for (i=0; i<numband; i++)
 639    {
 640       gridsize += (int)(2*r*griddensity*(bands[2*i+1] - bands[2*i]) + .5);
 641    }
 642    if (symmetry == NEGATIVE)
 643    {
 644       gridsize--;
 645    }
 646
 647 /*
 648  * Dynamically allocate memory for arrays with proper sizes
 649  */
 650    Grid = (double *)malloc(gridsize * sizeof(double));
 651    D = (double *)malloc(gridsize * sizeof(double));
 652    W = (double *)malloc(gridsize * sizeof(double));
 653    E = (double *)malloc(gridsize * sizeof(double));
 654    Ext = (int *)malloc((r+1) * sizeof(int));
 655    taps = (double *)malloc((r+1) * sizeof(double));
 656    x = (double *)malloc((r+1) * sizeof(double));
 657    y = (double *)malloc((r+1) * sizeof(double));
 658    ad = (double *)malloc((r+1) * sizeof(double));
 659
 660 /*
 661  * Create dense frequency grid
 662  */
 663    CreateDenseGrid(r, numtaps, numband, bands, des, weight,
 664                    gridsize, Grid, D, W, symmetry, griddensity);
 665    InitialGuess(r, Ext, gridsize);
 666
 667 /*
 668  * For Differentiator: (fix grid)
 669  */
 670    if (type == DIFFERENTIATOR)
 671    {
 672       for (i=0; i<gridsize; i++)
 673       {
 674 /* D[i] = D[i]*Grid[i]; */
 675          if (D[i] > 0.0001)
 676             W[i] = W[i]/Grid[i];
 677       }
 678    }
 679
 680 /*
 681  * For odd or Negative symmetry filters, alter the
 682  * D[] and W[] according to Parks McClellan
 683  */
 684    if (symmetry == POSITIVE)
 685    {
 686       if (numtaps % 2 == 0)
 687       {
 688          for (i=0; i<gridsize; i++)
 689          {
 690             c = cos(Pi * Grid[i]);
 691             D[i] /= c;
 692             W[i] *= c;
 693          }
 694       }
 695    }
 696    else
 697    {
 698       if (numtaps % 2)
 699       {
 700          for (i=0; i<gridsize; i++)
 701          {
 702             c = sin(Pi2 * Grid[i]);
 703             D[i] /= c;
 704             W[i] *= c;
 705          }
 706       }
 707       else
 708       {
 709          for (i=0; i<gridsize; i++)
 710          {
 711             c = sin(Pi * Grid[i]);
 712             D[i] /= c;
 713             W[i] *= c;
 714          }
 715       }
 716    }
 717
 718 /*
 719  * Perform the Remez Exchange algorithm
 720  */
 721    for (iter=0; iter<MAXITERATIONS; iter++)
 722    {
 723       CalcParms(r, Ext, Grid, D, W, ad, x, y);
 724       CalcError(r, ad, x, y, gridsize, Grid, D, W, E);
 725       int err = Search(r, Ext, gridsize, E);
 726       if (err) return err;
 727       for(int i=0; i <= r; i++) assert(Ext[i]<gridsize);
 728       if (isDone(r, Ext, E))
 729          break;
 730    }
 731
 732    CalcParms(r, Ext, Grid, D, W, ad, x, y);
 733
 734 /*
 735  * Find the 'taps' of the filter for use with Frequency
 736  * Sampling.  If odd or Negative symmetry, fix the taps
 737  * according to Parks McClellan
 738  */
 739    for (i=0; i<=numtaps/2; i++)
 740    {
 741       if (symmetry == POSITIVE)
 742       {
 743          if (numtaps%2)
 744             c = 1;
 745          else
 746             c = cos(Pi * (double)i/numtaps);
 747       }
 748       else
 749       {
 750          if (numtaps%2)
 751             c = sin(Pi2 * (double)i/numtaps);
 752          else
 753             c = sin(Pi * (double)i/numtaps);
 754       }
 755       taps[i] = ComputeA((double)i/numtaps, r, ad, x, y)*c;
 756    }
 757
 758 /*
 759  * Frequency sampling design with calculated taps
 760  */
 761    FreqSample(numtaps, taps, h, symmetry);
 762
 763 /*
 764  * Delete allocated memory
 765  */
 766    free(Grid);
 767    free(W);
 768    free(D);
 769    free(E);
 770    free(Ext);
 771    free(x);
 772    free(y);
 773    free(ad);
 774    return iter<MAXITERATIONS?0:-1;
 775 }
 776
 777 //////////////////////////////////////////////////////////////////////////////
 778 //
 779 //                          GNU Radio interface
 780 //
 781 //////////////////////////////////////////////////////////////////////////////
 782
 783
 784 static void
 785 punt (const std::string msg)
 786 {
 787   std::cerr << msg << '\n';
 788   throw std::runtime_error (msg);
 789 }
 790
 791 std::vector<double>
 792 gr_remez (int order,
 793           const std::vector<double> &arg_bands,
 794           const std::vector<double> &arg_response,
 795           const std::vector<double> &arg_weight,
 796           const std::string filter_type,
 797           int grid_density
 798           ) throw (std::runtime_error)
 799 {
 800   int numtaps = order + 1;
 801   if (numtaps < 4)
 802     punt ("gr_remez: number of taps must be >= 3");
 803
 804   int numbands = arg_bands.size () / 2;
 805   LOCAL_BUFFER (double, bands, numbands * 2);
 806   if (numbands < 1 || arg_bands.size () % 2 == 1)
 807     punt ("gr_remez: must have an even number of band edges");
 808
 809   for (unsigned int i = 1; i < arg_bands.size (); i++){
 810     if (arg_bands[i] < arg_bands[i-1])
 811       punt ("gr_remez: band edges must be nondecreasing");
 812   }
 813
 814   if (arg_bands[0] < 0 || arg_bands[arg_bands.size () - 1] > 1)
 815     punt ("gr_remez: band edges must be in the range [0,1]");
 816
 817   // Divide by 2 to fit with the implementation that uses a
 818   // sample rate of [0, 0.5] instead of [0, 1.0]
 819   for (int i = 0; i < 2 * numbands; i++)
 820     bands[i] = arg_bands[i] / 2;
 821
 822   LOCAL_BUFFER (double, response, numbands * 2);
 823   if (arg_response.size () != arg_bands.size ())
 824     punt ("gr_remez: must have one response magnitude for each band edge");
 825
 826   for (int i = 0; i < 2 * numbands; i++)
 827     response[i] = arg_response[i];
 828
 829   LOCAL_BUFFER (double, weight, numbands);
 830   for (int i = 0; i < numbands; i++)
 831     weight[i] = 1.0;
 832
 833   if (arg_weight.size () != 0){
 834     if ((int) arg_weight.size () != numbands)
 835       punt ("gr_remez: need one weight for each band [=length(band)/2]");
 836     for (int i = 0; i < numbands; i++)
 837       weight[i] = arg_weight [i];
 838   }
 839
 840   int itype = 0;
 841   if (filter_type == "bandpass")
 842     itype = BANDPASS;
 843   else if (filter_type == "differentiator")
 844     itype = DIFFERENTIATOR;
 845   else if (filter_type == "hilbert")
 846     itype = HILBERT;
 847   else
 848     punt ("gr_remez: unknown ftype '" + filter_type + "'");
 849
 850   if (grid_density < 16)
 851     punt ("gr_remez: grid_density is too low; must be >= 16");
 852
 853   LOCAL_BUFFER (double, coeff, numtaps + 5);            // FIXME why + 5?
 854   int err = remez (coeff, numtaps, numbands,
 855                    bands, response, weight, itype, grid_density);
 856
 857   if (err == -1)
 858     punt ("gr_remez: failed to converge");
 859
 860   if (err == -2)
 861     punt ("gr_remez: insufficient extremals -- cannot continue");
 862
 863   if (err == -3)
 864     punt ("gr_remez: too many extremals -- cannot continue");
 865
 866   return std::vector<double> (&coeff[0], &coeff[numtaps]);
 867 }
 868
 869
 870
 871 #if 0
 872 /* == Octave interface starts here ====================================== */
 873
 874 DEFUN_DLD (remez, args, ,
 875   "b = remez(n, f, a [, w] [, ftype] [, griddensity])\n\
 876 Parks-McClellan optimal FIR filter design.\n\
 877 n gives the number of taps in the returned filter\n\
 878 f gives frequency at the band edges [ b1 e1 b2 e2 b3 e3 ...]\n\
 879 a gives amplitude at the band edges [ a(b1) a(e1) a(b2) a(e2) ...]\n\
 880 w gives weighting applied to each band\n\
 881 ftype is 'bandpass', 'hilbert' or 'differentiator'\n\
 882 griddensity determines how accurately the filter will be\n\
 883     constructed. The minimum value is 16, but higher numbers are\n\
 884     slower to compute.\n\
 885 \n\
 886 Frequency is in the range (0, 1), with 1 being the nyquist frequency")
 887 {
 888   octave_value_list retval;
 889   int i;
 890
 891   int nargin = args.length();
 892   if (nargin < 3 || nargin > 6) {
 893     print_usage("remez");
 894     return retval;
 895   }
 896
 897   int numtaps = NINT (args(0).double_value()) + 1; // #coeff = filter order+1
 898   if (numtaps < 4) {
 899     error("remez: number of taps must be an integer greater than 3");
 900     return retval;
 901   }
 902
 903   ColumnVector o_bands(args(1).vector_value());
 904   int numbands = o_bands.length()/2;
 905   OCTAVE_LOCAL_BUFFER(double, bands, numbands*2);
 906   if (numbands < 1 || o_bands.length()%2 == 1) {
 907     error("remez: must have an even number of band edges");
 908     return retval;
 909   }
 910   for (i=1; i < o_bands.length(); i++) {
 911     if (o_bands(i)<o_bands(i-1)) {
 912       error("band edges must be nondecreasing");
 913       return retval;
 914     }
 915   }
 916   if (o_bands(0) < 0 || o_bands(1) > 1) {
 917     error("band edges must be in the range [0,1]");
 918     return retval;
 919   }
 920   for(i=0; i < 2*numbands; i++) bands[i] = o_bands(i)/2.0;
 921
 922   ColumnVector o_response(args(2).vector_value());
 923   OCTAVE_LOCAL_BUFFER (double, response, numbands*2);
 924   if (o_response.length() != o_bands.length()) {
 925     error("remez: must have one response magnitude for each band edge");
 926     return retval;
 927   }
 928   for(i=0; i < 2*numbands; i++) response[i] = o_response(i);
 929
 930   std::string stype = std::string("bandpass");
 931   int density = 16;
 932   OCTAVE_LOCAL_BUFFER (double, weight, numbands);
 933   for (i=0; i < numbands; i++) weight[i] = 1.0;
 934   if (nargin > 3) {
 935     if (args(3).is_real_matrix()) {
 936       ColumnVector o_weight(args(3).vector_value());
 937       if (o_weight.length() != numbands) {
 938         error("remez: need one weight for each band [=length(band)/2]");
 939         return retval;
 940       }
 941       for (i=0; i < numbands; i++) weight[i] = o_weight(i);
 942     }
 943     else if (args(3).is_string())
 944       stype = args(3).string_value();
 945     else if (args(3).is_real_scalar())
 946       density = NINT(args(3).double_value());
 947     else {
 948       error("remez: incorrect argument list");
 949       return retval;
 950     }
 951   }
 952   if (nargin > 4) {
 953     if (args(4).is_string() && !args(3).is_string())
 954       stype = args(4).string_value();
 955     else if (args(4).is_real_scalar() && !args(3).is_real_scalar())
 956       density = NINT(args(4).double_value());
 957     else {
 958       error("remez: incorrect argument list");
 959       return retval;
 960     }
 961   }
 962   if (nargin > 5) {
 963     if (args(5).is_real_scalar()
 964         && !args(4).is_real_scalar()
 965         && !args(3).is_real_scalar())
 966       density = NINT(args(4).double_value());
 967     else {
 968       error("remez: incorrect argument list");
 969       return retval;
 970     }
 971   }
 972
 973   int itype;
 974   if (stype == "bandpass")
 975     itype = BANDPASS;
 976   else if (stype == "differentiator")
 977     itype = DIFFERENTIATOR;
 978   else if (stype == "hilbert")
 979     itype = HILBERT;
 980   else {
 981     error("remez: unknown ftype '%s'", stype.data());
 982     return retval;
 983   }
 984
 985   if (density < 16) {
 986     error("remez: griddensity is too low; must be greater than 16");
 987     return retval;
 988   }
 989
 990   OCTAVE_LOCAL_BUFFER (double, coeff, numtaps+5);
 991   int err = remez(coeff,numtaps,numbands,bands,response,weight,itype,density);
 992
 993   if (err == -1)
 994     warning("remez: -- failed to converge -- returned filter may be bad.");
 995   else if (err == -2) {
 996     error("remez: insufficient extremals--cannot continue");
 997     return retval;
 998   }
 999   else if (err == -3) {
1000     error("remez: too many extremals--cannot continue");
1001     return retval;
1002   }
1003
1004   ColumnVector h(numtaps);
1005   while(numtaps--) h(numtaps) = coeff[numtaps];
1006
1007   return octave_value(h);
1008 }
1009
1010 /*
1011 %!test
1012 %! b = [
1013 %!    0.0415131831103279
1014 %!    0.0581639884202646
1015 %!   -0.0281579212691008
1016 %!   -0.0535575358002337
1017 %!   -0.0617245915143180
1018 %!    0.0507753178978075
1019 %!    0.2079018331396460
1020 %!    0.3327160895375440
1021 %!    0.3327160895375440
1022 %!    0.2079018331396460
1023 %!    0.0507753178978075
1024 %!   -0.0617245915143180
1025 %!   -0.0535575358002337
1026 %!   -0.0281579212691008
1027 %!    0.0581639884202646
1028 %!    0.0415131831103279];
1029 %! assert(remez(15,[0,0.3,0.4,1],[1,1,0,0]),b,1e-14);
1030
1031  */
1032
1033 #endif