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   for (int i = 0; i < 2 * numbands; i++)
 818     bands[i] = arg_bands[i] / 2;                // FIXME why / 2?
 819
 820   LOCAL_BUFFER (double, response, numbands * 2);
 821   if (arg_response.size () != arg_bands.size ())
 822     punt ("gr_remez: must have one response magnitude for each band edge");
 823
 824   for (int i = 0; i < 2 * numbands; i++)
 825     response[i] = arg_response[i];
 826
 827   LOCAL_BUFFER (double, weight, numbands);
 828   for (int i = 0; i < numbands; i++)
 829     weight[i] = 1.0;
 830
 831   if (arg_weight.size () != 0){
 832     if ((int) arg_weight.size () != numbands)
 833       punt ("gr_remez: need one weight for each band [=length(band)/2]");
 834     for (int i = 0; i < numbands; i++)
 835       weight[i] = arg_weight [i];
 836   }
 837
 838   int itype = 0;
 839   if (filter_type == "bandpass")
 840     itype = BANDPASS;
 841   else if (filter_type == "differentiator")
 842     itype = DIFFERENTIATOR;
 843   else if (filter_type == "hilbert")
 844     itype = HILBERT;
 845   else
 846     punt ("gr_remez: unknown ftype '" + filter_type + "'");
 847
 848   if (grid_density < 16)
 849     punt ("gr_remez: grid_density is too low; must be >= 16");
 850
 851   LOCAL_BUFFER (double, coeff, numtaps + 5);            // FIXME why + 5?
 852   int err = remez (coeff, numtaps, numbands,
 853                    bands, response, weight, itype, grid_density);
 854
 855   if (err == -1)
 856     punt ("gr_remez: failed to converge");
 857
 858   if (err == -2)
 859     punt ("gr_remez: insufficient extremals -- cannot continue");
 860
 861   if (err == -3)
 862     punt ("gr_remez: too many extremals -- cannot continue");
 863
 864   return std::vector<double> (&coeff[0], &coeff[numtaps]);
 865 }
 866
 867
 868
 869 #if 0
 870 /* == Octave interface starts here ====================================== */
 871
 872 DEFUN_DLD (remez, args, ,
 873   "b = remez(n, f, a [, w] [, ftype] [, griddensity])\n\
 874 Parks-McClellan optimal FIR filter design.\n\
 875 n gives the number of taps in the returned filter\n\
 876 f gives frequency at the band edges [ b1 e1 b2 e2 b3 e3 ...]\n\
 877 a gives amplitude at the band edges [ a(b1) a(e1) a(b2) a(e2) ...]\n\
 878 w gives weighting applied to each band\n\
 879 ftype is 'bandpass', 'hilbert' or 'differentiator'\n\
 880 griddensity determines how accurately the filter will be\n\
 881     constructed. The minimum value is 16, but higher numbers are\n\
 882     slower to compute.\n\
 883 \n\
 884 Frequency is in the range (0, 1), with 1 being the nyquist frequency")
 885 {
 886   octave_value_list retval;
 887   int i;
 888
 889   int nargin = args.length();
 890   if (nargin < 3 || nargin > 6) {
 891     print_usage("remez");
 892     return retval;
 893   }
 894
 895   int numtaps = NINT (args(0).double_value()) + 1; // #coeff = filter order+1
 896   if (numtaps < 4) {
 897     error("remez: number of taps must be an integer greater than 3");
 898     return retval;
 899   }
 900
 901   ColumnVector o_bands(args(1).vector_value());
 902   int numbands = o_bands.length()/2;
 903   OCTAVE_LOCAL_BUFFER(double, bands, numbands*2);
 904   if (numbands < 1 || o_bands.length()%2 == 1) {
 905     error("remez: must have an even number of band edges");
 906     return retval;
 907   }
 908   for (i=1; i < o_bands.length(); i++) {
 909     if (o_bands(i)<o_bands(i-1)) {
 910       error("band edges must be nondecreasing");
 911       return retval;
 912     }
 913   }
 914   if (o_bands(0) < 0 || o_bands(1) > 1) {
 915     error("band edges must be in the range [0,1]");
 916     return retval;
 917   }
 918   for(i=0; i < 2*numbands; i++) bands[i] = o_bands(i)/2.0;
 919
 920   ColumnVector o_response(args(2).vector_value());
 921   OCTAVE_LOCAL_BUFFER (double, response, numbands*2);
 922   if (o_response.length() != o_bands.length()) {
 923     error("remez: must have one response magnitude for each band edge");
 924     return retval;
 925   }
 926   for(i=0; i < 2*numbands; i++) response[i] = o_response(i);
 927
 928   std::string stype = std::string("bandpass");
 929   int density = 16;
 930   OCTAVE_LOCAL_BUFFER (double, weight, numbands);
 931   for (i=0; i < numbands; i++) weight[i] = 1.0;
 932   if (nargin > 3) {
 933     if (args(3).is_real_matrix()) {
 934       ColumnVector o_weight(args(3).vector_value());
 935       if (o_weight.length() != numbands) {
 936         error("remez: need one weight for each band [=length(band)/2]");
 937         return retval;
 938       }
 939       for (i=0; i < numbands; i++) weight[i] = o_weight(i);
 940     }
 941     else if (args(3).is_string())
 942       stype = args(3).string_value();
 943     else if (args(3).is_real_scalar())
 944       density = NINT(args(3).double_value());
 945     else {
 946       error("remez: incorrect argument list");
 947       return retval;
 948     }
 949   }
 950   if (nargin > 4) {
 951     if (args(4).is_string() && !args(3).is_string())
 952       stype = args(4).string_value();
 953     else if (args(4).is_real_scalar() && !args(3).is_real_scalar())
 954       density = NINT(args(4).double_value());
 955     else {
 956       error("remez: incorrect argument list");
 957       return retval;
 958     }
 959   }
 960   if (nargin > 5) {
 961     if (args(5).is_real_scalar()
 962         && !args(4).is_real_scalar()
 963         && !args(3).is_real_scalar())
 964       density = NINT(args(4).double_value());
 965     else {
 966       error("remez: incorrect argument list");
 967       return retval;
 968     }
 969   }
 970
 971   int itype;
 972   if (stype == "bandpass")
 973     itype = BANDPASS;
 974   else if (stype == "differentiator")
 975     itype = DIFFERENTIATOR;
 976   else if (stype == "hilbert")
 977     itype = HILBERT;
 978   else {
 979     error("remez: unknown ftype '%s'", stype.data());
 980     return retval;
 981   }
 982
 983   if (density < 16) {
 984     error("remez: griddensity is too low; must be greater than 16");
 985     return retval;
 986   }
 987
 988   OCTAVE_LOCAL_BUFFER (double, coeff, numtaps+5);
 989   int err = remez(coeff,numtaps,numbands,bands,response,weight,itype,density);
 990
 991   if (err == -1)
 992     warning("remez: -- failed to converge -- returned filter may be bad.");
 993   else if (err == -2) {
 994     error("remez: insufficient extremals--cannot continue");
 995     return retval;
 996   }
 997   else if (err == -3) {
 998     error("remez: too many extremals--cannot continue");
 999     return retval;
1000   }
1001
1002   ColumnVector h(numtaps);
1003   while(numtaps--) h(numtaps) = coeff[numtaps];
1004
1005   return octave_value(h);
1006 }
1007
1008 /*
1009 %!test
1010 %! b = [
1011 %!    0.0415131831103279
1012 %!    0.0581639884202646
1013 %!   -0.0281579212691008
1014 %!   -0.0535575358002337
1015 %!   -0.0617245915143180
1016 %!    0.0507753178978075
1017 %!    0.2079018331396460
1018 %!    0.3327160895375440
1019 %!    0.3327160895375440
1020 %!    0.2079018331396460
1021 %!    0.0507753178978075
1022 %!   -0.0617245915143180
1023 %!   -0.0535575358002337
1024 %!   -0.0281579212691008
1025 %!    0.0581639884202646
1026 %!    0.0415131831103279];
1027 %! assert(remez(15,[0,0.3,0.4,1],[1,1,0,0]),b,1e-14);
1028
1029  */
1030
1031 #endif