68 #include <QtCore/qmath.h>
84 ParksMcClellan::ParksMcClellan()
91 ParksMcClellan::ParksMcClellan(
int NumTaps,
double OmegaC,
double BW,
double ParksWidth, TPassType PassType)
92 : ExchangeIndex(SMALL)
105 FirCoeff = RowVectorXd::Zero(NumTaps);
106 init(NumTaps, OmegaC, BW, ParksWidth, PassType);
112 ParksMcClellan::~ParksMcClellan()
118 void ParksMcClellan::init(
int NumTaps,
double OmegaC,
double BW,
double ParksWidth, TPassType PassType)
122 if(NumTaps > 256) NumTaps = 256;
123 if(NumTaps < 9) NumTaps = 9;
124 if( (PassType == HPF || PassType == NOTCH) && NumTaps % 2 == 0) NumTaps--;
133 if(Edge[2] < 0.01)Edge[2] = 0.01;
134 if(Edge[2] > 0.98)Edge[2] = 0.98;
135 Edge[3] = Edge[2] + ParksWidth;
136 if(Edge[3] > 0.99)Edge[3] = 0.99;
142 InitWeight[2] = 10.0;
150 if(Edge[3] > 0.99)Edge[3] = 0.99;
151 if(Edge[3] < 0.02)Edge[3] = 0.02;
152 Edge[2] = Edge[3] - ParksWidth;
153 if(Edge[2] < 0.01)Edge[2] = 0.01;
158 InitWeight[1] = 10.0;
166 Edge[3] = OmegaC - BW/2.0;
167 if(Edge[3] < 0.02)Edge[3] = 0.02;
168 Edge[2] = Edge[3] - ParksWidth;
169 if(Edge[2] < 0.01)Edge[2] = 0.01;
170 Edge[4] = OmegaC + BW/2.0;
171 if(Edge[4] > 0.98)Edge[4] = 0.98;
172 Edge[5] = Edge[4] + ParksWidth;
173 if(Edge[5] > 0.99)Edge[5] = 0.99;
179 InitWeight[1] = 10.0;
181 InitWeight[3] = 10.0;
184 if(PassType == NOTCH)
188 Edge[3] = OmegaC - BW/2.0;
189 if(Edge[3] < 0.02)Edge[3] = 0.02;
190 Edge[2] = Edge[3] - ParksWidth;
191 if(Edge[2] < 0.01)Edge[2] = 0.01;
192 Edge[4] = OmegaC + BW/2.0;
193 if(Edge[4] > 0.98)Edge[4] = 0.98;
194 Edge[5] = Edge[4] + ParksWidth;
195 if(Edge[5] > 0.99)Edge[5] = 0.99;
202 InitWeight[2] = 10.0;
207 for(j=1; j<=2*NumBands; j++) Edge[j] /= 2.0;
209 CalcParkCoeff2(NumBands, NumTaps);
215 void ParksMcClellan::CalcParkCoeff2(
int NumBands,
int TapCount)
217 int j, k, GridCount, GridIndex, BandIndex, NumIterations;
218 double LowFreqEdge, UpperFreq, TempVar, Change;
222 if(TapCount % 2)OddNumTaps =
true;
223 else OddNumTaps =
false;
225 HalfTapCount = TapCount/2;
226 if(OddNumTaps) HalfTapCount++;
229 LowFreqEdge = GridCount * HalfTapCount;
230 LowFreqEdge = 0.5 / LowFreqEdge;
234 while(BandIndex <= NumBands)
236 UpperFreq = Edge[k+1];
237 while(Grid[j] <= UpperFreq)
240 DesiredMag[j] = BandMag[BandIndex];
241 Weight[j] = InitWeight[BandIndex];
243 Grid[j] = TempVar + LowFreqEdge;
246 Grid[j-1] = UpperFreq;
247 DesiredMag[j-1] = BandMag[BandIndex];
248 Weight[j-1] = InitWeight[BandIndex];
251 if(BandIndex <= NumBands)Grid[j] = Edge[k];
255 if(!OddNumTaps && Grid[GridIndex] > (0.5-LowFreqEdge)) GridIndex--;
259 for(j=1; j<=GridIndex; j++)
261 Change = cos(M_PI * Grid[j] );
262 DesiredMag[j] = DesiredMag[j] / Change;
263 Weight[j] = Weight[j] * Change;
267 TempVar = (double)(GridIndex-1)/(double)HalfTapCount;
268 for(j=1; j<=HalfTapCount; j++)
270 ExchangeIndex[j] = (double)(j-1) * TempVar + 1.0;
272 ExchangeIndex[HalfTapCount+1] = GridIndex;
274 NumIterations = Remez2(GridIndex);
280 for(j=1; j<=HalfTapCount-1; j++)
282 Coeff[j] = 0.5 * Alpha[HalfTapCount+1-j];
284 Coeff[HalfTapCount] = Alpha[1];
288 Coeff[1] = 0.25 * Alpha[HalfTapCount];
289 for(j=2; j<=HalfTapCount-1; j++)
291 Coeff[j] = 0.25 * (Alpha[HalfTapCount+1-j] + Alpha[HalfTapCount+2-j]);
293 Coeff[HalfTapCount] = 0.5 * Alpha[1] + 0.25 * Alpha[2];
298 for(j=1; j<=HalfTapCount; j++) FirCoeff[j-1] = Coeff[j];
300 for(j=1; j<HalfTapCount; j++) FirCoeff[HalfTapCount+j-1] = Coeff[HalfTapCount-j];
302 for(j=1; j<=HalfTapCount; j++ )FirCoeff[HalfTapCount+j-1] = Coeff[HalfTapCount-j+1];
304 FirCoeff.conservativeResize(TapCount);
309 if(NumIterations <= 3)
324 int ParksMcClellan::Remez2(
int GridIndex)
326 int j, JET, K, k, NU, JCHNGE, K1, KNZ, KLOW, NUT, KUP;
327 int NUT1, LUCK, KN, NITER;
328 double Deviation, DNUM, DDEN, TempVar;
329 double DEVL, COMP, YNZ, Y1, ERR;
336 ExchangeIndex[HalfTapCount+2] = GridIndex + 1;
338 for(j=1; j<=HalfTapCount+1; j++)
340 TempVar = Grid[ ExchangeIndex[j] ];
341 CosOfGrid[j] = cos(TempVar * M_2PI);
344 JET = (HalfTapCount-1)/15 + 1;
345 for(j=1; j<=HalfTapCount+1; j++)
347 LeGrangeD[j] = LeGrangeInterp2(j,HalfTapCount+1,JET);
353 for(j=1; j<=HalfTapCount+1; j++)
355 k = ExchangeIndex[j];
356 DNUM += LeGrangeD[j] * DesiredMag[k];
357 DDEN += (double)K * LeGrangeD[j]/Weight[k];
360 Deviation = DNUM / DDEN;
363 if(Deviation > 0.0) NU = -1;
364 Deviation = -(double)NU * Deviation;
366 for(j=1; j<=HalfTapCount+1; j++)
368 k = ExchangeIndex[j];
369 TempVar = (double)K * Deviation/Weight[k];
370 DesPlus[j] = DesiredMag[k] + TempVar;
374 if(Deviation <= DEVL)
return(NITER);
378 K1 = ExchangeIndex[1];
379 KNZ = ExchangeIndex[HalfTapCount+1];
387 while(j<HalfTapCount+2)
389 KUP = ExchangeIndex[j+1];
390 k = ExchangeIndex[j] + 1;
392 if(j == 2) Y1 = COMP;
395 if(k < KUP && !ErrTest(k, NUT, COMP, &ERR))
398 COMP = (double)NUT * ERR;
401 if( ErrTest(k, NUT, COMP, &ERR) )
break;
402 COMP = (double)NUT * ERR;
405 ExchangeIndex[j] = k-1;
417 k = ExchangeIndex[j] + 1;
420 ExchangeIndex[j] = k-1;
430 if( ErrTest(k, NUT, COMP, &ERR) )
continue;
434 KLOW = ExchangeIndex[j];
440 if( ErrTest(k, NUT, COMP, &ERR) && JCHNGE <= 0)
goto L225;
442 if( ErrTest(k, NUT, COMP, &ERR) )
444 KLOW = ExchangeIndex[j];
449 COMP = (double)NUT * ERR;
452 for(k--; k>KLOW; k--)
454 if( ErrTest(k, NUT, COMP, &ERR) )
break;
455 COMP = (double)NUT * ERR;
458 KLOW = ExchangeIndex[j];
459 ExchangeIndex[j] = k + 1;
464 if(j == HalfTapCount+2) YNZ = COMP;
466 while(j <= HalfTapCount+2)
468 if(K1 > ExchangeIndex[1]) K1 = ExchangeIndex[1];
469 if(KNZ < ExchangeIndex[HalfTapCount+1]) KNZ = ExchangeIndex[HalfTapCount+1];
474 COMP = YNZ * 1.00001;
479 if( ErrTest(k, NUT, COMP, &ERR) )
continue;
487 if(LUCK == 1 || LUCK == 2)
491 if(COMP > Y1) Y1 = COMP;
492 K1 = ExchangeIndex[HalfTapCount+2];
500 for(k--; k>KLOW; k--)
502 if( ErrTest(k, NUT, COMP, &ERR) )
continue;
504 COMP = (double)NUT * ERR;
511 if(JCHNGE > 0 && NITER++ < ITRMAX)
goto TOP_LINE;
515 for(j=1; j<=HalfTapCount; j++)
517 ExchangeIndex[HalfTapCount+2 - j] = ExchangeIndex[HalfTapCount+1 - j];
519 ExchangeIndex[1] = K1;
520 if(NITER++ < ITRMAX)
goto TOP_LINE;
525 KN = ExchangeIndex[HalfTapCount+2];
526 for(j=1; j<=HalfTapCount; j++)
528 ExchangeIndex[j] = ExchangeIndex[j+1];
530 ExchangeIndex[HalfTapCount+1] = KN;
531 if(NITER++ < ITRMAX)
goto TOP_LINE;
539 double ParksMcClellan::LeGrangeInterp2(
int K,
int N,
int M)
548 if(j != K)Dee = 2.0 * Dee * (Q - CosOfGrid[j]);
550 if(fabs(Dee) < MIN_TEST_VAL )
552 if(Dee < 0.0)Dee = -MIN_TEST_VAL;
553 else Dee = MIN_TEST_VAL;
561 double ParksMcClellan::GEE2(
int K,
int N)
567 XF = cos(M_2PI * XF);
571 C = XF - CosOfGrid[j];
572 if(fabs(C) < MIN_TEST_VAL )
574 if(C < 0.0)C = -MIN_TEST_VAL;
575 else C = MIN_TEST_VAL;
577 C = LeGrangeD[j] / C;
579 P = P + C*DesPlus[j];
581 if(fabs(Dee) < MIN_TEST_VAL )
583 if(Dee < 0.0)Dee = -MIN_TEST_VAL;
584 else Dee = MIN_TEST_VAL;
592 bool ParksMcClellan::ErrTest(
int k,
int Nut,
double Comp,
double *Err)
594 *Err = GEE2(k, HalfTapCount+1);
595 *Err = (*Err - DesiredMag[k]) * Weight[k];
596 if((
double)Nut * *Err - Comp <= 0.0)
return(
true);
603 void ParksMcClellan::CalcCoefficients()
606 double GTempVar, OneOverNumTaps;
607 double Omega, TempVar, FreqN, TempX, GridCos;
608 double GeeArray[SMALL];
611 CosOfGrid[HalfTapCount+2] = -2.0;
612 OneOverNumTaps = 1.0 /(double)(2*HalfTapCount-1);
615 for(j=1; j<=HalfTapCount; j++)
617 FreqN = (double)(j-1) * OneOverNumTaps;
618 TempX = cos(M_2PI * FreqN);
620 GridCos = CosOfGrid[k];
623 while(TempX <= GridCos && (GridCos-TempX) >= MIN_TEST_VAL)
626 GridCos = CosOfGrid[k];
629 if(TempX <= GridCos || (TempX-GridCos) < MIN_TEST_VAL)
631 GeeArray[j] = DesPlus[k];
636 GeeArray[j] = GEE2(1, HalfTapCount+1);
642 for(j=1; j<=HalfTapCount; j++)
645 Omega = (double)(j-1) * M_2PI * OneOverNumTaps;
646 for(n=1; n<=HalfTapCount-1; n++)
648 TempVar += GeeArray[n+1] * cos(Omega * (
double)n);
650 TempVar = 2.0 * TempVar + GeeArray[1];
654 Alpha[1] = Alpha[1] * OneOverNumTaps;
655 for(j=2; j<=HalfTapCount; j++)
657 Alpha[j] = 2.0 * Alpha[j] * OneOverNumTaps;