#pragma rtGlobals=1		// Use modern global access method.
//Allan variance function for 1-D wave
//9/13/12 maj

//***********************************************//
// from Wikipedia "Allan variance"--									//
// "Allan variance is defined as one half of the time average of the 		//
// squares of the differences between successive readings of the 		//
// signal deviation sampled over the sampling period. The Allan	// 
// variance depends on the time period used between samples: 			//
// therefore it is a function of the sample period, commonly 			//
// denoted as tau, likewise the distribution being measured, and is		// 
// displayed as a graph rather than a single number. "					//
//***********************************************//

Menu "Allan variance"
	"Plot Allan Variance of a wave plotted on the top graph",AllanVariance()
	"Add lines for white noise", AllanWhite()
End

	
Function AllanVariance()

	String wName
	Variable nMax = 50
	Prompt wName, "Plot Allan variance of",popup, WaveList("*", ";","WIN:")
	Prompt nMax, "Maximum number of variances:"
	DoPrompt "Allan variance",wName, nMax
	if(V_flag)	//user canceled
		Abort "User canceled Allan variance."
	endif
	WAVE w = $wName
	Variable wd
	wd = WaveDims(w)
	If(wd !=1)
		Abort "Allan variance for 1-D waves only."
	endif
	String aaxis,avar,alertStr
	aaxis = wName + "_avx"  // x-axis for output wave
	avar = wName + "_avar" // output wave
	If((Exists(aaxis) == 1) %| (Exists(avar)==1))
	alertStr = "Wave(s) for Allan variance of \'"+wName
	alertStr += "\' already exist, overwrite them?"
		DoAlert 1, alertStr
		If(V_flag == 2)
			Abort "User aborted AllanVariance."
		endif
	endif	
	Variable npt,dX,n2,n3,n4
	npt = numpnts(w)
	dX = deltax(w)
	n2 = 2*floor(sqrt(npt))  // later limit n2 to nMax points
	Make/O/N=(n2) $aaxis
	WAVE aaxisW = $aaxis
	aaxisW = (p+1)*dX
	aaxisW[(n2/2+1),(n2-1)] = dX*round(npt/(n2-p+1))
	If(n2 > nMax)  // replace middle point by  exp spaced pts
		n3 = floor(nMax/3)
		n4 = n2-2*n3
		DeletePoints n3,n4,aaxisW
		InsertPoints n3,n3,aaxisW
		Variable j=0
		Variable h = (aaxisW[(2*n3)]/aaxisW[(n3-1)])^(1/(n3+1))
		Variable hh = h
		Do
			aaxisW[j+n3] = aaxisW[(n3-1)]*hh
			hh *=h
			j += 1
		While(j <= n3)
	endif
	Make/O/N=(numpnts(aaxisW)) $avar
	WAVE avarW = $avar
	avarW = FAllanVar(w, aaxisW[p])
	Display/K=1  avarW vs aaxisW
	AllanVarianceStyle()
	
	return 0
End


Function FAllanVar(inWave, xInt)
	WAVE inWave
	Variable xInt   // x interval for averaging
	
	Variable npt = numpnts(inWave)
	Variable dX = deltax(inWave)   //  x-spacing for inWave
	Variable numInt, ptsInt
	numInt = ceil(npt*dX/xInt)  // number of intervals
	ptsInt = ceil(npt/numInt)   // number of pts per interval
	Make/FREE/N=(numInt) tempAllan
	tempAllan = 0

	Variable i = 0,p1,p2
	Do
		p1 = pnt2x(inWave,i*ptsInt)
		p2 = pnt2x(inWave,(i+1)*ptsInt-1)
		tempAllan[i] = mean(inWave,p1,p2)
		i += 1
	While(i < numInt)
	tempAllan = tempAllan[p+1] - tempAllan[p]
	DeletePoints (numInt-1),1,tempAllan
	tempAllan = tempAllan*tempAllan/2
	If(numInt >2)
		WaveStats/Q tempAllan
		return V_avg
	else
		return tempAllan[0]
	endif
end
	
	
// Draw lines for white noise
Function AllanWhite()

	String wName
	Prompt wName, "Add white noise lines for",popup, WaveList("*", ";","WIN:")
	DoPrompt "Add white noise lines", wName
	if(V_flag)	//user canceled
		Abort "User canceled white noise lines."
	endif
	String line1,line2,line3, aaxis,avar,alertStr,wn
	Variable uchar= 0,start = 0, dX
	Do   // find last "_" in wave name
		start = uchar+1
		uchar = strsearch(wName, "_",start)
	While(uchar != -1)
	If(start==1)
		alertStr = "Cannot find Allan variance wave(s). "
		Abort alertStr
	endif
	If(cmpstr(wName[(start),(start+1)],"av") !=0)
		alertStr = "Cannot find Allan variance wave(s). "
		Abort alertStr
	endif
	wn = wName[0,(start-2)]
	WAVE w = $wn
	If(WaveExists(w) != 1)
		alertStr = "Cannot find \'"+wn+"\'"
		Abort alertStr
	endif
	dX = deltax(w)
	line1 = wn+"_l1"
	line2 = wn+"_l2"
	line3 = wn+"_l3"
	aaxis = wn+"_avx"
	avar = wn+"_avar"
	If((Exists(line1) == 1) %| (Exists(line2)==1) %| (Exists(line3)==1))
		alertStr = "White noise wave(s) for  \'"+wName
		alertStr += "\' already exist, overwrite them?"
		DoAlert 1, alertStr
		If(V_flag == 2)
			Abort "User aborted AllanWhite"
		endif
	endif
	If((Exists(aaxis) != 1) %| Exists(avar) != 1)
		alertStr = "Cannot find \'"+aaxis+"\' (or \'"+avar+"\')"
		Abort alertStr
	endif
	Make/O/N=(numpnts($aaxis)) $line1,$line2,$line3
	WAVE lineW1 = $line1
	WAVE lineW2 = $line2
	WAVE lineW3 = $line3
	WAVE aaxisW = $aaxis
	WAVE avarW = $avar
	lineW2 = avarW[0]/aaxisW[p]*dX
	Variable npt = numpnts($wn)
	lineW1 =avarW[0]* dX*sqrt(2/npt/aaxisW[p]) + lineW2[p]
	lineW3 = -avarW[0]*dX*sqrt(2/npt/aaxisW[p]) + lineW2[p]
	
//	GetAxis/Q left	//gets left axis limits
	AppendToGraph lineW1,lineW2,lineW3 vs aaxisW
	ModifyGraph lstyle($line1)=7,lstyle($line2)=1,lstyle($line3)=7
	ModifyGraph lsize($line1)=1,lsize($line2)=1,lsize($line3)=1
//	SetAxis left,V_min,V_max

	return 0
End 


Function AllanVarianceStyle()
	ModifyGraph/Z grid=2
	ModifyGraph/Z log=1
	ModifyGraph/Z mirror=1
	Label/Z left "Allan variance"
	Label/Z bottom "Averaging time"
End