A new release of AllanTools is now available on PyPi.

#### Tagallantools

The AllanTools functions mtotdev() and htotdev() are now almost 10-times faster, after an update to the code that more efficiently calculates moving averages.

The old code used **numpy.mean()** for each iteration of a loop:

1 2 3 4 |
for j in range(0, 6*m): # summation of the 6m terms. xmean1 = np.mean(xstar[j : j+m]) xmean2 = np.mean(xstar[j+m : j+2*m]) xmean3 = np.mean(xstar[j+2*m : j+3*m]) |

However this can be computed much faster by noticing that the new mean differs from the old (already computed!) mean by just two points, one at the start is dropped, and a new one at the end is added:

1 2 3 4 5 6 7 8 9 10 11 |
for j in range(0, 6*m): # summation of the 6m terms. if j == 0: # intialize the sum xmean1 = np.sum( xstar[0:m] ) xmean2 = np.sum( xstar[m:2*m] ) xmean3 = np.sum( xstar[2*m:3*m] ) else: # j>=1, subtract old point, add new point xmean1 = xmean1 - xstar[j-1] + xstar[j+m-1] # xmean2 = xmean2 - xstar[m+j-1] + xstar[j+2*m-1] # xmean3 = xmean3 - xstar[2*m+j-1] + xstar[j+3*m-1] # |

Recently merged plotting-functions for allantools with a look and feel of Stable32 - thanks to yxie.

I should probably run python3 to get special characters like tau and sigma to work..

Inspired by discussion on time-nuts, here's a revised noise-colour graph. There are a few updates: The PSDs (both phase and frequency) now cross at 1 Hz (with the relation between phase-PSD and frequency-PSD explicitly stated), and the ADEV/MDEV theoretical lines now include the formula for the pre-factor (the old graph only had 'proportional to' here).

Source: example_noise_slopes2.py

It is left as an exercise to the reader to properly explain the ADEV pre-factor for flicker-phase noise :).

By popular demand, AllanTools 2018.03 is now released under LGPL license.

Work with integrating noise-identification with downstream confidence-interval estimation and bias correction continues as time permits.

For AllanTools statistics both bias-correction and confidence interval calculation requires identifying the dominant power law noise in the input time-series.

The usual noise-types studied have phase PSD noise-slopes "**b**" ranging from 0 to -4 (or even -5 or -6), corresponding to frequency PSD noise-slopes "**a**" ranging from +2 to -2 (where **a=b+2**). These 'colors of noise' can be visualized like this:

I've implemented three noise-identification algorithms based on Stable32-documentation and other papers: B1, R(n), and Lag-1 autocorrelation.

**B1** (Howe 2000, Barnes1969) is defined as the ratio of the standard N-sample variance to the (2-sample) Allan variance. From the definitions one can derive an expected B1 ratio of (the length of the time-series is **N**)

where **mu** is the tau-exponent of Allan variance for the noise-slope defined by **b** (or **a**). Since **mu** is the same (-2) for both b=0 and b=-1 (red and green data) we can't use B1 to resolve between these noise-types. B1 looks like a good noise-identifier for b=[-2, -3, -4] where it resolves very well between the noise types at short tau, and slightly worse at longer tau.

**R(n)** can be used to resolve between b=0 and b=-1. It is defined as the ratio MVAR/AVAR, and resolves between noise types because MVAR and AVAR have different tau-slopes **mu**. For b=0 we expect mu(AVAR, b=0) = -2 while mu(MVAR, b=0)=-3 so we get mu(R(n), b=0)=-1 (red data points/line). For b=-1 (green) the usual tables predict the same mu for MVAR and AVAR, but there's a weak log(tau) dependence in the prefactor (see e.g. Dawkins2007, or IEEE1139). For the other noise-types b=[-2,-3,-4] we can't use R(n) because the predicted ratio is one for all these noise types. In contrast to B1 the noise identification using R(n) works best at large tau (and not at all at tau=tau0 or AF=1).

