Allan Variance by Mike Johnson

#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. "                  //
//***********************************************                       //
 
Menu "Allan variance"
    "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
    
    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()
    
    // 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
    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