# HG changeset patch
# User Daniele Nicolodi <daniele.nicolodi@obspm.fr>
# Date 1390394467 -3600
# Node ID 792ac7151f0f67e53a2b83436b7f74e778c5b5dd
# Parent  7b9cf3d4346e5cc8ae7a84f427cc9e1e3ed7a9ff
Fix Allan deviation linear drift contribution removal

diff -r 7b9cf3d4346e -r 792ac7151f0f Allan.c
--- a/Allan.c	Wed Jan 22 12:29:58 2014 +0100
+++ b/Allan.c	Wed Jan 22 13:41:07 2014 +0100
@@ -68,8 +68,6 @@
 			if (N > 1) {
 				Instance->AllanVar[i] = (Instance->AllanVar[i]*(N-2)
 										+0.5*pow(CurrentMean-LastMean,2))/(N-1) ;
-			} else { // ie if N=1, which is realized for the first completed block
-				Instance->FirstMean[i] = CurrentMean;
 			}
 			Instance->LastMean[i] = CurrentMean;
 			Instance->BlocksNumber[i] = N;
@@ -85,7 +83,6 @@
 static void Allan_Reset(Allan_Data *Instance)
 {
 	memset(Instance->AllanVar, 0, sizeof(Instance->AllanVar));
-	memset(Instance->FirstMean, 0, sizeof(Instance->FirstMean));
 	memset(Instance->LastMean, 0, sizeof(Instance->LastMean));
 	memset(Instance->BlocksNumber, 0, sizeof(Instance->BlocksNumber));
 	memset(Instance->CurrentAverage, 0, sizeof(Instance->CurrentAverage));
@@ -97,32 +94,29 @@
 	
 static void Allan_Display(Allan_Data *Instance)
 {
-	int i, N;
+	int i;
 	double X[ALLAN_NUM_DATAPOINTS];
 	double Y[ALLAN_NUM_DATAPOINTS];
-	double Normalizer = Instance->normalization;
-	double DriftTau, FirstFreq, LastFreq, Error;
+	double Error;
 	int dedrift;
 	
 	GetCtrlVal(Instance->AllanPanel, ALLANPANEL_DEDRIFT, &dedrift);
 	SetCtrlVal(Instance->AllanPanel, ALLANPANEL_DRIFT, Instance->Drift);
 	
 	for (i = 0; i < ALLAN_NUM_DATAPOINTS; i++) {
+		if (Instance->AllanVar[i] == 0.0)
+			break;
 		X[i] = pow(2, i);
+		/* remove  linear drift estimated contribution to the allan variance */
 		if (dedrift) {
-			DriftTau = Instance->Drift * X[i];
-			N = Instance->BlocksNumber[i];
-			if (N > 1) { // if enough data to calculate AVAR for tau = 2^i*tau0, and therefore, AllanVar[i] is meaningful
-				FirstFreq = Instance->FirstMean[i];
-				LastFreq = Instance->LastMean[i];
-				Y[i] = sqrt(Instance->AllanVar[i]+DriftTau*(DriftTau/2-(LastFreq-FirstFreq)/(N-1)))/Normalizer;
-				// Not totaly sure it works... to be checked thoroughly   
-				}
-			else {
-				Y[i] = 0;
-			}
+			double corr = 0.5 * (Instance->Drift * X[i]) * (Instance->Drift * X[i]);
+			/* uncertainty in the estimate of the drift rate may over correct
+			   the estimated allan variance and result in a negative value. in
+			   this case we take the absolute value of the corrected value as
+			   the best estimator of the allan variance */
+			Y[i] = sqrt(fabs(Instance->AllanVar[i] - corr)) / Instance->normalization;
 		} else {
-			Y[i] = sqrt(Instance->AllanVar[i])/Normalizer;
+			Y[i] = sqrt(Instance->AllanVar[i]) / Instance->normalization;
 		}
 	}
 	
diff -r 7b9cf3d4346e -r 792ac7151f0f Allan.h
--- a/Allan.h	Wed Jan 22 12:29:58 2014 +0100
+++ b/Allan.h	Wed Jan 22 13:41:07 2014 +0100
@@ -11,7 +11,6 @@
 	int AllanPanel;									/* display panel handle */
 	double X[ALLAN_NUM_DATAPOINTS];
 	double AllanVar[ALLAN_NUM_DATAPOINTS];			/* Allan variance computed at tau = 2^(N tau0) */
-	double FirstMean[ALLAN_NUM_DATAPOINTS];			/* mean of the first block of 2^N data for which averaging was completed */
 	double LastMean[ALLAN_NUM_DATAPOINTS];			/* mean of the last block of 2^N data for which averaging was completed */
 	int BlocksNumber[ALLAN_NUM_DATAPOINTS];			/* number of 2^N data points blocks used for computing AVAR(2^(N tau0)) */
 	double CurrentAverage[ALLAN_NUM_DATAPOINTS];  	/* sum of the last P data points retrieved. P < 2^N */