CSL 5.2

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

FIR.cpp
Go to the documentation of this file.
1 //
2 // FIR.cpp -- CSL FIR filter class
3 // See the copyright notice and acknowledgment of authors in the file COPYRIGHT
4 //
5
6 #include "FIR.h"
7 #include <stdlib.h>
8 #include <math.h>
9 #include <string.h>
10
11 using namespace csl;
12
13 // FilterSpecification class
14
15 // Constants used by the design algorithm
16
17 #define BANDPASS 1
18 #define DIFFERENTIATOR 2
19 #define HILBERT 3
20
21 // Constructors
22
23 FilterSpecification::FilterSpecification (unsigned numTaps, unsigned numBands, double *freqs, double *resps, double *weights)
24 : mNumTaps(numTaps), mNumBands(numBands), mFrequencies(NULL), mResponses(NULL), mWeights(NULL), mTapData(NULL) {
25 setNumTaps(numTaps);
26 setFrequencies(freqs);
27 setResponses(resps);
28 setWeights(weights);
29 if (freqs != NULL && resps != NULL && weights != NULL)
30 planFilter();
31 }
32
33 FilterSpecification::~FilterSpecification () {
34 delete[] mFrequencies;
35 delete[] mResponses;
36 delete[] mWeights;
37 delete[] mTapData;
38 }
39
40 // Accessors
41
42 void FilterSpecification::setFrequencies(double *frequencies) {
43 if (mFrequencies == NULL)
44 mFrequencies = new double[mNumBands * 2];
45 if (frequencies != NULL)
46 memcpy(mFrequencies, frequencies, (mNumBands * 2 * sizeof(double)));
47 }
48
49 void FilterSpecification::setResponses(double *responses) {
50 if (mResponses == NULL)
51 mResponses = new double[mNumBands];
52 if (responses != NULL)
53 memcpy(mResponses, responses, (mNumBands * sizeof(double)));
54 }
55
56 void FilterSpecification::setWeights(double *weights) {
57 if (mWeights == NULL)
58 mWeights = new double[mNumBands];
59 if (weights != NULL)
60 memcpy(mWeights, weights, (mNumBands * sizeof(double)));
61 }
62
63 void FilterSpecification::setNumTaps(unsigned numTaps) {
64 if (numTaps < 128) {
65 mNumTaps = numTaps;
66 if (mTapData != NULL)
67 delete[] mTapData;
68 mTapData = new double[mNumTaps];
69 } else logMsg(kLogError, "Too many taps. The use of convolution is faster if more taps needed.");
70
71 }
72
73 // method to plan the filter (execute the search/iterate algorithm)
74
75 void remez(double h[], int numtaps, int numband, double bands[], double des[], double weight[], int type);
76
77 void FilterSpecification::planFilter() {
78 bool shouldNormalize = false;
79 unsigned i; // check of the frequencies values are normalized (0.0 - 0.5 Fs) or not
80 for (i = 0; i < (mNumBands * 2); i++) {
81 if (mFrequencies[i] > 0.5) {
82 shouldNormalize = true;
83 break;
84 }
85 }
86 if (shouldNormalize) // normalize the frequencies if they're given in Hz
87 for (i = 0; i < (mNumBands * 2); i++)
88 mFrequencies[i] = mFrequencies[i] / CGestalt::frameRate();
89 // now call the iterative design function
90 remez(mTapData, mNumTaps, mNumBands, mFrequencies, mResponses, mWeights, BANDPASS);
91 }
92
93 // FIR class
94
95 FIR::FIR(UnitGenerator & in, unsigned numTaps, float * tapDelays) : Effect(in) {
96 mFilterSpec = new FilterSpecification;
97 setTaps(numTaps, tapDelays);
98
99 if (tapDelays == NULL) { // If passed num taps, but no IR, then just make a cheap lowpass averaging filter.
100 for (unsigned i = 0; i < numTaps; i++)
101 mFilterSpec->mTapData[i] = 0.5;
102 }
103 }
104
105 // Read in a tap data file
106 FIR::FIR(UnitGenerator & in, char * fileName) : Effect(in) {
107 mFilterSpec = new FilterSpecification;
108 readTaps(fileName);
109 }
110
111 // give it a filter specification object
112
113 FIR::FIR (FilterSpecification & fs) : mFilterSpec(&fs) {
114 resetDLine();
115 fs.planFilter();
116 }
117
118 FIR::FIR (UnitGenerator & in, FilterSpecification & fs) : Effect(in), mFilterSpec(&fs) {
119 resetDLine();
120 fs.planFilter();
121 }
122
123 FIR::~FIR () {
124 delete[] mDLine;
125 }
126
127 void FIR::resetDLine() {
128 mDLine = new sample[mFilterSpec->mNumTaps];
129 for (unsigned i = 0; i < mFilterSpec->mNumTaps; i++) // empty the delay line
130 mDLine[i] = 0.0;
131 mOffset = 0;
132 }
133
134 void FIR::setTaps(unsigned numTaps, float *tapDelays) {
135 mFilterSpec->setNumTaps(numTaps);
136 if (tapDelays != NULL)
137 memcpy(mFilterSpec->mTapData, tapDelays, (numTaps * sizeof(double)));
138 resetDLine();
139 mOffset = 0;
140 }
141
142 void FIR::readTaps(char * fileName) {
143 unsigned i, numTaps;
144 FILE * configData;
145
146 configData = fopen(fileName, "r"); // open config file
147 fscanf(configData, "%d\n", &numTaps); // read # of taps
148
149 float * theTapData = new float[numTaps];
150 for (i = 0; i < numTaps; i++) // read tap coefficients
151 fscanf(configData, "%f\n", &theTapData[i]);
152 setTaps(numTaps, theTapData);
153 fclose(configData);
154
155 #ifdef CSL_DEBUG
156 logMsg("Taps:"); // print out the tap data
157 for (i = 0; i < numTaps; i++)
158 logMsg("\t%g\n ", theTapData[i]);
159 #endif
160 resetDLine();
161 }
162
163 void FIR::nextBuffer(Buffer &outputBuffer, unsigned outBufNum) throw (CException) {
164 float sum;
165 int tap;
166 sample * outPtr = outputBuffer.buffer(outBufNum);
167 unsigned numFrames = outputBuffer.mNumFrames;
168 unsigned numTaps = mFilterSpec->mNumTaps;
169 double *tapDataPtr = mFilterSpec->mTapData;
170 #ifdef CSL_DEBUG
171 logMsg("FIR nextBuffer");
172 #endif
173 Effect::pullInput(numFrames); // get some input
174 sample *inputPtr = mInputPtr;
175
176 for (unsigned i = 0; i < numFrames; i++) { // sample block loop
177 sum = 0.0; // output sample accumulator
178 mDLine[mOffset] = inputPtr[i]; // write input sample into delay line
179 tap = mOffset; // initial reader tap position
180 for (unsigned j = 0; j < numTaps; j++) { // loop through delay line taps
181 sum += mDLine[tap] * tapDataPtr[j]; // sum scaled d_line data (previous sample * weighting)
182 tap--; // decrement tap offset
183 if (tap < 0) tap = numTaps - 1; // wrap around
184 } // end of tap-per-sample loop
185 outPtr[i] = sum; // write sum sample to output buffer
186 mOffset++; // increment and wrap offset
187 mOffset %= numTaps; // "%=" means modulo-assignment
188 } // end of sample block loop
189 return;
190 }
191
192 ///////////////////////////////////////////////////////////////////////////////////
193
194 /**************************************************************************
195 * Parks-McClellan algorithm for FIR filter design (C version)
196 *-------------------------------------------------
197 * Copyright (c) 1995,1998 Jake Janovetz (janovetz@uiuc.edu)
198 *
199 * This library is free software; you can redistribute it and/or
200 * modify it under the terms of the GNU Library General Public
201 * License as published by the Free Software Foundation; either
202 * version 2 of the License, or (at your option) any later version.
203 *
204 * This library is distributed in the hope that it will be useful,
205 * but WITHOUT ANY WARRANTY; without even the implied warranty of
206 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
207 * Library General Public License for more details.
208
209 * You should have received a copy of the GNU Library General Public
210 * License along with this library; if not, write to the Free
211 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
212 *
213 *************************************************************************/
214
215 #define NEGATIVE 0
216 #define POSITIVE 1
217
218 #define Pi 3.1415926535897932
219 #define Pi2 6.2831853071795865
220 #define GRIDDENSITY 16
221 #define MAXITERATIONS 40
222
223 /*******************
224 * CreateDenseGrid
225 *=================
226 * Creates the dense grid of frequencies from the specified bands.
227 * Also creates the Desired Frequency Response function (D[]) and
228 * the Weight function (W[]) on that dense grid
229 *
230 *
231 * INPUT:
232 * ------
233 * int r - 1/2 the number of filter coefficients
234 * int numtaps - Number of taps in the resulting filter
235 * int numband - Number of bands in user specification
236 * double bands[] - User-specified band edges [2*numband]
237 * double des[] - Desired response per band [numband]
238 * double weight[] - Weight per band [numband]
239 * int symmetry - Symmetry of filter - used for grid check
240 *
241 * OUTPUT:
242 * -------
243 * int gridsize - Number of elements in the dense frequency grid
244 * double Grid[] - Frequencies (0 to 0.5) on the dense grid [gridsize]
245 * double D[] - Desired response on the dense grid [gridsize]
246 * double W[] - Weight function on the dense grid [gridsize]
247 *******************/
248
249 void CreateDenseGrid(int r, int numtaps, int numband, double bands[],
250 double des[], double weight[], int *gridsize,
251 double Grid[], double D[], double W[],
252 int symmetry)
253 {
254 int i, j, k, band;
255 double delf, lowf, highf;
256
257 delf = 0.5/(GRIDDENSITY*r);
258
259 /*
260 * For differentiator, hilbert,
261 * symmetry is odd and Grid[0] = max(delf, band[0])
262 */
263
264 if ((symmetry == NEGATIVE) && (delf > bands[0]))
265 bands[0] = delf;
266
267 j=0;
268 for (band=0; band < numband; band++)
269 {
270 Grid[j] = bands[2*band];
271 lowf = bands[2*band];
272 highf = bands[2*band + 1];
273 k = (int)((highf - lowf)/delf + 0.5); /* .5 for rounding */
274 for (i=0; i<k; i++)
275 {
276 D[j] = des[band];
277 W[j] = weight[band];
278 Grid[j] = lowf;
279 lowf += delf;
280 j++;
281 }
282 Grid[j-1] = highf;
283 }
284
285 /*
286 * Similar to above, if odd symmetry, last grid point can't be .5
287 * - but, if there are even taps, leave the last grid point at .5
288 */
289 if ((symmetry == NEGATIVE) &&
290 (Grid[*gridsize-1] > (0.5 - delf)) &&
291 (numtaps % 2))
292 {
293 Grid[*gridsize-1] = 0.5-delf;
294 }
295 }
296
297
298 /********************
299 * InitialGuess
300 *==============
301 * Places Extremal Frequencies evenly throughout the dense grid.
302 *
303 *
304 * INPUT:
305 * ------
306 * int r - 1/2 the number of filter coefficients
307 * int gridsize - Number of elements in the dense frequency grid
308 *
309 * OUTPUT:
310 * -------
311 * int Ext[] - Extremal indexes to dense frequency grid [r+1]
312 ********************/
313
314 void InitialGuess(int r, int Ext[], int gridsize)
315 {
316 int i;
317
318 for (i=0; i<=r; i++)
319 Ext[i] = i * (gridsize-1) / r;
320 }
321
322
323 /***********************
324 * CalcParms
325 *===========
326 *
327 *
328 * INPUT:
329 * ------
330 * int r - 1/2 the number of filter coefficients
331 * int Ext[] - Extremal indexes to dense frequency grid [r+1]
332 * double Grid[] - Frequencies (0 to 0.5) on the dense grid [gridsize]
333 * double D[] - Desired response on the dense grid [gridsize]
334 * double W[] - Weight function on the dense grid [gridsize]
335 *
336 * OUTPUT:
337 * -------
338 * double ad[] - 'b' in Oppenheim & Schafer [r+1]
339 * double x[] - [r+1]
340 * double y[] - 'C' in Oppenheim & Schafer [r+1]
341 ***********************/
342
343 void CalcParms(int r, int Ext[], double Grid[], double D[], double W[],
344 double ad[], double x[], double y[])
345 {
346 int i, j, k, ld;
347 double sign, xi, delta, denom, numer;
348
349 /*
350 * Find x[]
351 */
352 for (i=0; i<=r; i++)
353 x[i] = cos(Pi2 * Grid[Ext[i]]);
354
355 /*
356 * Calculate ad[] - Oppenheim & Schafer eq 7.132
357 */
358 ld = (r-1)/15 + 1; /* Skips around to avoid round errors */
359 for (i=0; i<=r; i++)
360 {
361 denom = 1.0;
362 xi = x[i];
363 for (j=0; j<ld; j++)
364 {
365 for (k=j; k<=r; k+=ld)
366 if (k != i)
367 denom *= 2.0*(xi - x[k]);
368 }
369 if (fabs(denom)<0.00001)
370 denom = 0.00001;
371 ad[i] = 1.0/denom;
372 }
373
374 /*
375 * Calculate delta - Oppenheim & Schafer eq 7.131
376 */
377 numer = denom = 0;
378 sign = 1;
379 for (i=0; i<=r; i++)
380 {
381 numer += ad[i] * D[Ext[i]];
382 denom += sign * ad[i]/W[Ext[i]];
383 sign = -sign;
384 }
385 delta = numer/denom;
386 sign = 1;
387
388 /*
389 * Calculate y[] - Oppenheim & Schafer eq 7.133b
390 */
391 for (i=0; i<=r; i++)
392 {
393 y[i] = D[Ext[i]] - sign * delta/W[Ext[i]];
394 sign = -sign;
395 }
396 }
397
398
399 /*********************
400 * ComputeA
401 *==========
402 * Using values calculated in CalcParms, ComputeA calculates the
403 * actual filter response at a given frequency (freq). Uses
404 * eq 7.133a from Oppenheim & Schafer.
405 *
406 *
407 * INPUT:
408 * ------
409 * double freq - Frequency (0 to 0.5) at which to calculate A
410 * int r - 1/2 the number of filter coefficients
411 * double ad[] - 'b' in Oppenheim & Schafer [r+1]
412 * double x[] - [r+1]
413 * double y[] - 'C' in Oppenheim & Schafer [r+1]
414 *
415 * OUTPUT:
416 * -------
417 * Returns double value of A[freq]
418 *********************/
419
420 double ComputeA(double freq, int r, double ad[], double x[], double y[])
421 {
422 int i;
423 double xc, c, denom, numer;
424
425 denom = numer = 0;
426 xc = cos(Pi2 * freq);
427 for (i=0; i<=r; i++)
428 {
429 c = xc - x[i];
430 if (fabs(c) < 1.0e-7)
431 {
432 numer = y[i];
433 denom = 1;
434 break;
435 }
436 c = ad[i]/c;
437 denom += c;
438 numer += c*y[i];
439 }
440 return numer/denom;
441 }
442
443
444 /************************
445 * CalcError
446 *===========
447 * Calculates the Error function from the desired frequency response
448 * on the dense grid (D[]), the weight function on the dense grid (W[]),
449 * and the present response calculation (A[])
450 *
451 *
452 * INPUT:
453 * ------
454 * int r - 1/2 the number of filter coefficients
455 * double ad[] - [r+1]
456 * double x[] - [r+1]
457 * double y[] - [r+1]
458 * int gridsize - Number of elements in the dense frequency grid
459 * double Grid[] - Frequencies on the dense grid [gridsize]
460 * double D[] - Desired response on the dense grid [gridsize]
461 * double W[] - Weight function on the desnse grid [gridsize]
462 *
463 * OUTPUT:
464 * -------
465 * double E[] - Error function on dense grid [gridsize]
466 ************************/
467
468 void CalcError(int r, double ad[], double x[], double y[],
469 int gridsize, double Grid[],
470 double D[], double W[], double E[])
471 {
472 int i;
473 double A;
474
475 for (i=0; i<gridsize; i++)
476 {
477 A = ComputeA(Grid[i], r, ad, x, y);
478 E[i] = W[i] * (D[i] - A);
479 }
480 }
481
482 /************************
483 * Search
484 *========
485 * Searches for the maxima/minima of the error curve. If more than
486 * r+1 extrema are found, it uses the following heuristic (thanks
487 * Chris Hanson):
488 * 1) Adjacent non-alternating extrema deleted first.
489 * 2) If there are more than one excess extrema, delete the
490 * one with the smallest error. This will create a non-alternation
491 * condition that is fixed by 1).
492 * 3) If there is exactly one excess extremum, delete the smaller
493 * of the first/last extremum
494 *
495 *
496 * INPUT:
497 * ------
498 * int r - 1/2 the number of filter coefficients
499 * int Ext[] - Indexes to Grid[] of extremal frequencies [r+1]
500 * int gridsize - Number of elements in the dense frequency grid
501 * double E[] - Array of error values. [gridsize]
502 * OUTPUT:
503 * -------
504 * int Ext[] - New indexes to extremal frequencies [r+1]
505 ************************/
506
507 void Search(int r, int Ext[],
508 int gridsize, double E[])
509 {
510 int i, j, k, l, extra; /* Counters */
511 int up, alt;
512 int *foundExt; /* Array of found extremals */
513
514 /*
515 * Allocate enough space for found extremals.
516 */
517 foundExt = (int *)malloc((2*r) * sizeof(int));
518 k = 0;
519
520 /*
521 * Check for extremum at 0.
522 */
523 if (((E[0]>0.0) && (E[0]>E[1])) ||
524 ((E[0]<0.0) && (E[0]<E[1])))
525 foundExt[k++] = 0;
526
527 /*
528 * Check for extrema inside dense grid
529 */
530 for (i=1; i<gridsize-1; i++)
531 {
532 if (((E[i]>=E[i-1]) && (E[i]>E[i+1]) && (E[i]>0.0)) ||
533 ((E[i]<=E[i-1]) && (E[i]<E[i+1]) && (E[i]<0.0)))
534 foundExt[k++] = i;
535 }
536
537 /*
538 * Check for extremum at 0.5
539 */
540 j = gridsize-1;
541 if (((E[j]>0.0) && (E[j]>E[j-1])) ||
542 ((E[j]<0.0) && (E[j]<E[j-1])))
543 foundExt[k++] = j;
544
545
546 /*
547 * Remove extra extremals
548 */
549 extra = k - (r+1);
550
551 while (extra > 0)
552 {
553 if (E[foundExt[0]] > 0.0)
554 up = 1; /* first one is a maxima */
555 else
556 up = 0; /* first one is a minima */
557
558 l=0;
559 alt = 1;
560 for (j=1; j<k; j++)
561 {
562 if (fabs(E[foundExt[j]]) < fabs(E[foundExt[l]]))
563 l = j; /* new smallest error. */
564 if ((up) && (E[foundExt[j]] < 0.0))
565 up = 0; /* switch to a minima */
566 else if (( ! up) && (E[foundExt[j]] > 0.0))
567 up = 1; /* switch to a maxima */
568 else
569 {
570 alt = 0;
571 break; /* Ooops, found two non-alternating */
572 } /* extrema. Delete smallest of them */
573 } /* if the loop finishes, all extrema are alternating */
574
575 /*
576 * If there's only one extremal and all are alternating,
577 * delete the smallest of the first/last extremals.
578 */
579 if ((alt) && (extra == 1))
580 {
581 if (fabs(E[foundExt[k-1]]) < fabs(E[foundExt[0]]))
582 l = foundExt[k-1]; /* Delete last extremal */
583 else
584 l = foundExt[0]; /* Delete first extremal */
585 }
586
587 for (j=l; j<k; j++) /* Loop that does the deletion */
588 {
589 foundExt[j] = foundExt[j+1];
590 }
591 k--;
592 extra--;
593 }
594
595 for (i=0; i<=r; i++)
596 {
597 Ext[i] = foundExt[i]; /* Copy found extremals to Ext[] */
598 }
599
600 free(foundExt);
601 }
602
603
604 /*********************
605 * FreqSample
606 *============
607 * Simple frequency sampling algorithm to determine the impulse
608 * response h[] from A's found in ComputeA
609 *
610 *
611 * INPUT:
612 * ------
613 * int N - Number of filter coefficients
614 * double A[] - Sample points of desired response [N/2]
615 * int symmetry - Symmetry of desired filter
616 *
617 * OUTPUT:
618 * -------
619 * double h[] - Impulse Response of final filter [N]
620 *********************/
621 void FreqSample(int N, double A[], double h[], int symm)
622 {
623 int n, k;
624 double x, val, M;
625
626 M = (N-1.0)/2.0;
627 if (symm == POSITIVE)
628 {
629 if (N%2)
630 {
631 for (n=0; n<N; n++)
632 {
633 val = A[0];
634 x = Pi2 * (n - M)/N;
635 for (k=1; k<=M; k++)
636 val += 2.0 * A[k] * cos(x*k);
637 h[n] = val/N;
638 }
639 }
640 else
641 {
642 for (n=0; n<N; n++)
643 {
644 val = A[0];
645 x = Pi2 * (n - M)/N;
646 for (k=1; k<=(N/2-1); k++)
647 val += 2.0 * A[k] * cos(x*k);
648 h[n] = val/N;
649 }
650 }
651 }
652 else
653 {
654 if (N%2)
655 {
656 for (n=0; n<N; n++)
657 {
658 val = 0;
659 x = Pi2 * (n - M)/N;
660 for (k=1; k<=M; k++)
661 val += 2.0 * A[k] * sin(x*k);
662 h[n] = val/N;
663 }
664 }
665 else
666 {
667 for (n=0; n<N; n++)
668 {
669 val = A[N/2] * sin(Pi * (n - M));
670 x = Pi2 * (n - M)/N;
671 for (k=1; k<=(N/2-1); k++)
672 val += 2.0 * A[k] * sin(x*k);
673 h[n] = val/N;
674 }
675 }
676 }
677 }
678
679 /*******************
680 * isDone
681 *========
682 * Checks to see if the error function is small enough to consider
683 * the result to have converged.
684 *
685 * INPUT:
686 * ------
687 * int r - 1/2 the number of filter coeffiecients
688 * int Ext[] - Indexes to extremal frequencies [r+1]
689 * double E[] - Error function on the dense grid [gridsize]
690 *
691 * OUTPUT:
692 * -------
693 * Returns 1 if the result converged
694 * Returns 0 if the result has not converged
695 ********************/
696
697 short isDone(int r, int Ext[], double E[])
698 {
699 int i;
700 double min, max, current;
701
702 min = max = fabs(E[Ext[0]]);
703 for (i=1; i<=r; i++)
704 {
705 current = fabs(E[Ext[i]]);
706 if (current < min)
707 min = current;
708 if (current > max)
709 max = current;
710 }
711 if (((max-min)/max) < 0.0001)
712 return 1;
713 return 0;
714 }
715
716 /********************
717 * remez
718 *=======
719 * Calculates the optimal (in the Chebyshev/minimax sense)
720 * FIR filter impulse response given a set of band edges,
721 * the desired reponse on those bands, and the weight given to
722 * the error in those bands.
723 *
724 * INPUT:
725 * ------
726 * int numtaps - Number of filter coefficients
727 * int numband - Number of bands in filter specification
728 * double bands[] - User-specified band edges [2 * numband]
729 * double des[] - User-specified band responses [numband]
730 * double weight[] - User-specified error weights [numband]
731 * int type - Type of filter
732 *
733 * OUTPUT:
734 * -------
735 * double h[] - Impulse response of final filter [numtaps]
736 ********************/
737
738 void remez(double h[], int numtaps,
739 int numband, double bands[], double des[], double weight[],
740 int type)
741 {
742 double *Grid, *W, *D, *E;
743 int i, iter, gridsize, r, *Ext;
744 double *taps, c;
745 double *x, *y, *ad;
746 int symmetry;
747
748 if (type == BANDPASS)
749 symmetry = POSITIVE;
750 else
751 symmetry = NEGATIVE;
752
753 r = numtaps/2; /* number of extrema */
754 if ((numtaps%2) && (symmetry == POSITIVE))
755 r++;
756
757 /*
758 * Predict dense grid size in advance for memory allocation
759 * .5 is so we round up, not truncate
760 */
761 gridsize = 0;
762 for (i=0; i<numband; i++)
763 {
764 gridsize += (int)(2*r*GRIDDENSITY*(bands[2*i+1] - bands[2*i]) + .5);
765 }
766 if (symmetry == NEGATIVE)
767 {
768 gridsize--;
769 }
770
771 /*
772 * Dynamically allocate memory for arrays with proper sizes
773 */
774 Grid = (double *)malloc(gridsize * sizeof(double));
775 D = (double *)malloc(gridsize * sizeof(double));
776 W = (double *)malloc(gridsize * sizeof(double));
777 E = (double *)malloc(gridsize * sizeof(double));
778 Ext = (int *)malloc((r+1) * sizeof(int));
779 taps = (double *)malloc((r+1) * sizeof(double));
780 x = (double *)malloc((r+1) * sizeof(double));
781 y = (double *)malloc((r+1) * sizeof(double));
782 ad = (double *)malloc((r+1) * sizeof(double));
783
784 /*
785 * Create dense frequency grid
786 */
787 CreateDenseGrid(r, numtaps, numband, bands, des, weight,
788 &gridsize, Grid, D, W, symmetry);
789 InitialGuess(r, Ext, gridsize);
790
791 /*
792 * For Differentiator: (fix grid)
793 */
794 if (type == DIFFERENTIATOR)
795 {
796 for (i=0; i<gridsize; i++)
797 {
798 /* D[i] = D[i]*Grid[i]; */
799 if (D[i] > 0.0001)
800 W[i] = W[i]/Grid[i];
801 }
802 }
803
804 /*
805 * For odd or Negative symmetry filters, alter the
806 * D[] and W[] according to Parks McClellan
807 */
808 if (symmetry == POSITIVE)
809 {
810 if (numtaps % 2 == 0)
811 {
812 for (i=0; i<gridsize; i++)
813 {
814 c = cos(Pi * Grid[i]);
815 D[i] /= c;
816 W[i] *= c;
817 }
818 }
819 }
820 else
821 {
822 if (numtaps % 2)
823 {
824 for (i=0; i<gridsize; i++)
825 {
826 c = sin(Pi2 * Grid[i]);
827 D[i] /= c;
828 W[i] *= c;
829 }
830 }
831 else
832 {
833 for (i=0; i<gridsize; i++)
834 {
835 c = sin(Pi * Grid[i]);
836 D[i] /= c;
837 W[i] *= c;
838 }
839 }
840 }
841
842 /*
843 * Perform the Remez Exchange algorithm
844 */
845 for (iter=0; iter<MAXITERATIONS; iter++)
846 {
847 CalcParms(r, Ext, Grid, D, W, ad, x, y);
848 CalcError(r, ad, x, y, gridsize, Grid, D, W, E);
849 Search(r, Ext, gridsize, E);
850 if (isDone(r, Ext, E))
851 break;
852 }
853 if (iter == MAXITERATIONS)
854 {
855 printf("Reached maximum iteration count.\nResults may be bad.\n");
856 }
857
858 CalcParms(r, Ext, Grid, D, W, ad, x, y);
859
860 /*
861 * Find the 'taps' of the filter for use with Frequency
862 * Sampling. If odd or Negative symmetry, fix the taps
863 * according to Parks McClellan
864 */
865 for (i=0; i<=numtaps/2; i++)
866 {
867 if (symmetry == POSITIVE)
868 {
869 if (numtaps%2)
870 c = 1;
871 else
872 c = cos(Pi * (double)i/numtaps);
873 }
874 else
875 {
876 if (numtaps%2)
877 c = sin(Pi2 * (double)i/numtaps);
878 else
879 c = sin(Pi * (double)i/numtaps);
880 }
881 taps[i] = ComputeA((double)i/numtaps, r, ad, x, y)*c;
882 }
883
884 /*
885 * Frequency sampling design with calculated taps
886 */
887 FreqSample(numtaps, taps, h, symmetry);
888
889 /*
890 * Delete allocated memory
891 */
892 free(Grid);
893 free(W);
894 free(D);
895 free(E);
896 free(Ext);
897 free(x);
898 free(y);
899 free(ad);
900 }