Mercurial > hg > fxanalyse

changeset 73:1dfaf5ef0352

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor mean and slope measurement
author Daniele Nicolodi <daniele.nicolodi@obspm.fr>
date Thu, 08 Nov 2012 13:05:55 +0100
parents 8e6a2ec85169
children 9f81652432d0
files FXAnalyse.c
diffstat 1 files changed, 72 insertions(+), 130 deletions(-) [+]
line wrap: on
line diff
--- a/FXAnalyse.c	Mon Oct 29 16:13:27 2012 +0100
+++ b/FXAnalyse.c	Thu Nov 08 13:05:55 2012 +0100
@@ -69,8 +69,9 @@
 int Measuring_3 = FALSE;
 
 double FrequDDS1=110000000.0, FrequDDS4=110000000.0;
-double Slope_1=0.0,Slope_2=0.0,Slope_3=0.0,Beatslope_2=0.0; 
+double Slope_1=0.0,Slope_2=0.0,Slope_3=0.0,Beatslope_2=0.0;
 double SlopeTime1=40.0, SlopeTime2=40.0; SlopeTime3=40.0;
+double Ch4Slope = 0.0;
 
 double N_1=0.0, N_2=0.0, N_3=0.0;
 double DeltaT_1=20.0, DeltakHz_1=500.0, t1_1=0.0, t2_1=0.0, t3_1=0.0, Frepplus_1=0.0, Frepminus_1=0.0;
@@ -78,10 +79,6 @@
 double DeltaT_3=20.0, DeltakHz_3=500.0, t1_3=0.0, t2_3=0.0, t3_3=0.0, Frepplus_3=0.0, Frepminus_3=0.0;
 
 int n_1=0, n_2=0, n_3=0;
-double Frequ_slope_1=0.0,Moy_slope_1=0.0,Slope_slope_1=0.0,Frequ_slope_2=0.0,Moy_slope_2=0.0,Slope_slope_2=0.0,Frequ_slope_3=0.0,Moy_slope_3=0.0,Slope_slope_3=0.0;  
-int N_slope_1=0,N_slope_2=0,N_slope_3=0;
-double  Beat_slope_2=0.0 ,Moy_Beatslope_2=0.0,Slope_Beatslope_2=0.0; 
-double Ch4_slope=0.0,Moy_Ch4slope_1=0.0,Slope_Ch4slope_1=0.0,Ch4Slope=0.0;
 
 double FrequencyDDSBesInit = 0.0;
 double FrequencyDDS3Init = 0.0;
@@ -106,14 +103,8 @@
 int StopSlopeCancellingOnUnlocked = TRUE;
 double TimetoSlope = 60.0;
 double SlopeMeasuringTimeBegin = 0.0;
+double appliedSlope = 0.0; // currently applied frequency dedrifiting slope
 
-int Nsamples = 0;				// number of samples in current measurement
-double previousFreq = 0.0;		// previous frequency value
-double meanFreq = 0.0;			// mean value
-double measuredSlope = 0.0;		// measured slope
-double appliedSlope = 0.0;		// currently applied frequency dedrifiting slope
-
-double LimitToDelock=5.0;
 double limitotakoff=70.0;
 
 int ratio=10; //Recentre la frequence tous les ratios
@@ -138,6 +129,36 @@
 double CenteringTimeBegin10K = 0.0;
 
 
