allan variance – Numerickly

This is going to be a super quick and dirty writeup until I take the time to make it more presentable, but here we go!

For the time being this is mostly going to just be a link and short description of my allan variance github repo used to create synthetic IMU error models using a first order Gauss-Markov process: https://github.com/EVictorson/allan_variance

Also, I’m going to REALLY simplify the technical discussion here because this seems to be one of those areas where PhDs in control theory like to start waving their… boy parts… around to show who knows more about what, which is unnecessary.

Allan Variance Background

The Allan Variance (AVAR) is historically used as a measure of frequency stability in clocks, oscillators, and amplifiers. It is named after the one and only David W. Allan.

Mathematically the AVAR is is expressed as:
$1) \displaystyle \qquad {\sigma_{y}}^2(\tau)$

While the Allan Deviation is, as you may expect from basic statistics:
$2) \displaystyle \qquad {\sigma_{y}}(\tau)$

We often see the M-sample Allan variance, which is a measure of frequency stability using M samples, time T between measurements, and observation time tau. The M-sample Allan variance is expressed as:

$3) \displaystyle \qquad {\sigma_{y}}^2(M,T,\tau)$

To explain the AVAR as academically as you can in as few words as possible is to say that it is defined as one half of the time average of the squares of the differences between successive readings of the frequency 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. (Editors note: TL;DR) A low Allan variance is characteristic of a clock with good stability over the measured period. [1]

The results of AVAR or ADEV are related to five basic noise terms appropriate for inertial sensor data, namely:

Quantization noise
Angle Random Walk
Bias Instability
Rate Random Walk
Rate Ramp

Now, you may ask, why the hell are we using something that was designed to measure the frequency stability of clocks in characterizing errors of an IMU? The answer is, to put it simply, we typically separate the errors of inertial measurement devices into two groups: deterministic errors like bias and nonlinearity, and stochastic processes like the ones mentioned above. Within the field of system dynamics, everything can be decomposed into frequency domain, and we are therefore concerned about the contribution to errors from different frequency sources.
As wikipedia succinctly puts it: “The Allan variance is intended to estimate stability due to noise processes and not that of systematic errors or imperfections such as frequency drift or temperature effects. The Allan variance and Allan deviation describe frequency stability.” Basically, the easiest way to think about this, if you’ve taken an undergraduate controls course, this will be a piece-meal Bode plot with different error frequencies. I’ll expand on this in a minute, but that’s basically what we are dealing with.

Again, the Allan variance / Allan deviation is nothing more than a Bode plot that analyzes the magnitude of different noise sources with differing frequency power. This will cause different contributing noise terms to appear with differing slopes on the Allan deviation plot; these differing slopes making it easy to see that each noise process differs by an integrator (1/s).

The following error terms are the ones most commonly used, and can be distinguished on the Allan Variance / Deviation plot with their corresponding slopes.

Quantization Noise: slope -1
- Quantization noise is Gaussian and zero mean. It is strictly due to the digital nature of ring laser gyro outputs.
Angle Random Walk: slope = -1/2
- ARW is high frequency noise that can be obvserved as a short-term variation in the sensor output. After integration, it causes random error in angle with distribution that is proportional to the square root of the elapsed time. It is often thought to be from randomized dither and the spontaneous emission of photons for a ring laser gyro.
Bias Instability: slope = 0
- BI has an impact on long-term stability. It is a slowly fluctuating error so it appears in low frequencies as 1/f (pink, or flicker) noise. For a gyroscope, this measures how the bias of the sensor changes over a specified period of time at a constant temperature. For a ring laser gyro this noise originates in the discharge assembly, electronics, and other components susceptible to random flickering.
Rate Random Walk: slope = 1/2
- RRW basically characterizes the time scale over which changes in the bias offset occur and can give an idea of when recalibration of sensors should occur. This random process is of uncertain origin, possibly a limiting case of an exponentially correlated noise with a very long correlation time. Mechanical as well as rate biased laser gyros exhibit this noise term.
Rate Ramp: slope = 1
- Slow monotonic change of output over a long time period. This is more of a deterministic error rather than random noise. Its presence in the data may indicate a very slow monotonic change of the ring laser gyro intensity persisting over a very long period of time.

