Allan variance
CJanssen
Sat, 01/28/2012 - 02:43 am
Sat, 01/28/2012 - 02:43 am
(1) see the article in Wikipedia, http://en.wikipedia.org/wiki/Allan_variance . There is one time-domain definition of the Allan variance expressed in terms of something that looks very much like a finite-difference second derivative. This might be easily calculated with Igor's derivative operation (if you are careful about the end points).
(2) check the Igor mailing list archives (Message-Id: 5133, Date: 2000-03-24, Subject: RE: Allan covariance) for Igor functions offered by Mike Johnson - http://www.wavemetrics.com/search/viewmlid.php?mid=5133 . Perhaps his code is still available.
January 28, 2012 at 03:36 pm - Permalink
Thank you for the comments and the pointers. Meanwhile I have been contacted and provided with the Allan variance code Mike Johnson from Lawrence livermore wrote and alluded to in his email (id 5133) from 2000. The code works flawless and I am happy to mess around with it in order to better understand what it does in detail.
In 2000, Mike Johnson offered to provide his code to whoever is interested. Provided Mike Johnson agrees, the code could be made public on this site.
January 30, 2012 at 07:58 am - Permalink
Thanks.
April 16, 2013 at 05:45 pm - Permalink
April 18, 2013 at 02:54 pm - Permalink
I received a very quick (and positive!) reply from Mike. I'm posting this here rather than as a snippet, because someone should really go over it and modernize the code (at the least, turning macros into functions)
You can test the code with
AllanVariance("test",5000)
AllanWhite("test_avx")
From Mike:
"The attached example shows that the Allan variance of white noise is
equal to the variance of the individual points divided by the length
of the averaging interval. This is the familiar decrease in RMS noise
with the square root of the averaging time."
And here is the code itself:
//4/22/96 maj
Macro AllanVariance(waveName, nMax)
String waveName
Variable nMax=50
Prompt waveName, "Plot Allan variance of",popup, WaveList("*", ";","WIN:")
// (";"+WaveList("*", ";","")) // lists waves in top window first, then
// all waves in current folder
Prompt nMax, "Maximum number of variances:"
Silent 1; PauseUpdate
Variable dt = datetime
Variable wd
wd = WaveDims($waveName)
If(wd !=1)
Abort "Allan variance for 1-D waves only"
endif
String aaxis,avar,alertStr
aaxis = waveName + "_avx" // x-axis for output wave
avar = waveName + "_avar" // output wave
If((Exists(aaxis) == 1) %| (Exists(avar)==1))
alertStr = "Wave(s) for Allan variance of \'"+waveName
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($waveName)
dX = deltax($waveName)
n2 = 2*floor(sqrt(npt)) // later limit n2 to nMax points
Make/O/N=(n2) $aaxis
$aaxis = (p+1)*dX
$aaxis[(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,$aaxis
InsertPoints n3,n3,$aaxis
Variable j=0
Variable h = ($aaxis[(2*n3)]/$aaxis[(n3-1)])^(1/(n3+1))
Variable hh = h
Do
$aaxis[j+n3] = $aaxis[(n3-1)]*hh
hh *=h
j += 1
While(j <= n3)
endif
Make/O/N=(numpnts($aaxis)) $avar
$avar = FAllanVar($waveName,npt,dX,$aaxis[p])
Killwaves/Z tempAllan
Display $avar vs $aaxis
ModifyGraph log=1
dt = datetime-dt
// Printf "Calculation took %5.1f minutes\r", (dt/60)
end
Function FAllanVar(inWave,npt,dX,xInt)
WAVE inWave
Variable npt // number of pts in inWave
Variable dX // x-spacing for inWave
Variable xInt // x interval for averaging
Variable numInt, ptsInt
numInt = ceil(npt*dX/xInt) // number of intervals
ptsInt = ceil(npt/numInt) // number of pts per interval
Make/O/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
Macro AllanWhite(waveName)
String waveName
Prompt waveName, "Add white noise lines for",popup, WaveList("*", ";","WIN:")
Silent 1; PauseUpdate
String line1,line2,line3, aaxis,avar,alertStr
Variable uchar= 0,start = 0
Do // find last "_" in wave name
start = uchar+1
uchar = strsearch(waveName, "_",start)
While(uchar != -1)
If(start==1)
alertStr = "Cannot find Allan variance wave(s) "
Abort alertStr
endif
Print waveName, start
If(cmpstr(waveName[(start),(start+1)],"av") !=0)
alertStr = "Cannot find Allan variance wave(s) "
Abort alertStr
endif
String wn
wn = waveName[0,(start-2)]
If(Exists(wn) != 1)
alertStr = "Cannot find \'"+wn+"\'"
Abort alertStr
endif
Variable dX = deltax($wn)
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 \'"+waveName
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
$line2 = $avar[0]/$aaxis[p]*dX
Variable npt = numpnts($wn)
$line1 =$avar[0]* dX*sqrt(2/npt/$aaxis[p]) + $line2[p]
$line3 = -$avar[0]*dX*sqrt(2/npt/$aaxis[p]) + $line2[p]
GetAxis/Q left
append $line1,$line2,$line3 vs $aaxis
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
end
April 18, 2013 at 04:28 pm - Permalink
Here is the procedure file with functions for calculating the Allan
variance of a 1-D data wave. I have assumed that the data are in the form
of an equally spaced time series expressed as a single scaled Igor wave
with its own X values, although the actual interval between values could
be anything. I also assume that these data are already plotted on the top
graph in Igor.
The new menu item "Allan variance -> Plot Allan Variance of a wave plotted
on the top graph" will bring up a dialog box asking for the name of the
data wave and the approximate number of intervals of different lengths you
want to use in calculating the Allan variance. In principal this number of
intervals can be nearly as large as the number of data values, but this is
usually many more intervals than you want, so I try to give you the number
of intervals you want and vary their lengths in a sensible (approximately
logarithmic) way.
When you click on "Continue" in the dialog you will get a plot of the
Allan variance for your data. This graph is in the form of a log-log
"Ywave vs. Xwave" style of Igor plot. The Y axis (labeled "Allan
variance") has units equal to the square of the units of the Y axis of
your data. The X axis (labeled "Averaging time") has the same units as the
X axis of your data.
There is an additional menu item "Allan variance -> Add lines for white
noise" that shows what the Allan variance for white noise would look like
if that white noise had the the same individual-point variance as your
data. The center line is flanked by dashed lines showing the expected
ranges (I think it is the band of plus and minus one standard deviation,
but I don't exactly remember) of the plot for white noise. These lines are
useful if you want to see how much your data resemble white noise.
Here is a quick test of the functions (after the procedure file is
compiled):
Make/O/N=100000 test = gnoise(3) //white noise with individual-point
variance = 9
Display test //Data are on 1-second time intervals
//Select menu item "Allan variance -> Plot Allan Variance of a wave
plotted on the top graph"
//Use the default arguments for now. A log-log plot appears.
//Select menu item "Allan variance -> Add lines for white noise". Three
lines get added.
//Set the left axis
SetAxis left 1e-05,10
====================================
The newer code is attached
April 18, 2013 at 05:15 pm - Permalink