/* @(#) $Id: QuantizationCorrection.cpp,v 1.3 2008/12/19 14:56:43 ramestic Exp $ * * Copyright (C) 2005 * Associated Universities, Inc. Washington DC, USA. * * This library is free software; you can redistribute it and/or modify it * under the terms of the GNU Library General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at your * option) any later version. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * Correspondence concerning ALMA Software should be addressed as follows: * Internet email: alma-sw-admin@nrao.edu */ // // System stuff // #include // // ACS stuff // // // CORR stuff // // // Local stuff // #include "QuantizationCorrection.h" using namespace std; // // it seems that the most efficient way for evaluating a polynomial like this // // y = a x + b x^3 + c x^5 // // it is to factorize it like this: // // y = a x [ 1 + (b x^2 / a) (1 + (c x^2 / b))] // // This suggests then that instead of keeping in memory the values of // a, b and c it is better to store a, b/a and c/b. That way, the above // expression can be repeated without the nead for recalculation every time // the original coefficients fractions. // // The following macro implements that computation. // #define EVAL_POLY(y, x, a, b_a, c_b) \ do { \ double _x2; \ _x2 = (x) * (x); \ y = a * (x) * (1.0 + b_a * _x2 * (1.0 + c_b * _x2)); \ } while(0) // // used in many places (x to power of 2) // #define STAR2(x) ((x) * (x)) #define STAR7(x) ((x) * (x) * (x) * (x) * (x) * (x) * (x)) // // initialize non-integral static memeber // const double QuantizationCorrection::PolynomialAproximationLimit_m(0.22); //------------------------------------------------------------------------------ QuantizationCorrection::QuantizationCorrection() { l1_m = sdp_real(0.0); l2_m = sdp_real(0.0); norm_m = sdp_real(0.0); rho_m = vector(0); Rn_m = vector(0); spl_m = vector(0); a_m = b_m = c_m = 0.0; ap_m = bp_m = cp_m = 0.0; b_a_m = c_b_m = 0.0; bp_ap_m = cp_bp_m = 0.0; Vs_m = 1.0; } //------------------------------------------------------------------------------ QuantizationCorrection::QuantizationCorrection(unsigned int N, unsigned int n, sdp_real sigma1, sdp_real sigma2, sdp_real Vs) { // // keep record of the number of quantization levels the instance is being // built for. // N_m = N; // // keep record of Vs and signal levels // Vs_m = Vs; // // record signal level to which the output corrected lags are to // be normalized. // norm_m = sigma1 * sigma2; // // compute coefficients for polynomial approximation (see eqn. 8): // // rho = a_m * Rn + b_m * Rn^3 + c_m * Rn^5 // double _2_pi = 2.0 / M_PI; l1_m = 1.0 / sigma1; l2_m = 1.0 / sigma2; a_m = _2_pi * Ck(N, 1, l1_m) * Ck(N, 1, l2_m); b_m = _2_pi * Ck(N, 3, l1_m) * Ck(N, 3, l2_m) / 6.0; c_m = _2_pi * Ck(N, 5, l1_m) * Ck(N, 5, l2_m) / 120.0; ap_m = 1.0 / a_m; bp_m = - b_m / STAR2(STAR2(a_m)); cp_m = (3.0 * STAR2(b_m) - a_m * c_m) / STAR7(a_m); // // just for speeding up the polynomial evaluation keep record of b/a and c/b // b_a_m = b_m / a_m; c_b_m = c_m / b_m; bp_ap_m = bp_m / ap_m; cp_bp_m = cp_m / bp_m; // // rho sampling as an expanded Chebyshev points set // rho_m = expandedChebyshevPoints(n); Rn_m.resize(rho_m.size()); spl_m.resize(rho_m.size()); // // Build interpolation table based on van Vleck curves. // Only negative points are actually computed, to minimize // calls to Rn(). The for-loop will completely cover the vector // of points only when their number is even, but if they are // an odd number then the cental point will be missing, but we now // that is at the center of the interval, that is, at 0.0 and, therefore, // the Rn value at that point is zero too. We take care of that condition // after the for-loop is ready. // For small correlation coefficients a polynomial approximation is // enough, for big correlation coefficients Rn is used for estimating // the quantized correlation function value. // for( unsigned int i = 0; i < rho_m.size() / 2; i++) { if ( rho_m[i] > -PolynomialAproximationLimit_m ) { EVAL_POLY(Rn_m[i], rho_m[i], a_m, b_a_m, c_b_m); } else { Rn_m[i] = Rn(l1_m, l2_m, rho_m[i], N_m); } Rn_m[rho_m.size() - 1 - i] = -Rn_m[i]; } if ( (rho_m.size() % 2) != 0 ) { Rn_m[(rho_m.size() - 1) / 2] = 0.0; } // // Compute 2nd derivatives, for spline interpolator // spline(Rn_m, rho_m, NatSpline, NatSpline, spl_m); } //------------------------------------------------------------------------------ QuantizationCorrection::~QuantizationCorrection() { } //------------------------------------------------------------------------------ QuantizationCorrection &QuantizationCorrection::operator =( const QuantizationCorrection &rhs) { if( this != &rhs ) { l1_m = rhs.l1_m; l2_m = rhs.l2_m; norm_m = rhs.norm_m; rho_m = rhs.rho_m; Rn_m = rhs.Rn_m; spl_m = rhs.spl_m; a_m = rhs.a_m; b_m = rhs.b_m; c_m = rhs.c_m; ap_m = rhs.ap_m; bp_m = rhs.bp_m; cp_m = rhs.cp_m; b_a_m = rhs.b_a_m; c_b_m = rhs.c_b_m; bp_ap_m = rhs.bp_ap_m; cp_bp_m = rhs.cp_bp_m; Vs_m = rhs.Vs_m; } return *this; } //------------------------------------------------------------------------------ void QuantizationCorrection::normalize(sdp_real Nm, sdp_real Vs, vector &L) { // // error handling kind of weak... // if ( fabs(Vs) < 1.0 ) { L = vector(L.size(), 9999.99); } // // lag normalizarion as per eqn. (2) (Gianni Comoretto) // // Rn(tau) = Nm ( L(tau) - Vs ) / Vs // // where 'n' could be 4, 8 or 16, and Nm adapts accordingly. // for ( vector::iterator l = L.begin(); l != L.end(); l++ ) { *l = Nm * (*l - Vs) / Vs; } } //------------------------------------------------------------------------------ sdp_real QuantizationCorrection::correct(const sdp_real &Rn) { sdp_real tmp; EVAL_POLY(tmp, Rn, ap_m, bp_ap_m, cp_bp_m); if( tmp < -PolynomialAproximationLimit_m || tmp > PolynomialAproximationLimit_m ) { tmp = splint(Rn_m, rho_m, spl_m, rho_m.size(), Rn); } return tmp * norm_m; } //------------------------------------------------------------------------------ void QuantizationCorrection::correct(vector &Rn) { for ( vector::iterator rn = Rn.begin(); rn != Rn.end(); rn++ ) { *rn = correct(*rn); } } //------------------------------------------------------------------------------ void QuantizationCorrection::correct(bool isAutoCorr, std::vector spectra) { static const sdp_real a = 0.2542; static const sdp_real b = 0.1134; // // As in section 9.4 Fixed correction // Note: more code but less time! // if ( isAutoCorr ) { for ( vector::iterator s = spectra.begin(); s != spectra.end(); s++ ) { #ifdef SDP_USE_FFTW3 *s *= a - b; #else s->re *= a - b; s->im *= a; /* b is a real number */ #endif } } else { for ( vector::iterator s = spectra.begin(); s != spectra.end(); s++ ) { #ifdef SDP_USE_FFTW3 *s *= a; #else s->re *= a; s->im *= a; #endif } } } //------------------------------------------------------------------------------ void QuantizationCorrection::correct4levels(vector &data) { static double r; static double r2; static double r4; static double r6; static double A = 1.1329552; static double B = 3.1056902; static double C = 2.9296994; static double D = 0.90122460; static double E = 2.7056559; static double F = 2.5012473; static double G = 0.73985978; // per LRD (from my note of long ago) lag 0 does not need correction // but here I am working with the list of lags that has been mirrored // around lag 0 and I am applying the correction to every lag in the list. for ( unsigned int i = 0; i < data.size(); i++ ) { r = (double)data[i]; r2 = r*r; r4 = r2*r2; r6 = r2*r4; data[i] = (sdp_real)( (A - B*r2 + C*r4 - D*r6) / (1 - E*r2 + F*r4 - G*r6) ) * r; } } //------------------------------------------------------------------------------ sdp_real QuantizationCorrection::Rn(sdp_real rho) { if ( rho < -PolynomialAproximationLimit_m || rho > PolynomialAproximationLimit_m ) { return Rn(l1_m, l2_m, rho, N_m); } return rho * a_m + rho * rho * rho * b_m + rho * rho * rho * rho * rho * c_m; } //------------------------------------------------------------------------------ // Private //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ double QuantizationCorrection::sigmaNlevels(double zho, int n) { if (n == 2) return 1.0; // meaningless if (n == 4) return M_SQRT1_2/inverse_erf((9.0-zho)/8.0); double tmp = thresholdStep(zho, n); if (tmp == 0.0) return 0.0; else return 1.0/tmp; } /* * Generic integral for correlator output for a n level correlator * with equispaced levels (quantization step = v1x, v1y) and * offset mux, muy, all in units of the respective RMS level * * rho is the correlation coefficient. * Works only for n even (for now) * Weights are +/-1, +/-3, .. +/-(n-1) */ //------------------------------------------------------------------------------ double QuantizationCorrection::Rn(double l1, double l2, double rho, int n) { double sum; #if 1 // This case is optimized for the n = 4 case as provided by Fred Schwab sum = 2.0 * (-4.0 + asin(rho)*M_1_PI + 3.0 * erf(l1*M_SQRT1_2) + erf(l2*M_SQRT1_2) + 4.0 * (L(0,l1,rho) + L(0,l2,rho) + L(l1, -l2, rho) + L(l1,l2,rho)) ); #else int i,j; int n2 = n/2; int n1 = -n2+1; double nd = double(n); double jl2; sum = -(n-1)*(n-1); for (j = n1; j < n2; ++j) { jl2 = double(j)*l2; sum += 4.0*L(0.0, jl2, rho); for (i = 1; i < n2; ++i) sum += 8.0*L(i*l1, jl2, rho); } for (i = 1; i < n2; ++i) { sum += 2.0*(nd-1.) * erf(i*l1_M_SQRT1_2); } #endif return sum; } //------------------------------------------------------------------------------ double QuantizationCorrection::bvnd(double dh, double dk, double r ) { const double ZERO=0.0, TWOPI=2.0*M_PI; const double SQRT_TWO_PI = sqrt(TWOPI); /* Gauss Legendre Points (X) and Weightsi (W): * First index selects number of points used in approximation: * 3 cases are for N = 6, N = 10 , N = 20 */ const static double X[3][10] = { { -0.9324695142031522e+00, -0.6612093864662647e+00, -0.2386191860831970e+00 }, { -0.9815606342467191e+00, -0.9041172563704750e+00, -0.7699026741943050e+00, -0.5873179542866171e+00, -0.3678314989981802e+00, -0.1252334085114692e+00 }, { -0.9931285991850949e+00, -0.9639719272779138e+00, -0.9122344282513259e+00, -0.8391169718222188e+00, -0.7463319064601508e+00, -0.6360536807265150e+00, -0.5108670019508271e+00, -0.3737060887154196e+00, -0.2277858511416451e+00, -0.7652652113349733e-01} }; const static double W[3][10] = { { 0.1713244923791705e+00, 0.3607615730481384e+00, 0.4679139345726904e+00 }, { 0.4717533638651177e-01, 0.1069393259953183e+00, 0.1600783285433464e+00, 0.2031674267230659e+00, 0.2334925365383547e+00, 0.2491470458134029e+00 }, { 0.1761400713915212e-01, 0.4060142980038694e-01, 0.6267204833410906e-01, 0.8327674157670475e-01, 0.1019301198172404e+00, 0.1181945319615184e+00, 0.1316886384491766e+00, 0.1420961093183821e+00, 0.1491729864726037e+00, 0.1527533871307259e+00 } }; int I, IS, LG, NG; double AS, A, B, C, D, RS, XS, BVN; double SN, ASR, H, K, BS, HS, HK; double fabsr = fabs(r); if( fabsr < 0.3 ) { NG = 0; LG = 3; } else if ( fabsr < 0.75 ) { NG = 1; LG = 6; } else { NG = 2; LG = 10; } H = dh; K = dk; HK = H*K; BVN = 0.0; if( fabsr < 0.925 ) { HS = ( H*H + K*K )/2.0; ASR = asin(r); for( I=0; I -100.0 ) { BVN = A*exp(ASR)*( 1.0 - C*( BS - AS )*( 1.0 - D*BS/5.0 )/3.0 + C*D*AS*AS/5.0 ); } if ( -HK < 100.0 ) { B = sqrt(BS); BVN -= exp( -HK/2.0 )*SQRT_TWO_PI*phid(-B/A)*B * ( 1.0 - C*BS*( 1.0 - D*BS/5.0 )/3.0 ); } A = A/2.0; for( I=0; I -100.0 ) { BVN = BVN + A*W[NG][I]*exp( ASR ) *( exp( -HK*( 1.0 - RS )/( 2.0*( 1.0 + RS ) ) )/RS - ( 1.0 + C*XS*( 1.0 + D*XS ) ) ); } } } BVN = -BVN/TWOPI; } if ( r > 0.0) BVN = BVN + phid( -MAX( H, K ) ); if ( r < 0.0 ) BVN = -BVN + MAX( ZERO, phid(-H) - phid(-K) ); } return BVN; } //------------------------------------------------------------------------------ double QuantizationCorrection::phid( double Z) { static const double P0 = 220.2068679123761, P1 = 221.2135961699311, P2 = 112.0792914978709, P3 = 33.91286607838300, P4 = 6.373962203531650, P5 = .7003830644436881, P6 = .03526249659989109; static const double Q0 = 440.4137358247522, Q1 = 793.8265125199484, Q2 = 637.3336333788311, Q3 = 296.5642487796737, Q4 = 86.78073220294608, Q5 = 16.06417757920695, Q6 = 1.755667163182642, Q7 = .08838834764831844; static const double ROOTPI = 2.506628274631001, CUTOFF = 7.071067811865475; double EXPNTL, ZABS, P; ZABS = fabs(Z); /* * |Z| > 37 */ if ( ZABS > 37.0 ) { P = 0.0; } else { /* * |Z| <= 37 */ EXPNTL = exp( -STAR2(ZABS)/2.0 ); /* * |Z| < CUTOFF = 10/sqrt(2) */ if ( ZABS < CUTOFF ) { P = EXPNTL*( ( ( ( ( ( P6*ZABS + P5 )*ZABS + P4 )*ZABS + P3 )*ZABS + P2 )*ZABS + P1 )*ZABS + P0 ) /( ( ( ( ( ( ( Q7*ZABS + Q6 )*ZABS + Q5 )*ZABS + Q4 )*ZABS + Q3 )*ZABS + Q2 )*ZABS + Q1 )*ZABS + Q0 ); /* * |Z| >= CUTOFF. */ } else { P = EXPNTL/( ZABS + 1.0/( ZABS + 2.0/( ZABS + 3.0/( ZABS + 4.0/( ZABS + 0.65 ) ) ) ) )/ROOTPI; } } if ( Z > 0.0 ) P = 1.0 - P; return(P); } //------------------------------------------------------------------------------ void QuantizationCorrection::spline(const vector &x, const vector &y, double yp1, double ypn, vector &y2) { int n = x.size(); double sig, p, qn, un; int i, k, i1, i_1; double x1, x2; vector u(n); if(yp1 > 0.99e30) { y2[0] = 0.0; u[0] = 0.0; } else { y2[0] = -0.5; u[0] = (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1); } int n_2 = n - 2; for(i=1; i <= n_2; i++) { i1 = i + 1; i_1 = i - 1; x1 = x[i] - x[i_1]; x2 = x[i1] - x[i_1]; sig = x1 / x2; p = sig*y2[i_1] + 2.0; y2[i] = (sig-1.0)/p; u[i] = (6.0*((y[i1]-y[i])/(x[i1]-x[i])-(y[i]-y[i_1])/x1)/x2-sig*u[i_1])/p; } int n_1 = n - 1; if (ypn > 0.99e30) { qn = 0.0; un = 0.0; } else { qn = 0.5; un = (3.0/(x[n_1]-x[n_2]))*(ypn-(y[n_1]-y[n_2])/(x[n_1]-x[n_2])); } y2[n_1] = (un-qn*u[n_2])/(qn*y2[n_2]+1.0); for(k = n-2; k >= 0; k--) y2[k] = y2[k]*y2[k+1] + u[k]; } //------------------------------------------------------------------------------ double QuantizationCorrection::splint(const vector &xa, const vector &ya, const vector &y2a, int n, double x) { double h, a, b, y; int k, klo, khi; klo = 0; khi = n-1; while( (khi - klo) > 1) { // This shift instead of divide save ~10% k = (khi+klo) >> 1; // k = (khi+klo)/2; if(xa[k] > x) khi = k; else klo = k; } h = xa[khi] - xa[klo]; if(h == 0.0) { //fprintf(stderr, "WARN Bad interpolation in routine splint"); return(x); } a = (xa[khi] - x)/h; b = (x - xa[klo])/h; y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi])*(h*h)/6.0; return(y); } //------------------------------------------------------------------------------ double QuantizationCorrection::inverse_erf( double x ) { static const double p1[] = {-13.0959967422,26.785225760,-9.289057635}; static const double q1[] = {-12.0749426297,30.960614529,-17.149977991,1.00000000}; static const double p2[] = {-.12402565221,1.0688059574,-1.9594556078,.4230581357}; static const double q2[] = {-.08827697997,.8900743359,-2.1757031196,1.0000000000}; static const double p3[] = {.1550470003116,1.382719649631,.690969348887, -1.128081391617, .680544246825,-.16444156791}; static const double q3[] = {.155024849822,1.385228141995,1.000000000000}; double ax, t, ret; ax = fabs(x); /* This approximation, taken from Table 10 of Blair et al., is valid for |x|<=0.75 and has a maximum relative error of 4.47 x 10^-8. */ if( ax <= 0.75 ) { t = x*x - 0.5625; ret = x*(p1[0] + t*(p1[1] + t*p1[2]))/ (q1[0] + t*(q1[1] + t*(q1[2] + t*q1[3]))); } else if( ax >= 0.75 && ax <= 0.9375 ) { /* This approximation, taken from Table 29 of Blair et al., is valid for .75<=|x|<=.9375 and has a maximum relative error of 4.17 x 10^-8. */ t = x*x - 0.87890625; ret = x*(p2[0] + t*(p2[1] + t*(p2[2] + t*p2[3])))/ (q2[0] + t*(q2[1] + t*(q2[2] + t*q2[3]))); } else if( ax >= 0.9375 && ax <= (1.0-1.0e-100) ) { /* This approximation, taken from Table 50 of Blair et al., is valid for .9375<=|x|<=1-10^-100 and has a maximum relative error of 2.45 x 10^-8. */ t = 1.0/sqrt(-log(1.0-ax)); ret = sign(x) * (p3[0]/t + p3[1] + t*(p3[2] + t*(p3[3] + t*(p3[4] + t*p3[5]))))/ (q3[0] + t*(q3[1] + t*(q3[2]))); } else { /* Finally, MyInverseErf[-1.]=-$MaxMachineNumber MyInverseErf[1.]=$MaxMachineNumber */ // // log something, but remember that we do not the context from where // we are being invoked. // ret = HUGE; } return(ret); } //------------------------------------------------------------------------------ double QuantizationCorrection::sign(double x) { if( x >= 0.0 ) return(1.0); else return(-1.0); } //------------------------------------------------------------------------------ double QuantizationCorrection::sign( double a, double b ) { if( b < 0.0 ) return( -fabs(a)); else return( fabs(a)); } // // This routine computes the threshold step in units of sigma, // given the zero lag ACF, for a generic N level quantization scheme. // N must be even (for now) and weights are +/-1, +/-3, .. +/-(N/2-1) // //------------------------------------------------------------------------------ double QuantizationCorrection::thresholdStep(double zerolag, int n) { double x, tol, f, fp, deltax, kd; double SQRT_05_PI = sqrt(2.0/M_PI); int itmax, k, i; x = 0.0; if( n % 2 == 0) x = 1.0; itmax = 30; tol = 1.0e-8; for( i = 1; i <= itmax; i++ ) { f = zerolag; fp = 0.0; if( n%2 == 1) { for( k = 1; k <= (n-1)/2; k++ ) { kd = (double)k; f = f-(2.0*kd-1.0)*erfc((2.0*kd-1.0)*x/M_SQRT2); fp = fp+SQRT_05_PI*STAR2(2.0*kd-1.0)*exp(-0.5*STAR2((2.0*kd-1.0)*x)); } } else { f = f - 1.0; for( k = 1; k <= (n-1)/2; k++ ) { kd = (double)k; f = f - 8*kd*erfc(kd*x/M_SQRT2); fp = fp + 8.0*STAR2(kd)*SQRT_05_PI*exp(-0.5*STAR2(kd*x)); } } deltax = -f/fp; deltax = sign(1.0,deltax)*MIN(0.5,fabs(deltax)); x += deltax; if (n % 2 == 1) x = MAX(0.0,x); if( fabs(deltax/x) < tol ) break; } return x; } //------------------------------------------------------------------------------ double QuantizationCorrection::inverse_cerf( double x ) { static const double p4[] = {.1550470003116,1.382719649631,.690969348887,-1.128081391617, .680544246825,-.16444156791}; static const double q4[] = {.155024849822,1.385228141995,1.000000000000}; static const double p5[] = {.00980456202915,.363667889171,.97302949837,-.5374947401}; static const double q5[] = {.00980451277802,.363699971544,1.000000000000}; double ret, t; // if( (x >= 0.0625) && (x <2.0) ) if( (0.0624 <= x) && (x< 2.0) ) { ret = inverse_erf( 1.0 - x ); } else if( (1.0e-100 <= x) && (x < 0.0625) ) { /* But for smaller values of x we'll use the approximations from Table 50 (again) of Blair et al., as well as Table 70. */ t = 1.0/sqrt(-log(x)); ret = (p4[0]/t + p4[1] + t*(p4[2] + t*(p4[3] + t*(p4[4] + t*p4[5]))))/ (q4[0] + t*(q4[1] + t*(q4[2]))); } //else if( x < 1.0e-100 && x > 1.0e-1000 ) else if( (1.0e-1000 < x) && (x < 1.0e-100) ) { /* This approximation, taken from Table 70 of Blair et al., is valid for 10^-1000<=|x|<=10^-100 and has a maximum relative error of 2.45 x 10^-8. */ t = 1.0/sqrt(-log(x)); ret = (p5[0]/t + p5[1] + t*(p5[2] + t*p5[3]))/(q5[0] + t*(q5[1] + t*(q5[2]))); } else { /* Finally, MyInverseErfc[0.]=$MaxMachineNumber MyInverseErfc[2.]=-$MaxMachineNumber */ ret = HUGE; } return ret; } //------------------------------------------------------------------------------ vector QuantizationCorrection::expandedChebyshevPoints(unsigned int n) { double PI_2FitPts = M_PI / ((double)2.0 * (double)(n + 1)); double cosPI_2FitPts = cos(PI_2FitPts); vector pts(n + 1); // // if n is even then take advantage of the symettry and the fact that // the central point is zero. Otherwise, compute analytically each point. // if ( n % 2 == 0 ) { for( unsigned int i = 0; i < (n + 1) / 2; i++) { pts[i]= -cos((2 * i + 1.0) * PI_2FitPts) / cosPI_2FitPts; pts[n - i] = -pts[i]; } pts[n / 2] = 0.0; } else { for ( unsigned int i = 0; i < n + 1; i++ ) { pts[i] = -cos((2 * i + 1) * PI_2FitPts) / cosPI_2FitPts; } } return pts; } //------------------------------------------------------------------------------ double QuantizationCorrection::Ck(int N, int k, double l) { if ( k != 1 && k != 3 && k != 5 ) { return 0.0; } int lim = N / 2 - 1; double c = 0; switch ( k ) { case 1: for ( int i = 1; i <= lim; i++ ) { c += exp(-STAR2(i * l) / 2); } c *= 2.0; c += 1.0; break; case 3: for ( int i = 1; i <= lim; i++ ) { c += (1.0 - STAR2(i * l)) * exp(-STAR2(i * l) / 2); } c *= 2.0; c += 1.0; break; case 5: for ( int i = 1; i <= lim; i++ ) { c += (3.0 - 6.0 * STAR2(i * l) + STAR2(STAR2(i * l))) * exp(-STAR2(i * l) / 2); } c *= 2.0; c += 3.0; break; default: return 0.0; } return c; } /*___oOo___*/