The **lag-1 autocorrelation** method (Riley, Riley & Greenhall 2004) is the newest, and uses the predicted lag-1 autocorrelation for (WPM b=0, FPM b=-1, WFM b=-2) to identify noise. For other noise types we differentiate the time-series, which adds +2 to the noise slope, until we recognize the noise type.

Here are three figures for ACF, B1, and R(n) noise identification where a simulated time series with known power law noise is first generated using the Kasdin&Walter algorithm, and then we try to identify the noise slope.

For Lag-1 ACF when we decimate the phase time-series for AF>1 there seems to be a bias to the predicted **a** (alpha) for b=-1, b=-3, b=-4 which I haven't seen described in the papers or understand that well. Perhaps an aliasing effect(??).

I've added confidence interval estimation to allantools, based on a 2004 paper by Greenhall & Riley: "Uncertainty of stability variances based on finite differences"

So far not much is automated so you have to run everything manually. After a normal call to **allantools.adev()** to calculate the ADEV we loop through each (tau, adev) pair and first call **allantools.edf_greenhall()** to get the number of equivalent-degrees-of-freedom (EDF), and then evaluate a confidence interval with **allantools.confidence_interval()**. A knowledge or estimate of the noise type "alpha" is required for **edf_greenhall()** - here we just assume alpha=0.

This example is on github at: https://github.com/aewallin/allantools/blob/master/examples/ci_demo.py

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 |
import numpy import matplotlib.pyplot as plt import allantools as at # this demonstrates how to calculate confidence intervals for ADEV # using the algorithms from Greenhall2004 data_file = '../tests/phasedat/PHASE.DAT' def read_datafile(filename): p = [] with open(filename) as f: for line in f: if not line.startswith("#"): # skip comments p.append(float(line)) return p # read input data from file phase = read_datafile(data_file) # normal ADEV computation, giving naive 1/sqrt(N) errors (taus,devs,errs,ns) = at.adev(phase, taus='octave') # Confidence-intervals for each (tau,adev) pair separately. cis=[] for (t,dev) in zip(taus,devs): # Greenhalls EDF (Equivalent Degrees of Freedom) # alpha +2,...,-4 noise type, either estimated or known # d 1 first-difference variance, 2 allan variance, 3 hadamard variance # we require: alpha+2*d >1 (is this ever false?) # m tau/tau0 averaging factor # N number of phase observations edf = at.edf_greenhall( alpha=0, d=2, m=t, N=len(phase), overlapping = False, modified=False ) # with the known EDF we get CIs # for 1-sigma confidence we set # ci = scipy.special.erf(1/math.sqrt(2)) = 0.68268949213708585 (lo,hi) = at.confidence_intervals( dev=dev, ci=0.68268949213708585, edf=edf ) cis.append( (lo,hi) ) # now we are ready to print and plot the results print "Tau\tmin Dev\t\tDev\t\tMax Dev" for (tau,dev,ci) in zip(taus,devs,cis): print "%d\t%f\t%f\t%f" % (tau, ci[0], dev, ci[1] ) """ output is Tau min Dev Dev Max Dev 1 0.285114 0.292232 0.299910 2 0.197831 0.205102 0.213237 4 0.141970 0.149427 0.158198 8 0.102541 0.110135 0.119711 16 0.056510 0.062381 0.070569 32 0.049153 0.056233 0.067632 64 0.027109 0.032550 0.043536 128 0.026481 0.033855 0.055737 256 0.007838 0.010799 0.031075 """ plt.figure(figsize=(12,8)) plt.gca().set_yscale('log') plt.gca().set_xscale('log') err_lo = [ d-ci[0] for (d,ci) in zip(devs,cis)] err_hi = [ ci[1]-d for (d,ci) in zip(devs,cis)] plt.errorbar(taus, devs, yerr=[ err_lo, err_hi ] ,fmt='o') plt.grid() plt.xlabel('Tau (s)') plt.ylabel('ADEV') plt.title('AllanTools 2016.11 - now with Confidence Intervals!') # just to check plot the intervals as dots also plt.plot(taus, [ci[0] for ci in cis],'r.') plt.plot(taus, [ci[1] for ci in cis],'g.') plt.show() |

