/***************************************************************************/ /* */ /* CombinedEDF() */ /* */ /* Function to calc the estimated # of chi-squared degrees of freedom */ /* for the normal and overlapping, modified and unmodified, Allan and */ /* Hadamard variances of a power law noise process. */ /* */ /* Parameters: int N = # phase data points at tau0 */ /* int m = averaging factor = tau/tau0 */ /* int a = alpha (-4 to +2) */ /* int d = order of phase difference */ /* 1 = first difference (not used) */ /* 2 = Allan variance */ /* 3 = Hadamard variance */ /* int F = filter factor */ /* 1 = modified */ /* m = unmodified */ /* int S = stride factor */ /* 1 = long (tau) */ /* m = short (tau0) */ /* int v = calc version type (0=simple, 1=full) */ /* */ /* Return: float = edf (or -1 if error) */ /* */ /* References: 1. C. Greenhall, "Normalized Uncertainty of */ /* Difference Variances", June 4, 2003 */ /* (private communication via e-mail) */ /* */ /* 2. C. Greenhall and W. Riley, "Uncertainty of */ /* Stability Variances Based on Finite */ /* Differences", Proc. 2003 PTTI Meeting, */ /* December 2003 (to be published). */ /* */ /* Notes: 1. Applies for alpha between -4 and 2 (RR FM to */ /* W PM) for long or short stride (tau or tau0), */ /* modified or unmodified Allan or Hadamard */ /* variances. Does not apply to Total variances */ /* or Thêo1. */ /* 2. Two algorithms are used, a simplified version */ /* using a truncated summation, and a full version */ /* using 4 cases and 3 tables. */ /* */ /* Revision Record: */ /* 06/08/03 Created */ /* 06/17/03 Running */ /* 06/23/03 Added full version */ /* 06/24/03 Full version running */ /* 06/26/03 Looking pretty good */ /* 06/27/03 Added calc version type parameter */ /* Added calc parameter recording */ /* 07/03/03 Made corrections per C. Greenhall e-mail of 7/1/03 */ /* 07/14/03 Brought into Stable_r.c from edf4.c */ /* 07/16/03 Made corrections per C/ Greenhall e-mail of 7/9/03 */ /* 07/18/03 This version created to accompany PTTI'03 paper */ /* 11/20/03 Replaced Numerical Recipes code with new bico() */ /* */ /* (c) Copyright 2003 W. J. Riley, Hamilton Technical Services, */ /* 195 Woodbury St., S. Hamilton, MA 01982 USA, edf@wriley.com */ /* */ /* This source code may be freely used without restriction provided that */ /* this copyright notice is included with it. No warranty of any kind is */ /* made. Use this code at your own risk. */ /* */ /* This source file should compile portably without errors as-is under */ /* all standard C compilers It is suggested that a simmple command line */ /* test driver suitable for your system be written to exercise this */ /* function and verify its proper operation before use. A good first */ /* test is to check it against the example in the paper. */ /* */ /***************************************************************************/ // Includes #include // Macros #ifndef min #define min(a,b) (((a)<(b))?(a):(b)) #endif // Combined EDF function float CombinedEDF(int N, int m, int a, int d, int F, int S, int v) { // Notes re function inputs: // N is Greenhall Nx // m is Greenhall m // a is Greenhall alpha // d is Greenhall d // F is Greenhall F // S is Greenhall S // Supporting function prototypes // Notes re functions: // Function names are capitalized // Functions Sw(), Sx(), and Sz() are same as Greehall names // Function Sum is Greenhall BasicSum() double Sw(double t, int a); double Sx(double t, double F, int a); double Sz(double t, double F, int a, int d); double Sum(int J, double M, double S, double F, int a, int d); // Log function to handle zero double Log(double x); // Binomial coefficient function double bico(int n, int k); // Local variables int L; // Same as Greenhall L double M; // Same as Greenhall M int J; // Summation limit int Jmax=100; // Max J to use summation int k; // Summation index int K; // ceil(r) double r; // Same as Greenhall r=M/S double s=(double)S; // double version of S double a0; // Same as Greenhall a0 (from Table 2) double a1; // Same as Greenhall a1 (from table 2) double b0; // Same as Greenhall b0 (from Table 2) double b1; // Same as Greenhall b1 (from table 2) double edf; // Same as Greenhall edf // Lookup tables // Table 1: Coefficients for modified variances // a0 & a1 as a function of a (row) & d (column) // d=1 d=2 d=3 // a0 a1 a0 a1 a0 a1 double T1[7][6]={ {0.667, 0.333, 0.778, 0.500, 0.880, 0.667}, // a=+2 {0.840, 0.345, 0.997, 0.616, 1.141, 0.843}, // a=+1 {1.079, 0.368, 1.033, 0.607, 1.184, 0.848}, // a= 0 {0.000, 0.000, 1.048, 0.534, 1.180, 0.816}, // a=-1 {0.000, 0.000, 1.302, 0.535, 1.175, 0.777}, // a=-2 {0.000, 0.000, 0.000, 0.000, 1.194, 0.703}, // a=-3 {0.000, 0.000, 0.000, 0.000, 1.489, 0.702}}; // a=-4 // Access as T1[2-a][i], where a=alpha // and i=2*(d-1)+j and j=a# // Table 2: Coefficients for unmodified variances // a0 & a1 as a function of a (row) & d (column) // d=1 d=2 d=3 // a0 a1 a0 a1 a0 a1 double T2[7][6]={ {1.500, 0.500, 1.944, 1.000, 2.310, 1.500}, // a=+2 {78.60, 25.20, 790.0, 410.0, 9950., 6520.}, // a=+1 {0.667, 0.167, 0.667, 0.333, 0.778, 0.500}, // a= 0 {0.000, 0.000, 0.852, 0.375, 0.997, 0.617}, // a=-1 {0.000, 0.000, 1.079, 0.368, 1.033, 0.607}, // a=-2 {0.000, 0.000, 0.000, 0.000, 1.053, 0.553}, // a=-3 {0.000, 0.000, 0.000, 0.000, 1.302, 0.535}}; // a=-4 // Access as T2[2-a][i], where a=alpha // and i=2*(d-1)+j and j=a# // Table 3: Coefficients for log denominator, // unmodified variances, F FM // b0 & b1 as a function of b (row) & d (column) // d=1 d=2 d=3 // b0 b1 b0 b1 b0 b1 double T3[6]= {6.000, 4.000, 15.23, 12.00, 47.80, 40.00}; // Access as T3[2*(d-1)+j], where j=b# // Check arguments // d must be 1, 2 or 3 if( (d<1) || (d>3)) { return -1; // Error } // F must be either 1 or m if( (F!=1) && (F!=m) ) { return -1; // Error } // m must be >=1 if(m<1) { return -1; // Error } // a must be between -4 and +2 if( (a<-4) || (a>2) ) { return -1; // Error } // S must be either 1 or m if( (S!=1) && (S!=m) ) { return -1; // Error } // Check for legal alpha: a+2d>1 if(a+2*d<=1) { return -1; // Error - Illegal alpha } // Find # summands L=(m/F)+m*d; // L is always an int (F is either 1 or m) // Calc M - not necessarily an int M=(double)(S); M*=(double)(N-L); M/=(double)m; M=floor(M); M+=1.0; J=min((int)M, (d+1)*S); // Check # data points if(N=d+1) { // Get a0 & a1 from Table 1 a0=T1[2-a][2*(d-1)+0]; a1=T1[2-a][2*(d-1)+1]; // Calc edf edf=((1.0/r)*(a0-(a1/r))); edf=(1.0/edf); } else { // Calc edf edf=(Sz(0, 1, a, d)*Sz(0, 1, a, d)*Jmax/ Sum(Jmax, Jmax, (double)Jmax/r, 1, a, d)); } } return (float)edf; // EDF for Case 1 } else // Unmodified variances: F=m { if(a<=0) // Case 2. W FM to RR FM { if(J<=Jmax) { if(m*(d+1)<=Jmax) // m'=m; { // Calc edf edf=(Sz(0, m, a, d)*Sz(0, m, a, d)*M/ Sum(J, M, S, m, a, d)); } else // m'=infinity, use F=m=0 as flag { edf=(Sz(0, 0, a, d)*Sz(0, 0, a, d)*M/ Sum(J, M, S, 0, a, d)); } } else // J>Jmax { if(r>=d+1) { // Get a0 & a1 from Table 2 a0=T2[2-a][2*(d-1)+0]; a1=T2[2-a][2*(d-1)+1]; // Calc edf edf=((1.0/r)*(a0-(a1/r))); edf=(1.0/edf); } else // r=d+1) { // Get a0 & a1 from Table 2 a0=T2[2-a][2*(d-1)+0]; a1=T2[2-a][2*(d-1)+1]; // Get b0 & b1 from Table 3 b0=T3[2*(d-1)+0]; b1=T3[2*(d-1)+1]; // Calc edf edf=(((b0+b1*Log(m))*(b0+b1*Log(m))*r)/ (a0-(a1/r))); } else { // Get b0 & b1 from Table 3 b0=T3[2*(d-1)+0]; b1=T3[2*(d-1)+1]; // Calc edf edf=(((b0+b1*Log(m))*(b0+b1*Log(m))*Jmax)/ Sum(Jmax, Jmax, (double)Jmax/r, (double)Jmax/r, 1, d)); } } return (float)edf; // EDF for Case 3 } else // Case 4. W PM a=2 { K=(int)ceil(r); if(K<=d) { // Use a0 and a1 as working variables a0=bico(2*d, d); // Calc sum edf=0.0; for(k=1; k<=K-1; k++) { a1=bico(2*d, d-k); edf+=((1.0-((double)k/r))*a1*a1); } // Complete edf calc edf*=(2.0/(a0*a0)); edf+=1.0; edf=M/edf; } else { // Get a0 & a1 for a=2 from Table 2 a0=T2[0][2*(d-1)+0]; a1=T2[0][2*(d-1)+1]; // Calc edf edf=(a0-(a1/r)); edf=M/edf; } return (float)edf; // EDF for Case 4 } } } // Supporting functions double Sw(double t, int a) { // Local variables double b; double sw; // Calc abs value b=fabs(t); // Calc sw switch(a) { case 2: { sw=-b; } break; case 1: { sw=t*t*Log(b); } break; case 0: { sw=b*b*b; } break; case -1: { sw=-t*t*t*t*Log(b); } break; case -2: { sw=-b*b*b*b*b; } break; case -3: { sw=t*t*t*t*t*t*Log(b); } break; case -4: { sw=b*b*b*b*b*b*b; } break; } return sw; } double Sx(double t, double F, int a) { // Local variables double sx; // Calc sx if(F>0) { sx=F*F*(2*Sw(t, a) - Sw(t-1.0/F, a) - Sw(t+1.0/F, a)); } else // F=0 is flag for F=infinity { // Note: This case applies only for a<=0 sx=Sw(t, a+2); } return sx; } double Sz(double t, double F, int a, int d) { // Local variables double sz; // Calc sz switch(d) { case 1: { sz=2*Sx(t, F, a) - Sx(t-1, F, a) - Sx(t+1, F, a); } break; case 2: { sz=6*Sx(t, F, a) - 4*Sx(t-1, F, a) - 4*Sx(t+1, F, a) + Sx(t-2, F, a) + Sx(t+2, F, a); } break; case 3: { sz=20*Sx(t, F, a) - 15*Sx(t-1, F, a) - 15*Sx(t+1, F, a) + 6*Sx(t-2, F, a) + 6*Sx(t+2, F, a) - Sx(t-3, F, a) - Sx(t+3, F, a); } break; } return sz; } double Sum(int J, double M, double S, double F, int a, int d) { // Local variables int j; double sum; double z; // Initialize sum sum=Sz(0, F, a, d); sum*=sum; // Calc sum for(j=1; j<=J-1; j++) { z=Sz(j/S, F, a, d); z*=z; z*=(1-j/M); sum+=2*z; } // Add last term z=Sz(J/S, F, a, d); z*=z; z*=(1-J/M); sum+=z; return sum; } // Wrapper function for log(x) modified to return 0 for x=0 double Log(double x) { if(x) { return log(x); } else { return 0; } } // Function to find binomial coefficients // OK for small integers // Ref: C. Greenhall e-mail of 10/14/03 double bico(int n, int k) { // Local variables int i; // Index double b; // Binominal coefficient // Compute binomial coefficient b=(n,k) b=1.0; for(i=0; i