# Allan Variance by Mike Johnson

Mike Johnson posted his Allan Variance code http://www.igorexchange.com/node/2757
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  http://www.igorexchange.com/node/2757
// original file http://www.igorexchange.com/files/AllanVariance.ipf
// original author, Mike Johnson http://www.igorexchange.com/user/3700

//***********************************************                       //
// 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. "                  //
//***********************************************                       //

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

// 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."
endif
endif

WAVE w = \$wName
Variable wd
wd = WaveDims(w)
If(wd !=1)
Abort "Allan variance is available for 1-D waves only."
endif

if (!ParamIsDefault(p0) || !ParamIsDefault(p1))
make /free/d/n=(p1-p0+1) wtrimmed
wtrimmed= w[p+p0]
wave w= wtrimmed
endif

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
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()

// get also white noise
AllanWhite(wName=aaxis)

//  // convert to deviation? [future feature]
//  if (getDev)
//      avarW= sqrt(avarW)
//      Label/Z left "Allan deviation" // AllanVarianceStyle() had labelled it before
//  endif

return 0
End

// 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
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([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
endif

if(V_flag)  //user canceled
Abort "User canceled white noise lines."
endif
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). "
endif
If(cmpstr(wName[(start),(start+1)],"av") !=0)
alertStr = "Cannot find Allan variance wave(s). "
endif
wn = wName[0,(start-2)]
WAVE w = \$wn
If(WaveExists(w) != 1)
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
If(V_flag == 2)
Abort "User aborted AllanWhite"
endif
endif
If((Exists(aaxis) != 1) %| Exists(avar) != 1)
alertStr = "Cannot find \'"+aaxis+"\' (or \'"+avar+"\')"
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

Forum

Support

Gallery