#ifndef QUANTIZATION_CORRECTION_H #define QUANTIZATION_CORRECTION_H /* @(#) $Id: QuantizationCorrection.h,v 1.3 2009/12/11 22:22:45 ramestic Exp $ * * Copyright (C) 2004 * 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 */ /** 3 and 9 level van vleck correction based on Fred Schwab's excellent paper: DRAFT, Van Vleck Correction for the GBT correlator, Feb 11, 2002 Fred Schwab (NRAO) These routines were proven by Bill Sisk (NAIC) bsisk@naic.edu These routines were converted to C by Jeff Hagen (NAIC) jeffh@naic.edu Modification for N levels by G. Comoretto (Arcetri) comore@arcetri.astro.it >From G. Comorreto's email of Nov 10, 2004: Dear Jim, I had a look at the software from Schwab. I send you the C++ file, and the include file with the function prototypes (plus those for the 3 bit routines by Schwab). I wrote a generic routine for Van Vleck correction for any quantization scheme with an even number of levels, and tested it against mine. As the original Schwab routine, it allows for the quantized signal to have nonzero mean and arbitrary RMS. The output INCLUDES multiplication for sigma1*sigma2, so that, after correction, zero delay lag is sigma1*sigma2, not 1.0. I wrote a small routine to derive the RMS value from the zero delay lag. I have tested it, and the overall accuracy is ~1e-7 for sigma in the range (0.2-2). I am now working on accuracy tests. I am trying direct tests on simulated data, before and after quantization. The overall accuracy of my routine is better, especially at rho near 1, but both are around 10-6. I do not know which is faster, I am testing this also. In the process I doscovered a small error on my routine, if you plan to use it I will send you the correct version. My routine, however, does not include at the moment correction for nonzero signal mean. I would therefore use the new one, even if I do not know how to determine the signal mean. The following fragment of code can show you how to use the class. stdev1=VanVleck::sigmaNlev(acf1[0]); stdev2=VanVleck::sigmaNlev(acf2[0]); VanVleck corr(4, stdev1,stdev2); while (correction_is_valid) { acquire(cf_data); corr.Correct(cf_data, npts); .... } */ #include #include #include // // Local stuff // #include "FFTDefs.h" #define L(x,y,z) bvnd(x,y,z) class QuantizationCorrection { public: /** Default ctor, everything defaulted to nothing */ QuantizationCorrection(); /** Parameterized ctor ** This ctor combines G. Comoretto's calls to sigmaNlev(acf1) & ** sigmaNlev(acf2) with the constructor eliminating the extra external ** calls. ** Note: according to F. Schwab, 16 intepolating intervals should provide ** sufficient number of points. Originally, this was 64. This decreases ** the ctor time by a factor of ~4. ** @param N quantizer levels ** @param n number of interpolating intervals over the van Vleck curve. ** @param sigma1 antenna i RMS of N-bit quantized signal ** @param sigma2 antenna j RMS of N-bit quantized signal ** @param Vs correlation offset */ QuantizationCorrection(unsigned int quantizerLevels, unsigned int interpolationIntervals, sdp_real sigma1, sdp_real sigma2, sdp_real Vs); /// copy constructor QuantizationCorrection( const QuantizationCorrection &toCopy) { *this = toCopy; } /// assignment operator QuantizationCorrection &operator =( const QuantizationCorrection &); /// destructor ~QuantizationCorrection(); /** Normalize lags as per equation (2) in Comoretto's paper. ** @param Nm scale factor, a function of the bits mode (9 for 2x2 and ** 225 for 3x3 and 4x4 bits modes. ** @param Vs correlation offset. ** @param raw lags as read at correlator output. */ static void normalize(sdp_real Nm, sdp_real Vs, std::vector &L); /** Correct quantization effects on received correlation function. ** The algorithm applies a polynomial approximation for small ** correlation values and a spline interpolation over van Vleck ** points when the correlation coefficient is big. */ sdp_real correct(const sdp_real &Rn); void correct(std::vector &cf); /** Correct quantization effects introduced by the analog sample. */ static void correct(bool isAutoCorr, std::vector spectra); /** The two bit, four level quantization correction used comes from ** VLBA Memo No. 75, by Fred Schwab, November 19, 1986. ** ** Equation (10) is used, for the case of weight = 3.3358750 and sampler ** threshold = 0.98159883. ** ** The equation applied to each "lag" (the results of the inverse transform ** of the power spectra) is: ** ** A B C D ** 1.1329552 - 3.1056902r^2 + 2.9296994r^4 - 0.90122460r^6 ** corrected = ------------------------------------------------------- * r ** 1 - 2.7056559r^2 + 2.5012473r^4 - 0.73985978r^6 ** E F G ** ** where 'r' is the kag value. */ static void correct4levels(std::vector &data); /** Returns the number of quantization levels the class was built for. */ unsigned int getLevels() { return N_m; } /** Returns the number of interpolation intervals applied over the van Vleck ** curve. */ unsigned int getVanVleckInterpolatingIntervals() { return rho_m.size(); } /** ** Useful function: derive sigma from zero lag ACF ** Converts zho (zero lag ACF) into an unbiased estimate of input RMS (sigma) ** output is unitless: (sigma*vstep) ** 4 level case is set apart for efficiency ** n must be > 2 and even */ static double sigmaNlevels(double zho, int n); /* * Compute threshold step in units of sigma, given zero lag ACF * Used by sigmaNlev; */ static double thresholdStep(double zerolag, int n); sdp_real Rn(sdp_real rho); private: // // Limit between polynomial approximation and spline interpolation. // static const double PolynomialAproximationLimit_m; // // ??? // static const double NatSpline = 1.0e30; // // Number of levels to correct for. This variable is set at construction // time and it reflects the number of levels the spline table was built // for. // unsigned int N_m; // // normalized quantization steps // sdp_real l1_m, l2_m; // // TFB processing must be such that the quantizer correction should output // lags normalized to the signal amplitude (s1 * s2), such that further // processing across sub-bands render coherent results. The following // variable records the s1*s2 value as computed from lag zero using equation // (4) to compute sigma. // sdp_real norm_m; // // Rn and pho as defined in Comoretto's paper // Note: rho is then filled with a sequence of expanded Chebyshev points, // assuming the standard interval [-1, 1] and FitPts_m points, or // FitPts_m-1 intervals. // std::vector Rn_m; std::vector rho_m; // // polynomial aproximation of the relation between Rn and rho. // a, b and c as in section 5.2.1, the specific expresion is right // before equation (5); ap_m, bp_m and cp_m ('p' for prime) represent // the polynomial coefficients in equation (8). // b_a, c_b, bp_ap and cp_bp are the quotients b/a and c/b, respectively. // Recording these coefficient quotients is convinient because evaluating // the polynomial becomes more efficient when expresed the following way: // // x = a x + b x^3 + c x5 = ax (1 + (b/a) x^2 (1 + (c/b) x^2)) // double a_m, b_m, c_m; double b_a_m, c_b_m; double ap_m, bp_m, cp_m; double bp_ap_m, cp_bp_m; // // spline ... // std::vector spl_m; sdp_real Vs_m; //! Vs /** * Generic integral for correlator output for a n level correlator * with equispaced levels (quantization step = v1x, v1y) and * offset mux, muy * v1x, v1y, mux, muy are expressed in units of the signal RMS * rho is the correlation coefficient. * Works only for n even (for now) * Weights are +/-1, +/-3, .. +/-(n-1) */ double Rn(double l1, double l2, double rho, int n); /** * A function for computing bivariate normal probabilities. * * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenzwsu.edu * * This function is based on the method described by * Drezner, Z and G.O. Wesolowsky, (1989), * On the computation of the bivariate normal inegral, * Journal of Statist. Comput. Simul. 35, pp. 101-107, * with major modifications for double precision, and for |R| close to 1. * * BVND - calculate the probability that X is larger than DH and Y is * larger than DK. * * @param dh DOUBLE PRECISION, integration limit * @param dk DOUBLE PRECISION, integration limit * @param r DOUBLE PRECISION, correlation coefficient * */ static double bvnd(double dh, double dk, double r ); /** * Normal distribution probabilities accurate to 1.e-15. * Z = no. of standard deviations from the mean. * * Based upon algorithm 5666 for the error function, from: * Hart, J.F. et al, 'Computer Approximations', Wiley 1968 * * Programmer: Alan Miller * Latest revision - 30 March 1986 */ static double phid( double Z); /** ** This subroutine calculates the second derivatives of the ** interpolating cubic spline function. It was taken from ** "Numerical Recipes," by Press et al. 1986, p. 88. ** removed fortran stupidity - now 0 zero based arrays ** @param x, y are the arrays of input interpolation values; ** @param yp1, ypn are the first derivatives at the first and last points, resp. They can be set to 1.E30 for a natural spline; ** @param y2 is the output array of second derivatives to be used by SPLINT. */ static void spline(const std::vector &x, const std::vector &y, double yp1, double ypn, std::vector &y2); /** ** This subroutine performs a cubic spline interpolation. It was ** taken from "Numerical Recipes," by Press et al. 1986, p. 88. ** remove fortran stupidity: return double as answer and make arrays zero based ** Modified 20 Oct 89 to change x**3 to x*x*x ** Mod. 25 Oct 89 to take out PAUSE and replace with call to WARNING ** @param xa, ya are the arrays of input interpolation data; ** @param y2a is the array of second derivatives produced by SPLINE; ** @param n is the number of input interpolation points; ** @param x is the X value for which we wish an interpolation; ** @param y is the output interpolation value. */ static double splint(const std::vector &xa, const std::vector &ya, const std::vector &y2a, int n, double x); /** inverse cerf from Fred Schwab The following approximations to the inverse of the error function are taken from J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev Approximations for the Inverse of the Error Function", Mathematics of Computation, 30 (1976) 827-830 + microfiche appendix. */ static double inverse_erf( double x ); /** Now, for the inverse of the complementary error function we'll generally use the relation InverseErfc[x]=InverseErf[1-x] for larger values of x, those satisfying .0625<=x<2 : JAP - I don't think that this is used */ static double inverse_cerf( double x ); /// return sign (+1 or -1) of x. static double sign(double x); /// transport sign of b on a static double sign(double,double); /** Build an expanded Chebyshev sequence of points. ** The general expression for the expanded Chebyshev set of points is: ** ** (a+b) / 2 + ((a-b)/2) cos((2i+1) pi / (2n+2)) / cos(pi / (2n+2)) ** ** in our case the inteval is always [-1,1], that is, b=1 and a=-1. ** @param n number of intervals, therefore, n+1 is the number of points. */ static std::vector expandedChebyshevPoints(unsigned int n); /** Compute polynomial approximation coeffient. ** This method computes Ck(l) as introduced in equations (10), (11) ** and (12) of Comoretto's paper. The idea is that for low correlation ** coefficients (fabs(rho) < 0.2 ) then a 5th order polynom (odd powers ** only) is enough for approximating representing the R4(l1, l2, rho) ** function. ** @param N quantization levels ** @param k polynomial power (1, 3 or 5) ** @param l normalized threshold spacing (l = 1 / sigma) */ static double Ck(int N, int k, double l); }; #endif