+struct stat {
+	int samples;
+	double mean;
+	double slope;
+	double previous;
+};
+
+void stat_zero(struct stat *s)
+{
+	s->samples = 0;
+	s->mean = 0.0;
+	s->slope = 0.0;
+	s->previous = 0.0;
+}
+
+void stat_accumulate(struct stat *s, double value)
+{
+	s->samples += 1;
+	
+	if (s->samples > 1)
+		s->slope = (s->slope * (s->samples - 2) + 6 * (value - s->mean) / s->samples) / (s->samples + 1);
+					
+	s->mean = ((s->samples - 1) * s->mean + value) / s->samples;
+	
+	s->previous = value;
+}
+
+struct stat stat_math1, stat_ch2, stat_ch4, freq;
+
+
 muParserHandle_t initMathParser() 
 {
 	muParserHandle_t parser = mupCreate();
@@ -673,11 +694,9 @@
 						GetCtrlVal(MainPanel, PANEL_DDS2, &FrequencyDDSBesInit);
 						t2_1 = t3_1 = 0.0;
 						t1_1 = utc;
-						Frequ_slope_1 = Math1;
-						Moy_slope_1 = Frequ_slope_1;
-						Ch4_slope = Ch4;
-						Moy_Ch4slope_1 = Ch4_slope;
-						N_slope_1 = 1;
+						
+						stat_zero(&stat_math1);
+						stat_zero(&stat_ch4);
 						
 						// next step
 						Measuring_1 += 1;
@@ -686,27 +705,14 @@
 					case N_MEASUREMENT_SLOPE:
 						// slope measurement
 						
-						N_slope_1 = N_slope_1 + 1; 
-						Frequ_slope_1 = Math1;
-						Ch4_slope = Ch4; 
-						Moy_slope_1 = ((N_slope_1-1)*Moy_slope_1 + Frequ_slope_1)/N_slope_1;
-						Moy_Ch4slope_1 = ((N_slope_1-1)*Moy_Ch4slope_1 + Ch4_slope)/N_slope_1;
-						Slope_slope_1 = (Slope_slope_1*(N_slope_1-2) + 6*(Frequ_slope_1-Moy_slope_1)/N_slope_1)/(N_slope_1+1);
-						Slope_Ch4slope_1 = (Slope_Ch4slope_1*(N_slope_1-2) + 6*(Ch4_slope-Moy_Ch4slope_1)/N_slope_1)/(N_slope_1+1);
+						stat_accumulate(&stat_math1, Math1);
+						stat_accumulate(&stat_ch4, Ch4);
 						
 						if ((utc - t1_1) > SlopeTime1) {
-							Slope_1 = Slope_slope_1;
-							Ch4Slope =  Slope_Ch4slope_1;
+							Slope_1 = stat_math1.slope;
+							Ch4Slope =  stat_ch4.slope;
 							SetCtrlVal(CalcN1Panel, CALCN1_SLOPE, Slope_1);
 							
-							N_slope_1 = 0;
-							Frequ_slope_1 = 0.0;
-							Moy_slope_1 = 0.0;
-							Slope_slope_1 = 0.0;
-							Ch4_slope = 0.0;
-							Moy_Ch4slope_1 = 0.0;
-							Slope_Ch4slope_1 = 0.0;
-							
 							// frep positive step
 							DDS4xAD9912_FrequencyRampe(&DDS4xAD9912,1, FrequDDS1,(FrequDDS1+DeltakHz_1*1000), Step1/Ndiv);
 							SetCtrlVal(MainPanel, PANEL_DDS1, (FrequDDS1+DeltakHz_1*1000));
@@ -834,11 +840,10 @@
 						GetCtrlVal(MainPanel, PANEL_DDS2, &FrequencyDDSBesInit);
 						GetCtrlVal(MainPanel, PANEL_DDS3, &FrequencyDDS3Init);
 						t1_2 = utc;
-						Frequ_slope_2 = Math1;
-						Beat_slope_2 = Ch2;
-						Moy_slope_2 = Frequ_slope_2;
-						Moy_Beatslope_2 = Beat_slope_2;
-						N_slope_2 = 1;
+						
+						stat_zero(&stat_math1);
+						stat_zero(&stat_ch2);
+						
 						Nu1 = N1 * (250000000 + Math1);
 						
 						// next step
@@ -848,26 +853,13 @@
 					case N_MEASUREMENT_SLOPE:
 						// slope measurement
 						
-						N_slope_2 = N_slope_2 + 1; 
-						Frequ_slope_2 = Math1;
-						Beat_slope_2 = Ch2;  
-						Moy_slope_2 = ((N_slope_2-1)*Moy_slope_2 + Frequ_slope_2)/N_slope_2;
-						Moy_Beatslope_2 = ((N_slope_2-1)*Moy_Beatslope_2 + Beat_slope_2)/N_slope_2;
-						Slope_slope_2 = (Slope_slope_2*(N_slope_2-2) + 6*(Frequ_slope_2-Moy_slope_2)/N_slope_2)/(N_slope_2+1);
-						Slope_Beatslope_2 = (Slope_Beatslope_2*(N_slope_2-2) + 6*(Beat_slope_2-Moy_Beatslope_2)/N_slope_2)/(N_slope_2+1);
+						stat_accumulate(&stat_math1, Math1);
+						stat_accumulate(&stat_ch2, Ch2);
 
 						if ((utc - t1_2) > SlopeTime2) {
-							Slope_2 = Slope_slope_2;
-							Beatslope_2 = Slope_Beatslope_2;
+							Slope_2 = stat_math1.slope;
+							Beatslope_2 = stat_ch2.slope;
 							SetCtrlVal(CalcN2Panel, CALCN2_SLOPE, Beatslope_2);
-
-							N_slope_2 = 0;
-							Frequ_slope_2 = 0.0;
-							Moy_slope_2 = 0.0;
-							Slope_slope_2 = 0.0;
-							Moy_Beatslope_2 = 0.0;
-							Slope_Beatslope_2 = 0.0;
-							Beat_slope_2 = 0.0;
 							
 							// frep positive step
 							double fDDS1 = FrequDDS1 + DeltakHz_2 * 1000;
@@ -1017,7 +1009,8 @@
 						settling = 3;
 						
 						t1_3 = utc;
-						N_slope_3 = 0;
+						stat_zero(&stat_ch2);
+						
 						// record current DDS3 frequency
 						GetCtrlVal(MainPanel, PANEL_DDS3, &FrequencyDDS3Init);
 						
@@ -1033,20 +1026,14 @@
 							break;
 						}
 						
-						N_slope_3++;
-						Frequ_slope_3 = Ch2;
-						Moy_slope_3 = ((N_slope_3-1)*Moy_slope_3 + Frequ_slope_3)/N_slope_3;
-						Slope_slope_3 = (Slope_slope_3*(N_slope_3-2) + 6*(Frequ_slope_3-Moy_slope_3)/N_slope_3)/(N_slope_3+1);
+						stat_accumulate(&stat_ch2, Ch2);
 						
 						if (utc - t1_3 > SlopeTime3) {
 							// slope measurement
-							Slope_3 = Slope_slope_3;
+							Slope_3 = stat_ch2.slope;
 							
 							t2_3 = utc;
-							N_slope_3 = 0;
-							Frequ_slope_3 = 0.0;
-							Moy_slope_3 = 0.0;
-							Slope_slope_3 = 0.0;
+							stat_zero(&stat_ch2);
 							
 							// frep positive step
 							SetCtrlVal(MainPanel, PANEL_DDS4, FrequDDS4 + DeltakHz_3 * 1000);
@@ -1193,20 +1180,20 @@
 				}
 				
 				// select reference
-				double currentFreq = 0.0;
+				double current = 0.0;
 				switch (slopeReference) {
 					case SLOPE_REFERENCE_MICROWAVE:
-						currentFreq = Math2;
+						current = Math2;
 						break;
 					case SLOPE_REFERENCE_HG_CAVITY:
-						currentFreq = Ch2 * 1062.5 / 1542.2;
+						current = Ch2 * 1062.5 / 1542.2;
 						break;
 				}
 				
 				// stop slope cancelling if the comb is not locked
 				if (SlopeMeasuring & StopSlopeCancellingOnUnlocked
-					& (previousFreq != 0.0)
-					& (fabs(currentFreq - previousFreq) > limitotakoff)) {
+					& (freq.previous != 0.0)
+					& (fabs(current - freq.previous) > limitotakoff)) {
 						
 					if (! KeepSlope) {
 						appliedSlope = 0.0;
@@ -1216,8 +1203,9 @@
 					if (! KeepFrequ) {
 						DDSFox_Set(&DDS1xAD9956, DEDRIFT_DDS_FREQUENCY, appliedSlope);
 					}
-					measuredSlope = 0.0;
-					SetCtrlVal(MainPanel, PANEL_SLOPE_MEASURED, measuredSlope);
+					
+					stat_zero(&freq);
+					SetCtrlVal(MainPanel, PANEL_SLOPE_MEASURED, freq.slope);
 					SlopeMeasuring = FALSE;
 					SetCtrlVal(MainPanel, PANEL_MEASURE_SLOPE, 0);
 				}
@@ -1226,30 +1214,25 @@
 				if (SlopeMeasuring)
 				{
 					// update slope measurement
-					Nsamples = Nsamples + 1;
-					if (Nsamples > 1)
-						measuredSlope = (measuredSlope * (Nsamples - 2) + 6 * (currentFreq - meanFreq) / Nsamples) / (Nsamples + 1);
-					else
-						measuredSlope = 0.0;
-					meanFreq = ((Nsamples - 1) * meanFreq + currentFreq) / Nsamples;
-					previousFreq = currentFreq;
+					stat_accumulate(&freq, current);
+					
 					// update indicator
-					SetCtrlVal(MainPanel, PANEL_SLOPE_MEASURED, measuredSlope);
+					SetCtrlVal(MainPanel, PANEL_SLOPE_MEASURED, freq.slope);
 					
 					// update applied slope
 					if ((utc - SlopeMeasuringTimeBegin) > TimetoSlope) {
 						
 						if (invertSlopeSign)
-							appliedSlope = appliedSlope - measuredSlope;
+							appliedSlope = appliedSlope - freq.slope;
 						else
-							appliedSlope = appliedSlope + measuredSlope;
+							appliedSlope = appliedSlope + freq.slope;
 						
 						if (FrequCorrec) {
 							// proportional correction
 							
 							Nratio += 1;
 							if (Nratio >= 1) {
-								MoyMath2 = MoyMath2 + meanFreq;
+								MoyMath2 = MoyMath2 + freq.mean;
 							}
 							if (Nratio == 1 && CenterFrequencyCh2ToDetermine == TRUE) {
 								CenterFrequencyCh2 = MoyMath2;
@@ -1266,9 +1249,7 @@
 						SetCtrlVal(MainPanel, PANEL_SLOPE_APPLIED, appliedSlope);
 						DDSFox_SetSweepRate(&DDS1xAD9956, appliedSlope);
 						
-						Nsamples = 0;
-						meanFreq = 0.0;
-						measuredSlope = 0.0;
+						stat_zero(&freq);
 						SlopeMeasuringTimeBegin = utc;
 					}
 				}
@@ -1942,18 +1923,6 @@
 					HidePanel(CalcN1Panel);
 				
 				Measuring_1 = FALSE;
-				Frepminus_1=0.0;
-				Frepplus_1=0.0;
-			  	t1_1=0.0;
-			 	t2_1=0.0;
-			 	t3_1=0.0; 
-				N_slope_1=0;
-				Frequ_slope_1=0.0;
-				Moy_slope_1=0.0;
-				Slope_slope_1 =0.0;
-				Ch4_slope=0.0;
-				Moy_Ch4slope_1=0.0;
-				Slope_Ch4slope_1=0.0;
 					
 				SetCtrlVal(MainPanel, PANEL_DDS1, FrequDDS1) ;
 				DDS4xAD9912_SetFrequency(&DDS4xAD9912, 1, FrequDDS1);
@@ -1965,22 +1934,7 @@
 				if (PanelIsVisible)
 					HidePanel(CalcN2Panel);
 				
-				Measuring_2=FALSE;
-				Frepminus_2=0.0;
-				Delta10K_Minus=0.0;
-				Frepplus_2=0.0;
-				Delta10K_Plus=0.0; 
-				DeltaDDS3=0.0; 
-			 	t1_2=0.0;
-			  	t2_2=0.0;
-			  	t3_2=0.0; 
-				N_slope_2=0;
-				Frequ_slope_2=0.0;
-				Moy_slope_2=0.0;
-				Beat_slope_2=0.0;
-				Moy_Beatslope_2=0.0;
-				Slope_Beatslope_2 =0.0;
-				Slope_slope_2 =0.0;
+				Measuring_2 = FALSE;
 				
 				SetCtrlVal(MainPanel, PANEL_DDS1, FrequDDS1) ;
 				DDS4xAD9912_SetFrequency(&DDS4xAD9912, 1, FrequDDS1);
@@ -1995,15 +1949,6 @@
 					HidePanel(CalcN3Panel);
 				
 				Measuring_3 = FALSE;
-				Frepminus_3 = 0.0;
-				Frepplus_3 = 0.0;
-			   	t1_3 = 0.0;
-			  	t2_3 = 0.0;
-			   	t3_3 = 0.0; 
-				N_slope_3 = 0;
-				Frequ_slope_3 = 0.0;
-				Moy_slope_3 = 0.0;
-				Slope_slope_3 = 0.0;
 			}
 		break;
 	}
@@ -2196,10 +2141,7 @@
 				
 				SlopeMeasuringTimeBegin = utc;
 				
-				Nsamples = 0;
-				previousFreq = 0.0;
-				meanFreq = 0.0;
-				measuredSlope = 0.0;
+				stat_zero(&freq);
 				
 				Nratio = -1;
 				MoyMath2 = 0.0; 
@@ -2215,8 +2157,8 @@
 				if (! KeepFrequ) {
 					DDSFox_Set(&DDS1xAD9956, DEDRIFT_DDS_FREQUENCY, appliedSlope);
 				}
-				measuredSlope = 0.0;
-				SetCtrlVal(panel, PANEL_SLOPE_MEASURED, measuredSlope);
+				stat_zero(&freq);
+				SetCtrlVal(panel, PANEL_SLOPE_MEASURED, freq.slope);
 			}
 			break;
 	}