A Software for Converting Spectral Density to Modified Allan Deviation

Twigg
Science Advisor
Gold Member
Messages
893
Reaction score
483
TL;DR Summary
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! :headbang:). 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!
 
Physics news on Phys.org
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.
 
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:
Code:
pip install allantools==2019.9
python
>>> import allantools
>>> allantools.psd2allan()

and I get as output:
AttributeError: module 'allantools' has no attribute 'psd2allan'
 
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]
 
  • Informative
Likes Twigg
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]
 
  • Like
  • Informative
Likes Mark44 and Twigg
Caveat emptor on that code as there is likely some hidden issue preventing it from being released.
 
  • Like
Likes 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!
 
  • Like
Likes jedishrfu
Back
Top