So, what do one of these bad boys look like you may ask? Well, to steal, er borrow, an image from Freescale’s application note AN5087 “Allan Variance: Noise analysis for gyroscopes” [2] we see something like this:

Allan Deviation Plot, courtesy Freescale AN5087

In the above plot, and with all AVAR / ADEV plots, the horizontal axis is the averaging time, an the vertical axis is either the AVAR or ADEV, ADEV in the case of the plot above.

Empirical Data Acquisition and Testing

For an in-depth tutorial on how best to collect data to perform AVAR, see Section B.4 of [4]. To summarize their recommendations, you can select the sample rate to be at least twice the sample bandwidth. The record length should be at least several times the required performance interval. The recommended approach is to limit the time/frequency domain dynamic range to about 3 orders of magnitude. In practice this is done by dividing the total range into overlapping intervals of geometrically increasing lengths. Thus, the high frequency data is acquired for a short period of time. Lower frequency data is filtered (integrated) and acquired for a longer period; e.g. 100 us data is collected for 0.1s, 0.01s data for 10s, 1s data for 1000s, and 100s data for 10^5 s. The appropriate sampling rates/ record lengths should be chosen to overap about one decade of frequency (time).

If we are considering a gyroscope, the Allan variance can be defined in terms of the output rate, ${\Omega}(t)$ , or the output angle as

$4) \displaystyle \qquad {\theta}(t) = \int^{t}{\Omega}(t')dt'$

The lower integration limit is not present in equation 4 because only angle differences are employed in all of the definitions.

The angle measurements that are fed into equation 4 are made at times $t = k{\tau}_0$ , where k varies from 1 to N (the number of samples), and the sample period is $\tau_0$ . For a discrete set of samples, a cumulative sum can be used to give N values of $\theta$ . In this case the cumulative sum of the gyro output samples at each $k\tau_0$ is taken and each sum obtained is then multiplied by the sample period $\tau_0$ to give N values of $\theta$ . [2]

A general heuristic to use is to set the averaging time to be $\tau = m{\tau}_0$ , where m is the averaging factor. The averaging factor can be chosen arbitrarily as some value according the inequality $m < (N - 1) /2$ . [2]

Cool, now that we’ve got our values of $\theta$ calculated, we compute the Allan variance with the following equation:

$5) \displaystyle \qquad {\sigma}^2(\tau) = {\frac {1}{2\tau^2}} < (\theta_{k+2m} - 2\theta_{k+m} + \theta_k)^2 >$

Equation 5 represents Allan variance as a function of $\tau$ and < > is the ensemble average. Upon expanding the ensemble average in equation 5 we obtain the following estimate of the Allan variance:

$6) \displaystyle \qquad {\sigma}^2(\tau) = {\frac {1}{2\tau^2 (N - 2m)}} \sum_{k=1}^{N-2m} (\theta_{k+2m} - 2\theta_{k+m} + \theta_k)^2$

For those who care, and if you want to swing your boy parts around in public, the Allan variance is related to the two-sided power spectral density, $S_{\Omega}(f)$ by:

$7) \displaystyle \qquad {\sigma}^2(\tau) = 4 \int_0^\infty S_{\Omega}(f) \frac{sin^{4}(\pi f \tau)}{(\pi f \tau)^2} df$

Cool, so hopefully if you do this correctly you’ll have a plot that looks something like this:

Simulation and Validation

Alrighty, so we can get a plot from our time series data easily enough, what next? Well, if we presume to already know the noise coefficients we can simulate the Allan variance plot and compare the two, or we can try to characterize the noise coefficients from the Allan variance plot itself.

For simulation and modelling purposes, typically a first order Gauss-Markov process is used to model these kinds of errors, which is just an exponentially correlated, noise driven stochastic process characterized by the following differential equation:

$8) \displaystyle \qquad \frac{dx}{dt} = \frac{-1}{\tau} x + N(\mu, \sigma)$

Which is ultimately just an exponentially decaying autoregressive term with driven white noise. The discrete time version of this is:

