Software for Converting Spectral Density to Modified Allan Deviation

[Twigg]

TL;DR

I am looking for a program or programming library that can take the power spectral density (PSD) of a signal as an input and return it's modified allan deviation as an output.

Mathematically, you can convert between a power spectral density (PSD) and the modified allan variance as follows:
$$\sigma_y^2 (\tau) = \int_0^{\infty} \frac{G_\nu(f)}{\nu^2} \times 32 \frac{(\sin(\pi f \tau/2))^4 \times |\sin(\pi f \tau)|^2}{(\pi \tau f)^4} df$$
I was wondering if anyone knew of a piece of software or a programming library that implements this conversion. I'm working on a project and trying to minimize the amount of code I personally need to maintain.

Definitions:
$\nu$: the center frequency of the signal being studied
$\sigma_y(\tau)$: the modified allan deviation of the fractional frequency deviation, $y = \frac{\delta \nu}{\nu}$, as a function of time $\tau$
$G_{\nu}(f)$: the power spectral density of frequency deviations ($\delta \nu$, not $y$)I have already looked into the program stable32 and python's allantools module. As far as I could tell from reading the documentation, stable32 doesn't do a psd-to-allan feature. The allantools module does have a psd2allan function in it's docs, but whenever I pip install the module and import it, that function simply doesn't appear in the module (all the other documented functions do! ). I think it must not be implemented yet or something.

Note to mods: I realize this is a programming question, but I felt it belonged more in the statistics subforum since it's so specialized. Feel free to move it if necessary!

 

[jedishrfu]

One place to check would be the Matlab forums at the math works website. You might get a lead on where else to look or even some implementation that you could port.

 

[jedishrfu]

I found this reference on github

https://github.com/aewallin/allantools/search?q=psd2allan

So they have a testcase which infers the function is embedded in the allantools.py source script.

It may be that your pip pulled an earlier version of the code.

 

[Twigg]

I did notice that, however I already tried installing each of the last 3 versions of allantools (2019.9, 2019.7, and 2018.3), which all had the same effect.

I'm on windows, so I run in the terminal:

pip install allantools==2019.9
python
>>> import allantools
>>> allantools.psd2allan()

and I get as output:

AttributeError: module 'allantools' has no attribute 'psd2allan'

 

[jedishrfu]

Curiously, the test case I referenced has the psd2allan lines commented out so they must be having trouble with the code or it's not in the allantools.py script.

In the allantools.py there's a comment on unreleased code and psd2allen is listed so I guess you'll have to roll your own or find another library that implements it.

[CODE lang="python" title="allantools.py"]“””
Allan deviation tools
=====================
**Author:** Anders Wallin (anders.e.e.wallin "at" gmail.com)
Version history
---------------
**unreleased**
- ITU PRC, PRTC, ePRTC masks for TDEV and MTIE in new file mask.py
- psd2allan() - convert PSD to ADEV/MDEV
- GCODEV
**2019.09** 2019 September
- packaging changes, for conda package
(see https://anaconda.org/conda-forge/allantools)

...

“””[/CODE]

 

[jedishrfu]

I looked in the allantools.py and found they do have code for it so maybe you could use that as a basis and roll a better one.

[CODE lang="python" title="Psd2allan()"]
def psd2allan(S_y, f=1.0, kind='adev', base=2):
""" Convert a given (one-sided) power spectral density S_y(f) to Allan
deviation or modified Allan deviation
For ergodic noise, the Allan variance or modified Allan variance
is related to the power spectral density :math:`S_y` of the fractional
frequency deviation:
.. math::
\\sigma^2_y(\\tau) = 2 \\int_0^\\infty S_y(f)
\\left|sin(\\pi*f*\\tau).^(k+1)./(\\pi*f*tau).^k).^2\\right| df,
where :math:`f` is the Fourier frequency and :math:`\\tau` the averaging
time. The exponent :math:`k` is 1 for the Allan variance and 2 for the
modified Allan variance.
NIST [SP1065]_ eqs (65-66), page 73.
psd2allan() implements the integral by discrete numerical integration via
a sum.
Parameters
----------
S_y: np.array
Single-sided power spectral density (PSD) of fractional frequency
deviation S_y in 1/Hz^2. First element is S_y(f=0).
f: np.array or scalar numeric (float or int)
if np.array: Spectral frequency vector in Hz
if numeric scalar: Spectral frequency step in Hz
default: Spectral frequency step 1 Hz
kind: {'adev', 'mdev'}
Which kind of Allan deviation to compute. Defaults to 'adev'
base: float
Base for logarithmic spacing of tau values. E.g. base= 10: decade,
base= 2: octave, base <= 1: all
Returns
-------
(taus_used, ad): tuple
Tuple of 2 values
taus_used: np.array
tau values for which ad computed
ad: np.array
Computed Allan deviation of requested kind for each tau value
"""

"""
# determine taus from df
# first oversample S_y by a factor of 10 in order to avoid numerical
# problem at tau > 1/2df
if isinstance(S_y, np.ndarray):
if isinstance(f, np.ndarray): # f is frequency vector
df = f[1]-f[0]
elif np.isscalar(f):
# assume that f is the frequency step, not frequency vector
df = f
else:
raise ValueError(np.ndarray, float, int)
else:
raise ValueError(np.ndarray) # raise error
oversamplingfactor = 4
df0 = oversamplingfactor * df
f0 = np.arange(S_y.size * df0, step=df0)
f = np.arange(df, (S_y.size - 1) * df0 + df, df)
S_y = interpolate.interp1d(f0, S_y, kind='cubic')(f)
f = f / oversamplingfactor

tau0 = 1/np.max(f) # minimum tau derived from the given frequency vector
n = 1/df/tau0/2
if base > 1:
m = np.unique(np.round(np.append(base**np.arange(
np.floor(np.log(n)/np.log(base))), n)))
else:
m = np.arange(1, n)
taus_used = m*tau0

# TODO: In principle, this approach can be extended to the other kinds of
# Allan deviations, we just need to determine the respective transfer
# function in the frequency domain.

if kind[0].lower() == 'a': # for ADEV
exponent = 1.0
elif kind[0].lower() == 'm': # for modADEV
exponent = 2.0

integrand = np.array([
S_y *
np.abs(np.sin(np.pi * f * taus_used[idx])**(exponent + 1.0)
/ (np.pi * f * taus_used[idx])**exponent)**2.0
for idx, mj in enumerate(m)])
integrand = np.insert(integrand, 0, 0.0, axis=1)
f = np.insert(f, 0, 0.0)
ad = np.sqrt(2.0 * simps(integrand, f))
return taus_used, ad[/CODE]

 

[jedishrfu]

Caveat emptor on that code as there is likely some hidden issue preventing it from being released.

 

[Twigg]

Thanks for this!

I ended up writing something quick and dirty from scratch as a stand-in until I find more reliable code, but I think this code you found will have some good hints. I had a few artifacts for large tau that maybe this code addresses. I'll check it out tomorrow.

Thanks again!