Allan Variance by Mike Johnson

Average rating
(0 votes)

Mike Johnson posted his Allan Variance code
and I started making minor usability changes to it. I created this project so that others can benefit and collaborate on this code.
Usage: save as an .ipf, load into Igor, plot the data you will analyze, continue through the AllanVariance menu.

#pragma rtGlobals=1		// Use modern global access method.
//Allan variance function for 1-D wave
//9/13/12 maj
//2017-09-08 jcc, a few usability tweaks
// originally posted at  <a href="<br />
//" title="<br />
//" rel="nofollow"><br />
//</a> original file <a href="<br />
//" title="<br />
//" rel="nofollow"><br />
//</a> original author, Mike Johnson <a href="</p>
<p>//***********************************************" title="</p>
<p>//***********************************************" rel="nofollow"></p>
<p>//*******************************...</a>						//
// 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()
// creates two new waves containing allan variance and averaging window
Function AllanVariance([wName,nMax,p0,p1])
	String wName	// wave to use
	Variable nMax	// maximum npnts to average for Allan variance
	variable p0,p1	// optionally use only wName[p0,p1] for the calculation[jcc]
//	variable getDev // optionally return the deviation instead of variance (sqrt(var)) [future feature]
	If (ParamIsDefault(nMax) || ParamIsDefault(wName))
		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."
	WAVE w = $wName
	Variable wd
	wd = WaveDims(w)
	If(wd !=1)
		Abort "Allan variance is available for 1-D waves only."
	if (!ParamIsDefault(p0) || !ParamIsDefault(p1))
		make /free/d/n=(p1-p0+1) wtrimmed
		wtrimmed= w[p+p0]
		wave w= wtrimmed
	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."
	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
			aaxisW[j+n3] = aaxisW[(n3-1)]*hh
			hh *=h
			j += 1
		While(j <= n3)
	Make/O/N=(numpnts(aaxisW)) $avar
	WAVE avarW = $avar
	avarW = FAllanVar(w, aaxisW[p])
	Display/K=1  avarW vs aaxisW
	// get also white noise
//	// convert to deviation? [future feature]
//	if (getDev)
//		avarW= sqrt(avarW)
//		Label/Z left "Allan deviation" // AllanVarianceStyle() had labelled it before
//	endif
	return 0
// this function actually performs the Allan variance analysis
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
		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
		return tempAllan[0]
// Draw lines for white noise
Function AllanWhite([wName])
	String wName	// wave to use
	if (ParamisDefault(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."
	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)
		alertStr = "Cannot find Allan variance wave(s). "
		Abort alertStr
	If(cmpstr(wName[(start),(start+1)],"av") !=0)
		alertStr = "Cannot find Allan variance wave(s). "
		Abort alertStr
	wn = wName[0,(start-2)]
	WAVE w = $wn
	If(WaveExists(w) != 1)
		alertStr = "Cannot find \'"+wn+"\'"
		Abort alertStr
	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"
	If((Exists(aaxis) != 1) %| Exists(avar) != 1)
		alertStr = "Cannot find \'"+aaxis+"\' (or \'"+avar+"\')"
		Abort alertStr
	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
Function AllanVarianceStyle()
	ModifyGraph/Z grid=2
	ModifyGraph/Z log=1
	ModifyGraph/Z mirror=1
	Label/Z left "Allan variance"
	Label/Z bottom "Averaging time"

Back to top