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.
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.
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.
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
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 *************************************************************************/
36 * This code was extracted from octave.sf.net, and wrapped with
51 #define LOCAL_BUFFER(T, buf, size) \
52 std::vector<T> buf ## _vector (size); \
53 T *buf = &(buf ## _vector[0])
59 #define DIFFERENTIATOR 2
65 #define Pi 3.14159265358979323846
68 #define GRIDDENSITY 16
69 #define MAXITERATIONS 40
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
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
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]
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)
105 double delf, lowf, highf, grid0;
107 delf = 0.5/(griddensity*r);
110 * For differentiator, hilbert,
111 * symmetry is odd and Grid[0] = max(delf, bands[0])
113 grid0 = (symmetry == NEGATIVE) && (delf > bands[0]) ? delf : bands[0];
116 for (band=0; band < numband; band++)
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 */
123 D[j] = des[2*band] + i*(des[2*band+1]-des[2*band])/(k-1);
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
136 if ((symmetry == NEGATIVE) &&
137 (Grid[gridsize-1] > (0.5 - delf)) &&
140 Grid[gridsize-1] = 0.5-delf;
145 /********************
148 * Places Extremal Frequencies evenly throughout the dense grid.
153 * int r - 1/2 the number of filter coefficients
154 * int gridsize - Number of elements in the dense frequency grid
158 * int Ext[] - Extremal indexes to dense frequency grid [r+1]
159 ********************/
162 InitialGuess (int r, int Ext[], int gridsize)
167 Ext[i] = i * (gridsize-1) / r;
171 /***********************
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]
186 * double ad[] - 'b' in Oppenheim & Schafer [r+1]
188 * double y[] - 'C' in Oppenheim & Schafer [r+1]
189 ***********************/
192 CalcParms (int r, int Ext[], double Grid[], double D[], double W[],
193 double ad[], double x[], double y[])
196 double sign, xi, delta, denom, numer;
202 x[i] = cos(Pi2 * Grid[Ext[i]]);
205 * Calculate ad[] - Oppenheim & Schafer eq 7.132
207 ld = (r-1)/15 + 1; /* Skips around to avoid round errors */
214 for (k=j; k<=r; k+=ld)
216 denom *= 2.0*(xi - x[k]);
218 if (fabs(denom)<0.00001)
224 * Calculate delta - Oppenheim & Schafer eq 7.131
230 numer += ad[i] * D[Ext[i]];
231 denom += sign * ad[i]/W[Ext[i]];
238 * Calculate y[] - Oppenheim & Schafer eq 7.133b
242 y[i] = D[Ext[i]] - sign * delta/W[Ext[i]];
248 /*********************
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.
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]
262 * double y[] - 'C' in Oppenheim & Schafer [r+1]
266 * Returns double value of A[freq]
267 *********************/
270 ComputeA (double freq, int r, double ad[], double x[], double y[])
273 double xc, c, denom, numer;
276 xc = cos(Pi2 * freq);
280 if (fabs(c) < 1.0e-7)
294 /************************
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[])
304 * int r - 1/2 the number of filter coefficients
305 * double ad[] - [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]
315 * double E[] - Error function on dense grid [gridsize]
316 ************************/
319 CalcError (int r, double ad[], double x[], double y[],
320 int gridsize, double Grid[],
321 double D[], double W[], double E[])
326 for (i=0; i<gridsize; i++)
328 A = ComputeA(Grid[i], r, ad, x, y);
329 E[i] = W[i] * (D[i] - A);
333 /************************
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
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
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]
355 * int Ext[] - New indexes to extremal frequencies [r+1]
356 ************************/
358 Search (int r, int Ext[],
359 int gridsize, double E[])
361 int i, j, k, l, extra; /* Counters */
363 int *foundExt; /* Array of found extremals */
366 * Allocate enough space for found extremals.
368 foundExt = (int *)malloc((2*r) * sizeof(int));
372 * Check for extremum at 0.
374 if (((E[0]>0.0) && (E[0]>E[1])) ||
375 ((E[0]<0.0) && (E[0]<E[1])))
379 * Check for extrema inside dense grid
381 for (i=1; i<gridsize-1; i++)
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;
392 * Check for extremum at 0.5
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;
401 // PAK: we sometimes get not enough extremal frequencies
402 if (k < r+1) return -2;
406 * Remove extra extremals
413 if (E[foundExt[0]] > 0.0)
414 up = 1; /* first one is a maxima */
416 up = 0; /* first one is a minima */
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 */
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:
435 // if (fabs(E[foundExt[j]]) < fabs(E[foundExt[j-1]])) l=j;
437 break; /* Ooops, found two non-alternating */
438 } /* extrema. Delete smallest of them */
439 } /* if the loop finishes, all extrema are alternating */
442 * If there's only one extremal and all are alternating,
443 * delete the smallest of the first/last extremals.
445 if ((alt) && (extra == 1))
447 if (fabs(E[foundExt[k-1]]) < fabs(E[foundExt[0]]))
448 /* Delete last extremal */
450 // PAK: changed from l = foundExt[k-1];
452 /* Delete first extremal */
454 // PAK: changed from l = foundExt[0];
457 for (j=l; j<k-1; j++) /* Loop that does the deletion */
459 foundExt[j] = foundExt[j+1];
460 assert(foundExt[j]<gridsize);
468 assert(foundExt[i]<gridsize);
469 Ext[i] = foundExt[i]; /* Copy found extremals to Ext[] */
477 /*********************
480 * Simple frequency sampling algorithm to determine the impulse
481 * response h[] from A's found in ComputeA
486 * int N - Number of filter coefficients
487 * double A[] - Sample points of desired response [N/2]
488 * int symmetry - Symmetry of desired filter
492 * double h[] - Impulse Response of final filter [N]
493 *********************/
495 FreqSample (int N, double A[], double h[], int symm)
501 if (symm == POSITIVE)
510 val += 2.0 * A[k] * cos(x*k);
520 for (k=1; k<=(N/2-1); k++)
521 val += 2.0 * A[k] * cos(x*k);
535 val += 2.0 * A[k] * sin(x*k);
543 val = A[N/2] * sin(Pi * (n - M));
545 for (k=1; k<=(N/2-1); k++)
546 val += 2.0 * A[k] * sin(x*k);
556 * Checks to see if the error function is small enough to consider
557 * the result to have converged.
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]
567 * Returns 1 if the result converged
568 * Returns 0 if the result has not converged
569 ********************/
572 isDone (int r, int Ext[], double E[])
575 double min, max, current;
577 min = max = fabs(E[Ext[0]]);
580 current = fabs(E[Ext[i]]);
586 return (((max-min)/max) < 0.0001);
589 /********************
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.
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
608 * double h[] - Impulse response of final filter [numtaps]
609 * returns - true on success, false on failure to converge
610 ********************/
613 remez (double h[], int numtaps,
614 int numband, const double bands[],
615 const double des[], const double weight[],
616 int type, int griddensity)
618 double *Grid, *W, *D, *E;
619 int i, iter, gridsize, r, *Ext;
624 if (type == BANDPASS)
629 r = numtaps/2; /* number of extrema */
630 if ((numtaps%2) && (symmetry == POSITIVE))
634 * Predict dense grid size in advance for memory allocation
635 * .5 is so we round up, not truncate
638 for (i=0; i<numband; i++)
640 gridsize += (int)(2*r*griddensity*(bands[2*i+1] - bands[2*i]) + .5);
642 if (symmetry == NEGATIVE)
648 * Dynamically allocate memory for arrays with proper sizes
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));
661 * Create dense frequency grid
663 CreateDenseGrid(r, numtaps, numband, bands, des, weight,
664 gridsize, Grid, D, W, symmetry, griddensity);
665 InitialGuess(r, Ext, gridsize);
668 * For Differentiator: (fix grid)
670 if (type == DIFFERENTIATOR)
672 for (i=0; i<gridsize; i++)
674 /* D[i] = D[i]*Grid[i]; */
681 * For odd or Negative symmetry filters, alter the
682 * D[] and W[] according to Parks McClellan
684 if (symmetry == POSITIVE)
686 if (numtaps % 2 == 0)
688 for (i=0; i<gridsize; i++)
690 c = cos(Pi * Grid[i]);
700 for (i=0; i<gridsize; i++)
702 c = sin(Pi2 * Grid[i]);
709 for (i=0; i<gridsize; i++)
711 c = sin(Pi * Grid[i]);
719 * Perform the Remez Exchange algorithm
721 for (iter=0; iter<MAXITERATIONS; iter++)
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);
727 for(int i=0; i <= r; i++) assert(Ext[i]<gridsize);
728 if (isDone(r, Ext, E))
732 CalcParms(r, Ext, Grid, D, W, ad, x, y);
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
739 for (i=0; i<=numtaps/2; i++)
741 if (symmetry == POSITIVE)
746 c = cos(Pi * (double)i/numtaps);
751 c = sin(Pi2 * (double)i/numtaps);
753 c = sin(Pi * (double)i/numtaps);
755 taps[i] = ComputeA((double)i/numtaps, r, ad, x, y)*c;
759 * Frequency sampling design with calculated taps
761 FreqSample(numtaps, taps, h, symmetry);
764 * Delete allocated memory
774 return iter<MAXITERATIONS?0:-1;
777 //////////////////////////////////////////////////////////////////////////////
779 // GNU Radio interface
781 //////////////////////////////////////////////////////////////////////////////
785 punt (const std::string msg)
787 std::cerr << msg << '\n';
788 throw std::runtime_error (msg);
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,
798 ) throw (std::runtime_error)
800 int numtaps = order + 1;
802 punt ("gr_remez: number of taps must be >= 3");
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");
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");
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]");
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;
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");
826 for (int i = 0; i < 2 * numbands; i++)
827 response[i] = arg_response[i];
829 LOCAL_BUFFER (double, weight, numbands);
830 for (int i = 0; i < numbands; i++)
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];
841 if (filter_type == "bandpass")
843 else if (filter_type == "differentiator")
844 itype = DIFFERENTIATOR;
845 else if (filter_type == "hilbert")
848 punt ("gr_remez: unknown ftype '" + filter_type + "'");
850 if (grid_density < 16)
851 punt ("gr_remez: grid_density is too low; must be >= 16");
853 LOCAL_BUFFER (double, coeff, numtaps + 5); // FIXME why + 5?
854 int err = remez (coeff, numtaps, numbands,
855 bands, response, weight, itype, grid_density);
858 punt ("gr_remez: failed to converge");
861 punt ("gr_remez: insufficient extremals -- cannot continue");
864 punt ("gr_remez: too many extremals -- cannot continue");
866 return std::vector<double> (&coeff[0], &coeff[numtaps]);
872 /* == Octave interface starts here ====================================== */
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\
886 Frequency is in the range (0, 1), with 1 being the nyquist frequency")
888 octave_value_list retval;
891 int nargin = args.length();
892 if (nargin < 3 || nargin > 6) {
893 print_usage("remez");
897 int numtaps = NINT (args(0).double_value()) + 1; // #coeff = filter order+1
899 error("remez: number of taps must be an integer greater than 3");
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");
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");
916 if (o_bands(0) < 0 || o_bands(1) > 1) {
917 error("band edges must be in the range [0,1]");
920 for(i=0; i < 2*numbands; i++) bands[i] = o_bands(i)/2.0;
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");
928 for(i=0; i < 2*numbands; i++) response[i] = o_response(i);
930 std::string stype = std::string("bandpass");
932 OCTAVE_LOCAL_BUFFER (double, weight, numbands);
933 for (i=0; i < numbands; i++) weight[i] = 1.0;
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]");
941 for (i=0; i < numbands; i++) weight[i] = o_weight(i);
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());
948 error("remez: incorrect argument list");
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());
958 error("remez: incorrect argument list");
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());
968 error("remez: incorrect argument list");
974 if (stype == "bandpass")
976 else if (stype == "differentiator")
977 itype = DIFFERENTIATOR;
978 else if (stype == "hilbert")
981 error("remez: unknown ftype '%s'", stype.data());
986 error("remez: griddensity is too low; must be greater than 16");
990 OCTAVE_LOCAL_BUFFER (double, coeff, numtaps+5);
991 int err = remez(coeff,numtaps,numbands,bands,response,weight,itype,density);
994 warning("remez: -- failed to converge -- returned filter may be bad.");
995 else if (err == -2) {
996 error("remez: insufficient extremals--cannot continue");
999 else if (err == -3) {
1000 error("remez: too many extremals--cannot continue");
1004 ColumnVector h(numtaps);
1005 while(numtaps--) h(numtaps) = coeff[numtaps];
1007 return octave_value(h);
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);