$9) \displaystyle \qquad x(n+1) = x(n) e^{\frac{-1}{\tau}dt} + N(\mu, \sigma) dt$

The first order gauss-markov process may also be defined in terms of its autocorrelation function:

$10) \displaystyle \qquad R_{xx}= \sigma^2 e^{-\beta \tau}$

where $\sigma$ is the entire process variation (NOT JUST THE DRIVEN NOISE VARIATION), $\beta$ is the inverse of the correlation time, $T=\frac{1}{\beta}$ , and $\tau$ is the autocorrelation time lag.

For sufficiently long time series (much longer than the time constant), the autocorrelation plot of the first order Gauss-Markov process can verify the correlation time, which is the time lag at which $R_{xx} = \frac{\sigma^2} {e}$ , and the entire time series variance will appear as the peak at time lag = 0. Note that the driven noise standard deviation cannot be validated via autocorrelation. Also note that matlab has two different autocorrelation functions that will produce slightly different results if your time scales are short (autocorr and xcorr), as well as an autocovariance function, xcov. Autocorr and xcov are normalized by the process mean such that a constant process will have a steady autocorrelation (autocorrelation is just autocovariance scaled by the inverse of the process variance). xcorr, however, does not do this, so a constant process will have a decaying autocorrelation as the time lag increases, which results in a shortening of the summation by 1 element each iteration.

For process lengths on the order of magnitude of the time constant or less the process will be dominated by the integrated white noise term, resulting in a random walk. For time scales much longer than the time constant the exponentially correlated term will begin to be apparent. This can be see on an Allan deviation plot, where for sampling intervals much shorter than the time constant the Gauss-Markov Allan variance reduces to that of a singly integrated white noise process (rate random walk), whose slope is +1/2, and the noise magnitude (standard deviation) may be picked off by finding the intersection of the +1/2 slope line and sampling_interval = 3.
This can also be verified by differentiating this time series to obtain Gaussian noise, whose Allan deviation has a slope of -1/2, and the noise magnitude is interpreted as the intersection of this line with sampling_interval = 1.

For process time scales large enough to see the decay time constant, and Allan deviations with sampling intervals ranging from much less than the time constant to sampling intervals much larger than the time constant, the slope will transition from +1/2 for sampling intervals much shorter than the time constant, to 0 as the sampling interval approaches the time constant at it’s maxima, to -1/2 as the interval becomes much larger than the time constant.

When multiple gauss markov processes are added together it will be nearly impossible to pick apart any of these parameters, but when run individually the driven noise magnitude (qc) and time constant (Tc) can be found on the allan deviation plot as follows:
Tc = argmax_sampling_interval / 1.89, where argmax_sampling_interval is the allan deviation sampling interval that maximizes the allan deviation. The driven noise magnitude may then be found with the following relationship:

$11) \displaystyle \qquad q_c = \frac{\sigma_{max}} {0.437 * \sqrt{T_c}}$

where $\sigma_{max}$ is the Allan deviation maxima.

If the combined output of all of these Gauss-Markov processes are analyzed, the only easy thing to pick off will be the velocity random walk / angle random walk, which may be found by fitting a line to the segment with a slope of -1/2 and finding the intersection of this line with tau = 1.

Characterizing Error Sources

Angle Random Walk

Angle random walk appears on the AVAR plot as the following equation:

$12) \displaystyle \qquad \sigma^2(\tau) = \frac{N^2}{\tau}$

where N is the angle random walk coefficient. In a log-log plot of $\sigma(\tau)$ vs $\tau$ ARW has a slope of -1/2. The numerical value of N can be obtained by reading the -1/2 slope line at $\tau = 1$ .

Bias Instability

Bias instability appears on the AVAR plot as the following cumbersome equation:

$13) \displaystyle \qquad \sigma^2(\tau) = \frac{2B^2}{\pi}\left [ln 2 - \frac{sin^3(x)}{2x^2}\left (sin(x) + 4xcos(x)\right ) + Ci(2x) - Ci(4x) \right ]$

Where x is $\pi f_0 \tau$ ,
B is the bias instability coefficient
$f_0$ is the cutoff frequency, and
Ci is the cosine-integral function.