Allantools has a dependence on scipy because it uses scipy.stats.chi2.ppf() (inverse of the chi-squared cumulative distribution function) in the code for confidence intervals.

When testing scipy takes quite a long time to install and the whole package seems a bit overkill for just this one function.

So I tried implementing it using just numpy. It sort of works, but not for corner cases where p is very close to 1.0 or k is a big number. I get for example:

1 2 |
top = pow(x, s)*math.exp(-x)*pow(x,k) OverflowError: (34, 'Numerical result out of range') |

Maybe there's something trick when numerically evaluating lower_gamma() using the series expansion that I am missing...

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 |
import scipy.stats import math import numpy def bisection(function, k, p, a, b, tol): # http://code.activestate.com/recipes/578417-bisection-method-in-python/ assert (function(a,k)-p)*(function(b,k)-p) < 0 # a, b must bracket root c = (a+b)/2.0 while (b-a)/2.0 > tol: if (function(c, k)-p) == 0: return c elif (function(a,k)-p)*(function(c,k)-p) < 0: b = c else : a = c c = (a+b)/2.0 return c def lower_gamma(s,x): # lower incomplete gamma function # https://en.wikipedia.org/wiki/Incomplete_gamma_function#Evaluation_formulae g = 0 last_g = 1.0 done = False tol = 1.0e-6 k=0 while not done: top = pow(x, s)*math.exp(-x)*pow(x,k) bot = numpy.prod( [float(s+j) for j in range(k+1) ] ) dg = float(top)/float(bot) if dg == float("Inf"): break g += dg k += 1 if k>100: # get at least 100 terms in the sum if g==0: break delta = abs(dg/g) if delta == float("Inf"): break if delta < tol: done = True last_g = g return g def chi2_cdf(x, k): # chi-squared cumulative density function # cdf(x; k) = lower_gamma(k/2, x/2) / gamma(k/2) return lower_gamma(k/2.0, x/2.0) / math.gamma(k/2.0) def chi2_ppf(p, k): # chi-squared Percent point function (inverse of cdf percentiles). # look for x such that # p = chi2_cdf( x=chi2_ppf(p, k), k) tol = 1e-8 lolim = 0 hilim = k while (chi2_cdf(lolim,k)-p)*(chi2_cdf(hilim,k)-p) > 0: hilim *= 1.5 return bisection( chi2_cdf, k, p, lolim, hilim, tol) print "scipy cdf: ",scipy.stats.chi2.cdf(55, 33) print "own cdf: ",chi2_cdf(55, 33) print "scipy ppf ", scipy.stats.chi2.ppf(0.4, 33) print " own ppf ", chi2_ppf(0.4, 33) # test that we really found the inverse print scipy.stats.chi2.cdf(scipy.stats.chi2.ppf(0.4, 33), 33) print chi2_cdf( chi2_ppf(0.4, 33), 33 ) # try to check the scipy function against our own function # for some random input of (p, k) for n in range(100): k = numpy.random.randint(20, 200) p = numpy.random.random() print k, p, a=scipy.stats.chi2.ppf(p, k) b=chi2_ppf(p, k) ok = numpy.isclose(a, b) if ok: print ok else: print ok, a, b assert ok |

New shiny project badges for allantools!

- package on pypi: https://pypi.python.org/pypi/AllanTools/2016.4
- continuous testing on Travis-CI: https://travis-ci.org/aewallin/allantools
- documentation on readthedocs: http://allantools.readthedocs.io/en/latest/?badge=latest
- test coverage report on coveralls: https://coveralls.io/github/aewallin/allantools?branch=master

Update: now with the colours matching in all graphs:

Time-series generated with colorednoise (following Kasdin&Walter), power-spectral-densities and Allan deviations computed with allantools, and compared to theoretical predictions in IEEE1139-2008.

The PSD lines and MDEV lines seem spot-on, but are the ADEV lines systematically a bit low?

Code here: example_noise_slopes.py

© 2019 anderswallin.net

Theme by Anders Norén — Up ↑

## Recent Comments