Software for Converting Spectral Density to Modified Allan Deviation

[Twigg]

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! ). 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.

“””
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)
    
    ...
    
“””

 

[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.

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

 

[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!