As you can see it’s probably a bit easier to just pick off the part of the AVAR plot that has a slope of 0.

Rate Random Walk

Rate random walk appears on the AVAR plot with the following equation:

$14) \displaystyle \qquad \sigma^2(\tau) = \frac{K^2 \tau}{3}$

Where K is the rate random walk coefficient.

Because we are essentially dealing with Bode plots, the presence of $\tau$ in the numerator indicates we will be searching for a slope of +1/2. The magnitude of this noise can be read off the +1/2 slope line at $\tau = 3$ .

Rate Ramp

Rate ramp appears on the AVAR plot with the following equation:

$15) \displaystyle \qquad \sigma^2(\tau) = \frac{R^2 \tau^2}{2}$

Where R is the rate ramp coefficient.

With one more power of tau in the numerator we are therefore looking for a slope of +1. The magnitude of R can be obtained from the +1 slope line at $\tau = \sqrt{2}$ .

Quantization Noise

Quantization noise appears on the AVAR plot with the following equation:

$16) \displaystyle \qquad \sigma^2(\tau) = 3\frac{Q^2}{\tau^2}$

Where Q is the quantization noise coefficient.

Quantization noise appears with a slope of -1 on the log-log plot and can be read off the intersection of the line with slope -1 and $\tau = \sqrt{3}$ .

Exponentially Correlated (Markov) Noise

Markov noise appears on the AVAR plot as:

$17) \displaystyle \qquad \sigma^2(\tau) = \frac{(q_c T_c)^2}{\tau}\left [1- \frac{T_c}{2\tau}\left (3-4e^{-\frac{\tau}{T_c}}+e^{-\frac{2\tau}{T_c}}\right )\right ]$

Where qc is the noise amplitude, and Tc is the correlation time as noted above.

It is easier to view the behavior of the Markov noise at various limits of the equation. For cases with $\tau >> T_c$ we obtain:

$18) \displaystyle \qquad \sigma^2(\tau) = \frac{(q_c T_c)^2}{\tau}$

Which is the Allan variance for angle random walk where N = qcTc is the angle random walk coefficient. For $\tau << T_c$ we obtain:

$19) \displaystyle \qquad \sigma^2(\tau) = \frac{q_c ^2}{3}\tau$

Which is the Allan variance for rate random walk.

Sinusoidal Noise

Also, sinusoidal noise exists, but I won’t talk about it because it hasn’t concerned any of my work. Just know it exists.

Now, if we take and superimpose all of these noise sources on one plot we should see a piece-wise representation of the empirically derived plot like that shown in [8]:

References

IA state AERE432 Lecture Notes, retrieved 8/10/2018 http://home.engineering.iastate.edu/~shermanp/AERE432/lectures/Rate%20Gyros/Allan%20variance.pdf
Freescale appnote AN5087, retrieved 8/10/2018 http://cache.freescale.com/files/sensors/doc/app_note/AN5087.pdf
IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Interferometric Fiber Optic Gyros. IEEE std 952-1997 (R2008)
IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Laser Gyros. IEEE std 647-1995
Kirkko-Jaakkola, M., Collin, J., & Takala, J. (2012). Bias Prediction for MEMS Gyroscopes. IEEE Sensors Journal, 12(6), 2157-2163. DOI: 10.1109/JSEN.2012.2185692
Peng-Yu, C. (2015). Experimental Assessment of MEMS INS Stochastic Error Model. Graduate Thesis, The Ohio Sate University, Department of Civil Engineering.
Martin, V. MEMS Gyroscope Performance Comparison Using Allan Variance Method, Brno University of Technology.
ST Design Tip DT0064. Noise analysis and identification in MEMS sensors. Retreived 8/10/2018. https://www.st.com/content/ccc/resource/technical/document/design_tip/group0/bf/92/6b/31/e7/c1/4d/c8/DM00311184/files/DM00311184.pdf/jcr:content/translations/en.DM00311184.pdf

Tag: allan variance

Allan Variance and Its Use in Characterizing Inertial Measurement Unit Errors