MNE-CPP  beta 1.0
newparksmcclellan.cpp
1 /*
2 October 19, 2013
3 From: http://www.iowahills.com/A7ExampleCodePage.html
4 
5 The original fortran code came from the Parks McClellan article on Wikipedia.
6 http://en.wikipedia.org/wiki/Parks%E2%80%93McClellan_filter_design_algorithm
7 
8 This code is quite different from the original. The original code had 69 goto statements,
9 which made it nearly impossible to follow. And of course, it was Fortran code, so many changes
10 had to be made regardless of style.
11 
12 Our first step was to get a C version of the code working with as few changes as possible.
13 Then, in our desire to see if the code could be made more understandable, we decided to
14 remove as many goto statements as possible. We checked our work by comparing the coefficients
15 between this code and our original translation on more than 1000 filters while varying all the parameters.
16 
17 Ultimately, we were able to reduce the goto count from 69 to 7, all of which are in the Remez
18 function. Of the 7 remaining, 3 of these are at the very bottom of the function, and go
19 back to the very top of the function. These could have been removed, but our goal was to
20 clarify the code, not restyle it, and since they are clear, we let them be.
21 
22 The other 4 goto statements are intertwined in a rather nasty way. We recommend you print out
23 the Remez code, tape the sheets end to end, and trace out the goto's. It wasn't apparent to
24 us that they can be removed without an extensive study of the code.
25 
26 For better or worse, we also removed any code that was obviously related to Hilbert transforms
27 and Differentiators. We did this because we aren't interested in these, and we also don't
28 believe this algorithm does a very good job with them (far too much ripple).
29 
30 We added the functions CalcCoefficients() and ErrTest() as a way to simplify things a bit.
31 
32 We also found 3 sections of code that never executed. Two of the sections were just a few lines
33 that the goto's always went around. The third section represented nearly half of the CalcCoefficients()
34 function. This statement always tested the same, which never allowed the code to execute.
35 if(GRID[1] < 0.01 && GRID[NGRID] > 0.49) KKK = 1;
36 This may be due to the 0.01 minimum width limit we set for the bands.
37 
38 Note our use of MIN_TEST_VAL. The original code wasn't guarding against division by zero.
39 Limiting the return values as we have also helped the algorithm's convergence behavior.
40 
41 In an effort to improve readability, we made a large number of variable name changes and also
42 deleted a large number of variables. We left many variable names in tact, in part as an aid when
43 comparing to the original code, and in part because a better name wasn't obvious.
44 
45 This code is essentially straight c, and should compile with few, if any changes. Take note
46 of the 4 commented lines at the bottom of CalcParkCoeff2. These lines replace the code from the
47 original Ouch() function. They warn of the possibility of convergence failure, but you will
48 find that the iteration count NITER, isn't always an indicator of convergence problems when
49 it is less than 3, as stated in the original Fortran code comments.
50 
51 If you find a problem with this code, please leave us a note on:
52 http://www.iowahills.com/feedbackcomments.html
53 
54 */
55 
56 #ifdef __cplusplus
57 extern "C" {
58 #endif
59 
60 #include "newparksmcclellan.h"
61 #include <math.h>
62 #include <new> // defines new(std::nothrow)
63 
64 #include <cstddef>
65 
66 #define BIG 4096 // Used to define array sizes. Must be somewhat larger than 8 * MaxNumTaps
67 #define SMALL 256
68 #define M_2PI 6.28318530717958647692
69 #define ITRMAX 50 // Max Number of Iterations. Some filters require as many as 45 iterations.
70 #define MIN_TEST_VAL 1.0E-6 // Min value used in LeGrangeInterp and GEE
71 
72 int HalfTapCount, *ExchangeIndex;
73 double *LeGrangeD, *Alpha, *CosOfGrid, *DesPlus, *Coeff, *Edge, *BandMag, *InitWeight, *DesiredMag, *Grid, *Weight, *FirCoeff;
74 bool InitDone2 = false;
75 
76 // Using nothrow prevents an exception from being thrown. new will instead return NULL.
77 // These array are much larger than actually needed. See the notes in the orig fortran file.
78 void InitParksMcClellan2(void)
79 {
80  InitDone2 = true;
81 
82  ExchangeIndex = new(std::nothrow) int[SMALL];
83  if(ExchangeIndex == NULL)InitDone2 = false;
84 
85  LeGrangeD = new(std::nothrow) double[SMALL];
86  if(LeGrangeD == NULL)InitDone2 = false;
87 
88  Alpha = new(std::nothrow) double[SMALL];
89  if(Alpha == NULL)InitDone2 = false;
90 
91  CosOfGrid = new(std::nothrow) double[SMALL];
92  if(CosOfGrid == NULL)InitDone2 = false;
93 
94  DesPlus = new(std::nothrow) double[SMALL];
95  if(DesPlus == NULL)InitDone2 = false;
96 
97  Coeff = new(std::nothrow) double[SMALL];
98  if(Coeff == NULL)InitDone2 = false;
99 
100  Edge = new(std::nothrow) double[SMALL];
101  if(Edge == NULL)InitDone2 = false;
102 
103  BandMag = new(std::nothrow) double[SMALL];
104  if(BandMag == NULL)InitDone2 = false;
105 
106  InitWeight = new(std::nothrow) double[SMALL];
107  if(InitWeight == NULL)InitDone2 = false;
108 
109  DesiredMag = new(std::nothrow) double[BIG];
110  if(DesiredMag == NULL)InitDone2 = false;
111 
112  Grid = new(std::nothrow) double[BIG];
113  if(Grid == NULL)InitDone2 = false;
114 
115  Weight = new(std::nothrow) double[BIG];
116  if(Weight == NULL)InitDone2 = false;
117 
118  FirCoeff = new(std::nothrow) double[BIG];
119  if(FirCoeff == NULL)InitDone2 = false;
120 
121  //if(!InitDone2)MessageBox(NULL, L"Memory Error", L"Memory not Allocated in InitParksMcClellan", MB_OK);
122 
123 }
124 
125 //---------------------------------------------------------------------------
126 
127 
128 // NumTaps must be odd for high pass and notch filters. Max number of taps is 128.
129 // The arrays can handle up to 256 taps, but 128 is a good practical limit for convergence.
130 // The minimum number of taps is 9 (maybe < 9, I forget the exact lower limit or what sets it)
131 // OmegaC is the 3 dB corner freq for low pass and high pass filters.
132 // It is the center freq for band pass and notch filters.
133 // BW is the bandwidth for bandpass and notch filters (ignored on low and high pass).
134 // OmegaC and BW are in terms of Pi. e.g. OmegaC = 0.5 centers a BPF at Omega = Pi/2.
135 // The PM algorithm however uses frequencies in terms of 2Pi, so we need to to this: Edge[j] /= 2.0
136 // ParksWidth is the width of the transition bands. For simplicity, we only use one width,
137 // but the algorithm allows for unique values on every band edge.
138 // Practical limits for ParksWidth are 0.02 - 0.15 for BPF and Notch, 0.02 - 0.30 for LPF and HPF.
139 // TPassType is defined in the header file. LPF = Low Pass Filter, etc.
140 // You should note our 0.01 minimum width for each band. This limit works well for the algorithm.
141 // You will also find that OmegaC and BW need to be scaled a bit, depending ParksWidth, to get the
142 // 3 dB corner frequencies to come in on target.
143 
144 // e.g. NewParksMcClellan(33, 0.7, 0.2, 0.1, HPF);
145 // gives a 33 tap high pass filter with 3 dB corner at 0.7 with a transition bandwidth of 0.1
146 // The FIR coefficients are placed in FirCoeff, starting at index 0.
147 
148 void NewParksMcClellan(int NumTaps, double OmegaC, double BW, double ParksWidth, TPassType PassType)
149 {
150  int j, NumBands;
151 
152  if(!InitDone2)
153  {
154  InitParksMcClellan2();
155  if(!InitDone2)return;
156  }
157 
158  if(NumTaps > 128)NumTaps = 128;
159  if(NumTaps < 9)NumTaps = 9;
160  if( (PassType == HPF || PassType == NOTCH) && NumTaps % 2 == 0)NumTaps--;
161 
162  // It helps the algorithm a great deal if each band is at least 0.01 wide.
163  // The weights used here came from the orig PM code.
164  if(PassType == LPF)
165  {
166  NumBands = 2;
167  Edge[1] = 0.0; // Omega = 0
168  Edge[2] = OmegaC; // Pass band edge
169  if(Edge[2] < 0.01)Edge[2] = 0.01;
170  if(Edge[2] > 0.98)Edge[2] = 0.98;
171  Edge[3] = Edge[2] + ParksWidth; // Stop band edge
172  if(Edge[3] > 0.99)Edge[3] = 0.99;
173  Edge[4] = 1.0; // Omega = Pi
174 
175  BandMag[1] = 1.0;
176  BandMag[2] = 0.0;
177  InitWeight[1] = 1.0;
178  InitWeight[2] = 10.0;
179  }
180 
181  if(PassType == HPF)
182  {
183  NumBands = 2;
184  Edge[1] = 0.0; // Omega = 0
185  Edge[3] = OmegaC; // Pass band edge
186  if(Edge[3] > 0.99)Edge[3] = 0.99;
187  if(Edge[3] < 0.02)Edge[3] = 0.02;
188  Edge[2] = Edge[3] - ParksWidth; // Stop band edge
189  if(Edge[2] < 0.01)Edge[2] = 0.01;
190  Edge[4] = 1.0; // Omega = Pi
191 
192  BandMag[1] = 0.0;
193  BandMag[2] = 1.0;
194  InitWeight[1] = 10.0;
195  InitWeight[2] = 1.0;
196  }
197 
198  if(PassType == BPF)
199  {
200  NumBands = 3;
201  Edge[1] = 0.0; // Omega = 0
202  Edge[3] = OmegaC - BW/2.0; // Left pass band edge.
203  if(Edge[3] < 0.02)Edge[3] = 0.02;
204  Edge[2] = Edge[3] - ParksWidth; // Left stop band edge
205  if(Edge[2] < 0.01)Edge[2] = 0.01;
206  Edge[4] = OmegaC + BW/2.0; // Right pass band edge
207  if(Edge[4] > 0.98)Edge[4] = 0.98;
208  Edge[5] = Edge[4] + ParksWidth; // Right stop band edge
209  if(Edge[5] > 0.99)Edge[5] = 0.99;
210  Edge[6] = 1.0; // Omega = Pi
211 
212  BandMag[1] = 0.0;
213  BandMag[2] = 1.0;
214  BandMag[3] = 0.0;
215  InitWeight[1] = 10.0;
216  InitWeight[2] = 1.0;
217  InitWeight[3] = 10.0;
218  }
219 
220  if(PassType == NOTCH)
221  {
222  NumBands = 3;
223  Edge[1] = 0.0; // Omega = 0
224  Edge[3] = OmegaC - BW/2.0; // Left stop band edge.
225  if(Edge[3] < 0.02)Edge[3] = 0.02;
226  Edge[2] = Edge[3] - ParksWidth; // Left pass band edge
227  if(Edge[2] < 0.01)Edge[2] = 0.01;
228  Edge[4] = OmegaC + BW/2.0; // Right stop band edge
229  if(Edge[4] > 0.98)Edge[4] = 0.98;
230  Edge[5] = Edge[4] + ParksWidth; // Right pass band edge
231  if(Edge[5] > 0.99)Edge[5] = 0.99;
232  Edge[6] = 1.0; // Omega = Pi
233 
234  BandMag[1] = 1.0;
235  BandMag[2] = 0.0;
236  BandMag[3] = 1.0;
237  InitWeight[1] = 1.0;
238  InitWeight[2] = 10.0;
239  InitWeight[3] = 1.0;
240  }
241 
242  // Parks McClellan's edges are based on 2Pi, we are based on Pi.
243  for(j=1; j<=2*NumBands; j++)Edge[j] /= 2.0;
244 
245  CalcParkCoeff2(NumBands, NumTaps);
246 
247 }
248 //---------------------------------------------------------------------------
249 
250 
251 void CalcParkCoeff2(int NumBands, int TapCount)
252 {
253  int j, k, GridCount, GridIndex, BandIndex, NumIterations;
254  double LowFreqEdge, UpperFreq, TempVar, Change;
255  bool OddNumTaps;
256  GridCount = 16; // Grid Density
257 
258  if(TapCount % 2)OddNumTaps = true;
259  else OddNumTaps = false;
260 
261  HalfTapCount = TapCount/2;
262  if(OddNumTaps) HalfTapCount++;
263 
264  Grid[1] = Edge[1];
265  LowFreqEdge = GridCount * HalfTapCount;
266  LowFreqEdge = 0.5 / LowFreqEdge;
267  j = 1;
268  k = 1;
269  BandIndex = 1;
270  while(BandIndex <= NumBands)
271  {
272  UpperFreq = Edge[k+1];
273  while(Grid[j] <= UpperFreq)
274  {
275  TempVar = Grid[j];
276  DesiredMag[j] = BandMag[BandIndex];
277  Weight[j] = InitWeight[BandIndex];
278  j++;;
279  Grid[j] = TempVar + LowFreqEdge;
280  }
281 
282  Grid[j-1] = UpperFreq;
283  DesiredMag[j-1] = BandMag[BandIndex];
284  Weight[j-1] = InitWeight[BandIndex];
285  k+=2;
286  BandIndex++;
287  if(BandIndex <= NumBands)Grid[j] = Edge[k];
288  }
289 
290  GridIndex = j-1;
291  if(!OddNumTaps && Grid[GridIndex] > (0.5-LowFreqEdge)) GridIndex--;
292 
293  if(!OddNumTaps)
294  {
295  for(j=1; j<=GridIndex; j++)
296  {
297  Change = cos(M_PI * Grid[j] );
298  DesiredMag[j] = DesiredMag[j] / Change;
299  Weight[j] = Weight[j] * Change;
300  }
301  }
302 
303  TempVar = (double)(GridIndex-1)/(double)HalfTapCount;
304  for(j=1; j<=HalfTapCount; j++)
305  {
306  ExchangeIndex[j] = (double)(j-1) * TempVar + 1.0;
307  }
308  ExchangeIndex[HalfTapCount+1] = GridIndex;
309 
310  NumIterations = Remez2(GridIndex);
311  CalcCoefficients();
312 
313  // Calculate the impulse response.
314  if(OddNumTaps)
315  {
316  for(j=1; j<=HalfTapCount-1; j++)
317  {
318  Coeff[j] = 0.5 * Alpha[HalfTapCount+1-j];
319  }
320  Coeff[HalfTapCount] = Alpha[1];
321  }
322  else
323  {
324  Coeff[1] = 0.25 * Alpha[HalfTapCount];
325  for(j=2; j<=HalfTapCount-1; j++)
326  {
327  Coeff[j] = 0.25 * (Alpha[HalfTapCount+1-j] + Alpha[HalfTapCount+2-j]);
328  }
329  Coeff[HalfTapCount] = 0.5 * Alpha[1] + 0.25 * Alpha[2];
330  }
331 
332 
333  // Output section.
334  for(j=1; j<=HalfTapCount; j++)FirCoeff[j-1] = Coeff[j];
335  if(OddNumTaps)
336  for(j=1; j<HalfTapCount; j++)FirCoeff[HalfTapCount+j-1] = Coeff[HalfTapCount-j];
337  else
338  for(j=1; j<=HalfTapCount; j++)FirCoeff[HalfTapCount+j-1] = Coeff[HalfTapCount-j+1];
339 
340  // Parks2Label was on my application's main form.
341  // These replace the original Ouch() function
342  if(NumIterations <= 3)
343  {
344  //TopForm->Parks2Label->Font->Color = clRed;
345  //TopForm->Parks2Label->Caption = "Convergence ?"; // Covergence is doubtful, but possible.
346  }
347  else
348  {
349  //TopForm->Parks2Label->Font->Color = clBlack;
350  //TopForm->Parks2Label->Caption = UnicodeString(NumIterations) + " Iterations";
351  }
352 
353 }
354 
355 //---------------------------------------------------------------------------------------
356 int Remez2(int GridIndex)
357 {
358  int j, JET, K, k, NU, JCHNGE, K1, KNZ, KLOW, NUT, KUP;
359  int NUT1, LUCK, KN, NITER;
360  double Deviation, DNUM, DDEN, TempVar;
361  double DEVL, COMP, YNZ, Y1, ERR;
362 
363  LUCK = 0;
364  DEVL = -1.0;
365  NITER = 1; // Init this to 1 to be consistent with the orig code.
366 
367  TOP_LINE: // We come back to here from 3 places at the bottom.
368  ExchangeIndex[HalfTapCount+2] = GridIndex + 1;
369 
370  for(j=1; j<=HalfTapCount+1; j++)
371  {
372  TempVar = Grid[ ExchangeIndex[j] ];
373  CosOfGrid[j] = cos(TempVar * M_2PI);
374  }
375 
376  JET = (HalfTapCount-1)/15 + 1;
377  for(j=1; j<=HalfTapCount+1; j++)
378  {
379  LeGrangeD[j] = LeGrangeInterp2(j,HalfTapCount+1,JET);
380  }
381 
382  DNUM = 0.0;
383  DDEN = 0.0;
384  K = 1;
385  for(j=1; j<=HalfTapCount+1; j++)
386  {
387  k = ExchangeIndex[j];
388  DNUM += LeGrangeD[j] * DesiredMag[k];
389  DDEN += (double)K * LeGrangeD[j]/Weight[k];
390  K = -K;
391  }
392  Deviation = DNUM / DDEN;
393 
394  NU = 1;
395  if(Deviation > 0.0) NU = -1;
396  Deviation = -(double)NU * Deviation;
397  K = NU;
398  for(j=1; j<=HalfTapCount+1; j++)
399  {
400  k = ExchangeIndex[j];
401  TempVar = (double)K * Deviation/Weight[k];
402  DesPlus[j] = DesiredMag[k] + TempVar;
403  K = -K;
404  }
405 
406  if(Deviation <= DEVL)return(NITER); // Ouch
407 
408  DEVL = Deviation;
409  JCHNGE = 0;
410  K1 = ExchangeIndex[1];
411  KNZ = ExchangeIndex[HalfTapCount+1];
412  KLOW = 0;
413  NUT = -NU;
414  j = 1;
415 
416  //Search for the extremal frequencies of the best approximation.
417 
418  j=1;
419  while(j<HalfTapCount+2)
420  {
421  KUP = ExchangeIndex[j+1];
422  k = ExchangeIndex[j] + 1;
423  NUT = -NUT;
424  if(j == 2) Y1 = COMP;
425  COMP = Deviation;
426 
427  if(k < KUP && !ErrTest(k, NUT, COMP, &ERR))
428  {
429  L210:
430  COMP = (double)NUT * ERR;
431  for(k++; k<KUP; k++)
432  {
433  if( ErrTest(k, NUT, COMP, &ERR) )break; // for loop
434  COMP = (double)NUT * ERR;
435  }
436 
437  ExchangeIndex[j] = k-1;
438  j++;
439  KLOW = k - 1;
440  JCHNGE++;
441  continue; // while loop
442  }
443 
444  k--;
445 
446  L225: k--;
447  if(k <= KLOW)
448  {
449  k = ExchangeIndex[j] + 1;
450  if(JCHNGE > 0)
451  {
452  ExchangeIndex[j] = k-1;
453  j++;
454  KLOW = k - 1;
455  JCHNGE++;
456  continue; // while loop
457  }
458  else // JCHNGE <= 0
459  {
460  for(k++; k<KUP; k++)
461  {
462  if( ErrTest(k, NUT, COMP, &ERR) )continue; // for loop
463  goto L210;
464  }
465 
466  KLOW = ExchangeIndex[j];
467  j++;
468  continue; // while loop
469  }
470  }
471  // Can't use a do while loop here, it would outdent the two continue statements.
472  if( ErrTest(k, NUT, COMP, &ERR) && JCHNGE <= 0)goto L225;
473 
474  if( ErrTest(k, NUT, COMP, &ERR) )
475  {
476  KLOW = ExchangeIndex[j];
477  j++;
478  continue; // while loop
479  }
480 
481  COMP = (double)NUT * ERR;
482 
483  L235:
484  for(k--; k>KLOW; k--)
485  {
486  if( ErrTest(k, NUT, COMP, &ERR) ) break; // for loop
487  COMP = (double)NUT * ERR;
488  }
489 
490  KLOW = ExchangeIndex[j];
491  ExchangeIndex[j] = k + 1;
492  j++;
493  JCHNGE++;
494  } // end while(j<HalfTapCount
495 
496  if(j == HalfTapCount+2) YNZ = COMP;
497 
498  while(j <= HalfTapCount+2)
499  {
500  if(K1 > ExchangeIndex[1]) K1 = ExchangeIndex[1];
501  if(KNZ < ExchangeIndex[HalfTapCount+1]) KNZ = ExchangeIndex[HalfTapCount+1];
502  NUT1 = NUT;
503  NUT = -NU;
504  k = 0 ;
505  KUP = K1;
506  COMP = YNZ * 1.00001;
507  LUCK = 1;
508 
509  for(k++; k<KUP; k++)
510  {
511  if( ErrTest(k, NUT, COMP, &ERR) ) continue; // for loop
512  j = HalfTapCount+2;
513  goto L210;
514  }
515  LUCK = 2;
516  break; // break while(j <= HalfTapCount+2) loop
517  } // end while(j <= HalfTapCount+2)
518 
519  if(LUCK == 1 || LUCK == 2)
520  {
521  if(LUCK == 1)
522  {
523  if(COMP > Y1) Y1 = COMP;
524  K1 = ExchangeIndex[HalfTapCount+2];
525  }
526 
527  k = GridIndex + 1;
528  KLOW = KNZ;
529  NUT = -NUT1;
530  COMP = Y1 * 1.00001;
531 
532  for(k--; k>KLOW; k--)
533  {
534  if( ErrTest(k, NUT, COMP, &ERR) )continue; // for loop
535  j = HalfTapCount+2;
536  COMP = (double)NUT * ERR;
537  LUCK = 3; // last time in this if(LUCK == 1 || LUCK == 2)
538  goto L235;
539  }
540 
541  if(LUCK == 2)
542  {
543  if(JCHNGE > 0 && NITER++ < ITRMAX) goto TOP_LINE;
544  else return(NITER);
545  }
546 
547  for(j=1; j<=HalfTapCount; j++)
548  {
549  ExchangeIndex[HalfTapCount+2 - j] = ExchangeIndex[HalfTapCount+1 - j];
550  }
551  ExchangeIndex[1] = K1;
552  if(NITER++ < ITRMAX) goto TOP_LINE;
553  } // end if(LUCK == 1 || LUCK == 2)
554 
555 
556 
557  KN = ExchangeIndex[HalfTapCount+2];
558  for(j=1; j<=HalfTapCount; j++)
559  {
560  ExchangeIndex[j] = ExchangeIndex[j+1];
561  }
562  ExchangeIndex[HalfTapCount+1] = KN;
563  if(NITER++ < ITRMAX) goto TOP_LINE;
564 
565 
566  return(NITER);
567 
568 }
569 
570 //-----------------------------------------------------------------------
571 // Function to calculate the lagrange interpolation coefficients for use in the function gee.
572 double LeGrangeInterp2(int K, int N, int M) // D
573 {
574  int j, k;
575  double Dee, Q;
576  Dee = 1.0;
577  Q = CosOfGrid[K];
578  for(k=1; k<=M; k++)
579  for(j=k; j<=N; j+=M)
580  {
581  if(j != K)Dee = 2.0 * Dee * (Q - CosOfGrid[j]);
582  }
583  if(fabs(Dee) < MIN_TEST_VAL )
584  {
585  if(Dee < 0.0)Dee = -MIN_TEST_VAL;
586  else Dee = MIN_TEST_VAL;
587  }
588  return(1.0/Dee);
589 }
590 
591 //-----------------------------------------------------------------------
592 // Function to evaluate the frequency response using the Lagrange interpolation
593 // formula in the barycentric form.
594 double GEE2(int K, int N)
595 {
596  int j;
597  double P,C,Dee,XF;
598  P = 0.0;
599  XF = Grid[K];
600  XF = cos(M_2PI * XF);
601  Dee = 0.0;
602  for(j=1; j<=N; j++)
603  {
604  C = XF - CosOfGrid[j];
605  if(fabs(C) < MIN_TEST_VAL )
606  {
607  if(C < 0.0)C = -MIN_TEST_VAL;
608  else C = MIN_TEST_VAL;
609  }
610  C = LeGrangeD[j] / C;
611  Dee = Dee + C;
612  P = P + C*DesPlus[j];
613  }
614  if(fabs(Dee) < MIN_TEST_VAL )
615  {
616  if(Dee < 0.0)Dee = -MIN_TEST_VAL;
617  else Dee = MIN_TEST_VAL;
618  }
619  return(P/Dee);
620 }
621 
622 //-----------------------------------------------------------------------
623 
624 // This was added by IowaHills and is used in Remez() in 6 places.
625 bool ErrTest(int k, int Nut, double Comp, double *Err)
626 {
627  *Err = GEE2(k, HalfTapCount+1);
628  *Err = (*Err - DesiredMag[k]) * Weight[k];
629  if((double)Nut * *Err - Comp <= 0.0) return(true);
630  else return(false);
631 }
632 
633 //-----------------------------------------------------------------------
634 
635 // This was added by IowaHills and is called from CalcParkCoeff2().
636 // Calculation of the coefficients of the best approximation using the inverse discrete fourier transform.
637 void CalcCoefficients(void)
638 {
639  int j, k, n;
640  double GTempVar, OneOverNumTaps;
641  double Omega, TempVar, FreqN, TempX, GridCos;
642  double GeeArray[SMALL];
643 
644  GTempVar = Grid[1];
645  CosOfGrid[HalfTapCount+2] = -2.0;
646  OneOverNumTaps = 1.0 /(double)(2*HalfTapCount-1);
647  k = 1;
648 
649  for(j=1; j<=HalfTapCount; j++)
650  {
651  FreqN = (double)(j-1) * OneOverNumTaps;
652  TempX = cos(M_2PI * FreqN);
653 
654  GridCos = CosOfGrid[k];
655  if(TempX <= GridCos)
656  {
657  while(TempX <= GridCos && (GridCos-TempX) >= MIN_TEST_VAL) // MIN_TEST_VAL = 1.0E-6
658  {
659  k++;;
660  GridCos = CosOfGrid[k];
661  }
662  }
663  if(TempX <= GridCos || (TempX-GridCos) < MIN_TEST_VAL)
664  {
665  GeeArray[j] = DesPlus[k]; // Desired Response
666  }
667  else
668  {
669  Grid[1] = FreqN;
670  GeeArray[j] = GEE2(1, HalfTapCount+1);
671  }
672  if(k > 1) k--;
673  }
674 
675  Grid[1] = GTempVar;
676  for(j=1; j<=HalfTapCount; j++)
677  {
678  TempVar = 0.0;
679  Omega = (double)(j-1) * M_2PI * OneOverNumTaps;
680  for(n=1; n<=HalfTapCount-1; n++)
681  {
682  TempVar += GeeArray[n+1] * cos(Omega * (double)n);
683  }
684  TempVar = 2.0 * TempVar + GeeArray[1];
685  Alpha[j] = TempVar;
686  }
687 
688  Alpha[1] = Alpha[1] * OneOverNumTaps;
689  for(j=2; j<=HalfTapCount; j++)
690  {
691  Alpha[j] = 2.0 * Alpha[j] * OneOverNumTaps;
692  }
693 
694 }
695 
696 
697 #ifdef __cplusplus
698 }
699 #endif