Reference

Python module for continuous wavelet spectral analysis.

This module includes a collection of routines for wavelet transform and statistical analysis via FFT algorithm. In addition, the module also includes cross-wavelet transforms, wavelet coherence tests and sample scripts.

Functions

• cwt : Continuous wavelet transform.
• icwt: Inverse continuous wavelet transform.
• significance : Significance test for the one dimensional wavelet transform.
• xwt : cross-wavelet transform.
• wct : Wavelet coherence transform.
• wct_significance : Wavelet coherence significance using Monte Carlo simulations.

Classes

• Morlet : Morlet wavelet
• Paul : Paul wavelet
• DOG : Derivative of a Guassian wavelet family
• MexicanHat : Mexican hat wavelet

Disclaimer

This module is based on routines provided by C. Torrence and G. P. Compo available at http://paos.colorado.edu/research/wavelets/, on routines provided by A. Grinsted, J. Moore and S. Jevrejeva available at http://noc.ac.uk/using-science/crosswavelet-wavelet-coherence, and on routines provided by A. Brazhe available at http://cell.biophys.msu.ru/static/swan/.

This software is released under a BSD-style open source license. Please read the license file for further information. This routine is provided as is without any express or implied warranties whatsoever.

Acknowledgements

We would like to thank Christopher Torrence, Gilbert P. Compo, Aslak Grinsted, John Moore, Svetlana Jevrejevaand and Alexey Brazhe for their code and also Jack Ireland and Renaud Dussurget for their attentive eyes, feedback and debugging.

Authors

Sebastian Krieger, Nabil Freij, and contributors.

References

1. Torrence, C. and Compo, G. P.. A Practical Guide to Wavelet Analysis. Bulletin of the American Meteorological Society, American Meteorological Society, 1998, 79, 61-78. http://dx.doi.org/10.1175/1520-0477(1998)079\<0061:APGTWA>2.0.CO;2>
2. Torrence, C. and Webster, P. J.. Interdecadal changes in the ENSO-Monsoon system, Journal of Climate, 1999, 12(8), 2679-2690. http://dx.doi.org/10.1175/1520-0442(1999)012\<2679:ICITEM>2.0.CO;2>
3. Grinsted, A.; Moore, J. C. & Jevrejeva, S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes in Geophysics, 2004, 11, 561-566. http://dx.doi.org/10.5194/npg-11-561-2004
4. Mallat, S.. A wavelet tour of signal processing: The sparse way. Academic Press, 2008, 805.
5. Addison, P. S. The illustrated wavelet transform handbook: introductory theory and applications in science, engineering, medicine and finance. IOP Publishing, 2002. http://dx.doi.org/10.1201/9781420033397
6. Liu, Y., Liang, X. S. and Weisberg, R. H. Rectification of the bias in the wavelet power spectrum. Journal of Atmospheric and Oceanic Technology, 2007, 24, 2093-2102. http://dx.doi.org/10.1175/2007JTECHO511.1

DOG

Bases: object

Implements the derivative of a Guassian wavelet class.

Note that the input parameter f is the angular frequency and that for m=2 the DOG becomes the Mexican hat wavelet, which is then default.

Source code in src/pycwt/mothers.py

class DOG(object):
    """Implements the derivative of a Guassian wavelet class.

    Note that the input parameter f is the angular frequency and that
    for m=2 the DOG becomes the Mexican hat wavelet, which is then
    default.

    """

    def __init__(self, m=2):
        self._set_m(m)
        self.name = "DOG"

    def psi_ft(self, f):
        """Fourier transform of the DOG wavelet."""
        return (
            -(1j**self.m)
            / numpy.sqrt(gamma(self.m + 0.5))
            * f**self.m
            * numpy.exp(-0.5 * f**2)
        )

    def psi(self, t):
        """DOG wavelet as described in Torrence and Compo (1998).

        The derivative of a Gaussian of order `n` can be determined using
        the probabilistic Hermite polynomials. They are explicitly
        written as:
            Hn(x) = 2 ** (-n / s) * n! * sum ((-1) ** m) *
                    (2 ** 0.5 * x) ** (n - 2 * m) / (m! * (n - 2*m)!)
        or in the recursive form:
            Hn(x) = x * Hn(x) - nHn-1(x)

        Source: http://www.ask.com/wiki/Hermite_polynomials

        """
        p = hermitenorm(self.m)
        return (
            (-1) ** (self.m + 1)
            * numpy.polyval(p, t)
            * numpy.exp(-(t**2) / 2)
            / numpy.sqrt(gamma(self.m + 0.5))
        )

    def flambda(self):
        """Fourier wavelength as of Torrence and Compo (1998)."""
        return 2 * numpy.pi / numpy.sqrt(self.m + 0.5)

    def coi(self):
        """e-Folding Time as of Torrence and Compo (1998)."""
        return 1 / numpy.sqrt(2)

    def sup(self):
        """Wavelet support defined by the e-Folding time."""
        return 1 / self.coi

    def _set_m(self, m):
        # Sets the m derivative of a Gaussian, the degrees of freedom and the
        # empirically derived factors for the wavelet bases C_{\delta},
        # \gamma, \delta j_0 (Torrence and Compo, 1998, Table 2).
        self.m = m  # m-derivative
        self.dofmin = 1  # Minimum degrees of freedom
        if self.m == 2:
            self.cdelta = 3.541  # Reconstruction factor
            self.gamma = 1.43  # Decorrelation factor for time averaging
            self.deltaj0 = 1.40  # Factor for scale averaging
        elif self.m == 6:
            self.cdelta = 1.966
            self.gamma = 1.37
            self.deltaj0 = 0.97
        else:
            self.cdelta = -1
            self.gamma = -1
            self.deltaj0 = -1

coi()

e-Folding Time as of Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def coi(self):
    """e-Folding Time as of Torrence and Compo (1998)."""
    return 1 / numpy.sqrt(2)

flambda()

Fourier wavelength as of Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def flambda(self):
    """Fourier wavelength as of Torrence and Compo (1998)."""
    return 2 * numpy.pi / numpy.sqrt(self.m + 0.5)

psi(t)

DOG wavelet as described in Torrence and Compo (1998).

The derivative of a Gaussian of order n can be determined using the probabilistic Hermite polynomials. They are explicitly written as: Hn(x) = 2 ** (-n / s) * n! * sum ((-1) ** m) * (2 ** 0.5 * x) ** (n - 2 * m) / (m! * (n - 2*m)!) or in the recursive form: Hn(x) = x * Hn(x) - nHn-1(x)

Source: http://www.ask.com/wiki/Hermite_polynomials

Source code in src/pycwt/mothers.py

def psi(self, t):
    """DOG wavelet as described in Torrence and Compo (1998).

    The derivative of a Gaussian of order `n` can be determined using
    the probabilistic Hermite polynomials. They are explicitly
    written as:
        Hn(x) = 2 ** (-n / s) * n! * sum ((-1) ** m) *
                (2 ** 0.5 * x) ** (n - 2 * m) / (m! * (n - 2*m)!)
    or in the recursive form:
        Hn(x) = x * Hn(x) - nHn-1(x)

    Source: http://www.ask.com/wiki/Hermite_polynomials

    """
    p = hermitenorm(self.m)
    return (
        (-1) ** (self.m + 1)
        * numpy.polyval(p, t)
        * numpy.exp(-(t**2) / 2)
        / numpy.sqrt(gamma(self.m + 0.5))
    )

psi_ft(f)

Fourier transform of the DOG wavelet.

Source code in src/pycwt/mothers.py

def psi_ft(self, f):
    """Fourier transform of the DOG wavelet."""
    return (
        -(1j**self.m)
        / numpy.sqrt(gamma(self.m + 0.5))
        * f**self.m
        * numpy.exp(-0.5 * f**2)
    )

sup()

Wavelet support defined by the e-Folding time.

Source code in src/pycwt/mothers.py

def sup(self):
    """Wavelet support defined by the e-Folding time."""
    return 1 / self.coi

MexicanHat

Bases: DOG

Implements the Mexican hat wavelet class.

This class inherits the DOG class using m=2.

Source code in src/pycwt/mothers.py

class MexicanHat(DOG):
    """Implements the Mexican hat wavelet class.

    This class inherits the DOG class using m=2.

    """

    def __init__(self):
        self.name = "Mexican Hat"
        self._set_m(2)

Morlet

Bases: object

Implements the Morlet wavelet class.

Note that the input parameters f and f0 are angular frequencies. f0 should be more than 0.8 for this function to be correct, its default value is f0 = 6.

Source code in src/pycwt/mothers.py

class Morlet(object):
    """Implements the Morlet wavelet class.

    Note that the input parameters f and f0 are angular frequencies.
    f0 should be more than 0.8 for this function to be correct, its
    default value is f0 = 6.

    """

    def __init__(self, f0=6):
        self._set_f0(f0)
        self.name = "Morlet"

    def psi_ft(self, f):
        """Fourier transform of the approximate Morlet wavelet."""
        return (numpy.pi**-0.25) * numpy.exp(-0.5 * (f - self.f0) ** 2)

    def psi(self, t):
        """Morlet wavelet as described in Torrence and Compo (1998)."""
        return (numpy.pi**-0.25) * numpy.exp(1j * self.f0 * t - t**2 / 2)

    def flambda(self):
        """Fourier wavelength as of Torrence and Compo (1998)."""
        return (4 * numpy.pi) / (self.f0 + numpy.sqrt(2 + self.f0**2))

    def coi(self):
        """e-Folding Time as of Torrence and Compo (1998)."""
        return 1.0 / numpy.sqrt(2)

    def sup(self):
        """Wavelet support defined by the e-Folding time."""
        return 1.0 / self.coi

    def _set_f0(self, f0):
        # Sets the Morlet wave number, the degrees of freedom and the
        # empirically derived factors for the wavelet bases C_{\delta},
        # \gamma, \delta j_0 (Torrence and Compo, 1998, Table 2)
        self.f0 = f0  # Wave number
        self.dofmin = 2  # Minimum degrees of freedom
        if self.f0 == 6:
            self.cdelta = 0.776  # Reconstruction factor
            self.gamma = 2.32  # Decorrelation factor for time averaging
            self.deltaj0 = 0.60  # Factor for scale averaging
        else:
            self.cdelta = -1
            self.gamma = -1
            self.deltaj0 = -1

    def smooth(self, W, dt, dj, scales):
        """Smoothing function used in coherence analysis.

        Parameters
        ----------
        W :
        dt :
        dj :
        scales :

        Returns
        -------
        T :

        """
        # The smoothing is performed by using a filter given by the absolute
        # value of the wavelet function at each scale, normalized to have a
        # total weight of unity, according to suggestions by Torrence &
        # Webster (1999) and by Grinsted et al. (2004).
        m, n = W.shape

        # Filter in time.
        k = 2 * numpy.pi * fft.fftfreq(fft_kwargs(W[0, :])["n"])
        k2 = k**2
        snorm = scales / dt
        # Smoothing by Gaussian window (absolute value of wavelet function)
        # using the convolution theorem: multiplication by Gaussian curve in
        # Fourier domain for each scale, outer product of scale and frequency
        F = numpy.exp(-0.5 * (snorm[:, numpy.newaxis] ** 2) * k2)  # Outer product
        smooth = fft.ifft(
            F * fft.fft(W, axis=1, **fft_kwargs(W[0, :])),
            axis=1,  # Along Fourier frequencies
            **fft_kwargs(W[0, :], overwrite_x=True),
        )
        T = smooth[:, :n]  # Remove possibly padded region due to FFT

        if numpy.isreal(W).all():
            T = T.real

        # Filter in scale. For the Morlet wavelet it's simply a boxcar with
        # 0.6 width.
        wsize = self.deltaj0 / dj * 2
        win = rect(int(numpy.round(wsize)), normalize=True)
        T = convolve2d(T, win[:, numpy.newaxis], "same")  # Scales are "vertical"

        return T

coi()

e-Folding Time as of Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def coi(self):
    """e-Folding Time as of Torrence and Compo (1998)."""
    return 1.0 / numpy.sqrt(2)

flambda()

Fourier wavelength as of Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def flambda(self):
    """Fourier wavelength as of Torrence and Compo (1998)."""
    return (4 * numpy.pi) / (self.f0 + numpy.sqrt(2 + self.f0**2))

psi(t)

Morlet wavelet as described in Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def psi(self, t):
    """Morlet wavelet as described in Torrence and Compo (1998)."""
    return (numpy.pi**-0.25) * numpy.exp(1j * self.f0 * t - t**2 / 2)

psi_ft(f)

Fourier transform of the approximate Morlet wavelet.

Source code in src/pycwt/mothers.py

def psi_ft(self, f):
    """Fourier transform of the approximate Morlet wavelet."""
    return (numpy.pi**-0.25) * numpy.exp(-0.5 * (f - self.f0) ** 2)

smooth(W, dt, dj, scales)

Smoothing function used in coherence analysis.

Parameters

W : dt : dj : scales :

Returns

T :

Source code in src/pycwt/mothers.py

def smooth(self, W, dt, dj, scales):
    """Smoothing function used in coherence analysis.

    Parameters
    ----------
    W :
    dt :
    dj :
    scales :

    Returns
    -------
    T :

    """
    # The smoothing is performed by using a filter given by the absolute
    # value of the wavelet function at each scale, normalized to have a
    # total weight of unity, according to suggestions by Torrence &
    # Webster (1999) and by Grinsted et al. (2004).
    m, n = W.shape

    # Filter in time.
    k = 2 * numpy.pi * fft.fftfreq(fft_kwargs(W[0, :])["n"])
    k2 = k**2
    snorm = scales / dt
    # Smoothing by Gaussian window (absolute value of wavelet function)
    # using the convolution theorem: multiplication by Gaussian curve in
    # Fourier domain for each scale, outer product of scale and frequency
    F = numpy.exp(-0.5 * (snorm[:, numpy.newaxis] ** 2) * k2)  # Outer product
    smooth = fft.ifft(
        F * fft.fft(W, axis=1, **fft_kwargs(W[0, :])),
        axis=1,  # Along Fourier frequencies
        **fft_kwargs(W[0, :], overwrite_x=True),
    )
    T = smooth[:, :n]  # Remove possibly padded region due to FFT

    if numpy.isreal(W).all():
        T = T.real

    # Filter in scale. For the Morlet wavelet it's simply a boxcar with
    # 0.6 width.
    wsize = self.deltaj0 / dj * 2
    win = rect(int(numpy.round(wsize)), normalize=True)
    T = convolve2d(T, win[:, numpy.newaxis], "same")  # Scales are "vertical"

    return T

sup()

Wavelet support defined by the e-Folding time.

Source code in src/pycwt/mothers.py

def sup(self):
    """Wavelet support defined by the e-Folding time."""
    return 1.0 / self.coi

Paul

Bases: object

Implements the Paul wavelet class.

Note that the input parameter f is the angular frequency and that the default order for this wavelet is m=4.

Source code in src/pycwt/mothers.py

class Paul(object):
    """Implements the Paul wavelet class.

    Note that the input parameter f is the angular frequency and that
    the default order for this wavelet is m=4.

    """

    def __init__(self, m=4):
        self._set_m(m)
        self.name = "Paul"

    def psi_ft(self, f):
        """Fourier transform of the Paul wavelet."""
        return (
            2**self.m
            / numpy.sqrt(self.m * numpy.prod(range(2, 2 * self.m)))
            * f**self.m
            * numpy.exp(-f)
            * (f > 0)
        )

    def psi(self, t):
        """Paul wavelet as described in Torrence and Compo (1998)."""
        return (
            2**self.m
            * 1j**self.m
            * numpy.prod(range(2, self.m - 1))
            / numpy.sqrt(numpy.pi * numpy.prod(range(2, 2 * self.m + 1)))
            * (1 - 1j * t) ** (-(self.m + 1))
        )

    def flambda(self):
        """Fourier wavelength as of Torrence and Compo (1998)."""
        return 4 * numpy.pi / (2 * self.m + 1)

    def coi(self):
        """e-Folding Time as of Torrence and Compo (1998)."""
        return numpy.sqrt(2)

    def sup(self):
        """Wavelet support defined by the e-Folding time."""
        return 1 / self.coi

    def _set_m(self, m):
        # Sets the m derivative of a Gaussian, the degrees of freedom and the
        # empirically derived factors for the wavelet bases C_{\delta},
        # \gamma, \delta j_0 (Torrence and Compo, 1998, Table 2)
        self.m = m  # Wavelet order
        self.dofmin = 2  # Minimum degrees of freedom
        if self.m == 4:
            self.cdelta = 1.132  # Reconstruction factor
            self.gamma = 1.17  # Decorrelation factor for time averaging
            self.deltaj0 = 1.50  # Factor for scale averaging
        else:
            self.cdelta = -1
            self.gamma = -1
            self.deltaj0 = -1

coi()

e-Folding Time as of Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def coi(self):
    """e-Folding Time as of Torrence and Compo (1998)."""
    return numpy.sqrt(2)

flambda()

Fourier wavelength as of Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def flambda(self):
    """Fourier wavelength as of Torrence and Compo (1998)."""
    return 4 * numpy.pi / (2 * self.m + 1)

psi(t)

Paul wavelet as described in Torrence and Compo (1998).

Source code in src/pycwt/mothers.py

def psi(self, t):
    """Paul wavelet as described in Torrence and Compo (1998)."""
    return (
        2**self.m
        * 1j**self.m
        * numpy.prod(range(2, self.m - 1))
        / numpy.sqrt(numpy.pi * numpy.prod(range(2, 2 * self.m + 1)))
        * (1 - 1j * t) ** (-(self.m + 1))
    )

psi_ft(f)

Fourier transform of the Paul wavelet.

Source code in src/pycwt/mothers.py

def psi_ft(self, f):
    """Fourier transform of the Paul wavelet."""
    return (
        2**self.m
        / numpy.sqrt(self.m * numpy.prod(range(2, 2 * self.m)))
        * f**self.m
        * numpy.exp(-f)
        * (f > 0)
    )

sup()

Wavelet support defined by the e-Folding time.

Source code in src/pycwt/mothers.py

def sup(self):
    """Wavelet support defined by the e-Folding time."""
    return 1 / self.coi

ar1(x)

Lag-1 autocorrelation coefficient.

In an Allen and Smith AR(1) model,

\[ x(t) - <x> = \gamma(x(t-1) - <x>) + \alpha z(t) ,\]

where \(<x>\) is the process mean, \(\gamma\) and \(\alpha\) are process parameters and \(z(t)\) is a Gaussian unit-variance white noise.

Parameters

x : numpy.ndarray, list Univariate time series

Returns

g : float Estimate of the lag-one autocorrelation. a : float Estimate of the noise variance [var(x) ~= a2/(1-g2)] mu2 : float Estimated square on the mean of a finite segment of AR(1) noise, mormalized by the process variance.

References

[1] Allen, M. R. and Smith, L. A. Monte Carlo SSA: detecting irregular oscillations in the presence of colored noise. Journal of Climate, 1996, 9(12), 3373-3404. <http://dx.doi.org/10.1175/1520-0442(1996)009<3373:MCSDIO>2.0.CO;2> [2] http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html

Source code in src/pycwt/helpers.py

def ar1(x):
    r"""Lag-1 autocorrelation coefficient.

    In an Allen and Smith AR(1) model,

    $$  x(t) - <x> = \gamma(x(t-1) - <x>) + \alpha z(t) ,$$

    where $<x>$ is the process mean, $\gamma$ and $\alpha$ are process
    parameters and $z(t)$ is a Gaussian unit-variance white noise.

    Parameters
    ----------
    x : numpy.ndarray, list
        Univariate time series

    Returns
    -------
    g : float
        Estimate of the lag-one autocorrelation.
    a : float
        Estimate of the noise variance [var(x) ~= a**2/(1-g**2)]
    mu2 : float
        Estimated square on the mean of a finite segment of AR(1)
        noise, mormalized by the process variance.

    References
    ----------
    [1] Allen, M. R. and Smith, L. A. Monte Carlo SSA: detecting
        irregular oscillations in the presence of colored noise.
        *Journal of Climate*, **1996**, 9(12), 3373-3404.
        <http://dx.doi.org/10.1175/1520-0442(1996)009<3373:MCSDIO>2.0.CO;2>
    [2] http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html

    """
    x = numpy.asarray(x)
    N = x.size
    xm = x.mean()
    x = x - xm

    # Estimates the lag zero and one covariance
    c0 = x.transpose().dot(x) / N
    c1 = x[0 : N - 1].transpose().dot(x[1:N]) / (N - 1)

    # According to A. Grinsteds' substitutions
    B = -c1 * N - c0 * N**2 - 2 * c0 + 2 * c1 - c1 * N**2 + c0 * N
    A = c0 * N**2
    C = N * (c0 + c1 * N - c1)
    D = B**2 - 4 * A * C

    if D > 0:
        g = (-B - D**0.5) / (2 * A)
    else:
        raise Warning(
            "Cannot place an upperbound on the unbiased AR(1). "
            "Series is too short or trend is too large."
        )

    # According to Allen & Smith (1996), footnote 4
    mu2 = -1 / N + (2 / N**2) * (
        (N - g**N) / (1 - g) - g * (1 - g ** (N - 1)) / (1 - g) ** 2
    )
    c0t = c0 / (1 - mu2)
    a = ((1 - g**2) * c0t) ** 0.5

    return g, a, mu2

ar1_spectrum(freqs, ar1=0.0)

Lag-1 autoregressive theoretical power spectrum.

Parameters

freqs : numpy.ndarray, list Frequencies at which to calculate the theoretical power spectrum. ar1 : float Autoregressive lag-1 correlation coefficient.

Returns

Pk : numpy.ndarray Theoretical discrete Fourier power spectrum of noise signal.

References

[1] http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html

Source code in src/pycwt/helpers.py

def ar1_spectrum(freqs, ar1=0.0):
    """
    Lag-1 autoregressive theoretical power spectrum.

    Parameters
    ----------
    freqs : numpy.ndarray, list
        Frequencies at which to calculate the theoretical power
        spectrum.
    ar1 : float
        Autoregressive lag-1 correlation coefficient.

    Returns
    -------
    Pk : numpy.ndarray
        Theoretical discrete Fourier power spectrum of noise signal.

    References
    ----------
    [1] http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html

    """
    # According to a post from the MadSci Network available at
    # http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html,
    # the time-series spectrum for an auto-regressive model can be
    # represented as
    #
    # P_k = \frac{E}{\left|1- \sum\limits_{k=1}^{K} a_k \, e^{2 i \pi
    #   \frac{k f}{f_s} } \right|^2}
    #
    # which for an AR1 model reduces to
    #
    freqs = numpy.asarray(freqs)
    Pk = (1 - ar1**2) / numpy.abs(
        1 - ar1 * numpy.exp(-2 * numpy.pi * 1j * freqs)
    ) ** 2

    return Pk

cwt(signal, dt, dj=0.08333333333333333, s0=-1, J=-1, wavelet='morlet', freqs=None)

Continuous wavelet transform of the signal at specified scales.

Parameters

signal : numpy.ndarray, list Input signal array. dt : float Sampling interval. dj : float, optional Spacing between discrete scales. Default value is 1/12. Smaller values will result in better scale resolution, but slower calculation and plot. s0 : float, optional Smallest scale of the wavelet. Default value is 2dt. J : float, optional Number of scales less one. Scales range from s0 up to s0 * 2*(J * dj), which gives a total of (J + 1) scales. Default is J = (log2(N * dt / s0)) / dj. wavelet : instance of Wavelet class, or string Mother wavelet class. Default is Morlet wavelet. freqs : numpy.ndarray, optional Custom frequencies to use instead of the ones corresponding to the scales described above. Corresponding scales are calculated using the wavelet Fourier wavelength.

Returns

W : numpy.ndarray Wavelet transform according to the selected mother wavelet. Has (J+1) x N dimensions. sj : numpy.ndarray Vector of scale indices given by sj = s0 * 2**(j * dj), j={0, 1, ..., J}. freqs : array like Vector of Fourier frequencies (in 1 / time units) that corresponds to the wavelet scales. coi : numpy.ndarray Returns the cone of influence, which is a vector of N points containing the maximum Fourier period of useful information at that particular time. Periods greater than those are subject to edge effects. fft : numpy.ndarray Normalized fast Fourier transform of the input signal. fftfreqs : numpy.ndarray Fourier frequencies (in 1/time units) for the calculated FFT spectrum.

Example

mother = wavelet.Morlet(6.) wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(signal, 0.25, 0.25, 0.5, 28, mother)

fft_kwargs(signal, **kwargs)

Return next higher power of 2 for given signal to speed up FFT

Source code in src/pycwt/helpers.py

def fft_kwargs(signal, **kwargs):
    """Return next higher power of 2 for given signal to speed up FFT"""
    if _FFT_NEXT_POW2:
        return {"n": int(2 ** numpy.ceil(numpy.log2(len(signal))))}

find(condition)

Returns the indices where ravel(condition) is true.

Source code in src/pycwt/helpers.py

def find(condition):
    """Returns the indices where ravel(condition) is true."""
    (res,) = numpy.nonzero(numpy.ravel(condition))
    return res

get_cache_dir()

Returns the location of the cache directory.

Source code in src/pycwt/helpers.py

def get_cache_dir():
    """Returns the location of the cache directory."""
    # Sets cache directory according to user home path.
    cache_dir = "{}/.cache/pycwt/".format(expanduser("~"))
    # Creates cache directory if not existant.
    if not exists(cache_dir):
        makedirs(cache_dir)
    # Returns cache directory.
    return cache_dir

icwt(W, sj, dt, dj=0.08333333333333333, wavelet='morlet')

Inverse continuous wavelet transform.

Parameters

W : numpy.ndarray Wavelet transform, the result of the cwt function. sj : numpy.ndarray Vector of scale indices as returned by the cwt function. dt : float Sample spacing. dj : float, optional Spacing between discrete scales as used in the cwt function. Default value is 0.25. wavelet : instance of Wavelet class, or string Mother wavelet class. Default is Morlet

Returns

iW : numpy.ndarray Inverse wavelet transform.

Example

mother = wavelet.Morlet() wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(var, 0.25, 0.25, 0.5, 28, mother) iwave = wavelet.icwt(wave, scales, 0.25, 0.25, mother)

rednoise(N, g, a=1.0)

Red noise generator using filter.

Parameters

N : int Length of the desired time series. g : float Lag-1 autocorrelation coefficient. a : float, optional Noise innovation variance parameter.

Returns

y : numpy.ndarray Red noise time series.

Source code in src/pycwt/helpers.py

def rednoise(N, g, a=1.0):
    """
    Red noise generator using filter.

    Parameters
    ----------
    N : int
        Length of the desired time series.
    g : float
        Lag-1 autocorrelation coefficient.
    a : float, optional
        Noise innovation variance parameter.

    Returns
    -------
    y : numpy.ndarray
        Red noise time series.

    """
    if g == 0:
        yr = numpy.randn(N, 1) * a
    else:
        # Twice the decorrelation time.
        tau = int(numpy.ceil(-2 / numpy.log(numpy.abs(g))))
        yr = lfilter([1, 0], [1, -g], numpy.random.randn(N + tau, 1) * a)
        yr = yr[tau:]

    return yr.flatten()

significance(signal, dt, scales, sigma_test=0, alpha=None, significance_level=0.95, dof=-1, wavelet='morlet')

Significance test for the one dimensional wavelet transform.

Parameters

signal : array like, float Input signal array. If a float number is given, then the variance is assumed to have this value. If an array is given, then its variance is automatically computed. dt : float Sample spacing. scales : array like Vector of scale indices given returned by cwt function. sigma_test : int, optional Sets the type of significance test to be performed. Accepted values are 0 (default), 1 or 2. See notes below for further details. alpha : float, optional Lag-1 autocorrelation, used for the significance levels. Default is 0.0. significance_level : float, optional Significance level to use. Default is 0.95. dof : variant, optional Degrees of freedom for significance test to be set according to the type set in sigma_test. wavelet : instance of Wavelet class, or string Mother wavelet class. Default is Morlet

Returns

signif : array like Significance levels as a function of scale. fft_theor (array like): Theoretical red-noise spectrum as a function of period.

Notes

If sigma_test is set to 0, performs a regular chi-square test, according to Torrence and Compo (1998) equation 18.

If set to 1, performs a time-average test (equation 23). In this case, dof should be set to the number of local wavelet spectra that where averaged together. For the global wavelet spectra it would be dof=N, the number of points in the time-series.

If set to 2, performs a scale-average test (equations 25 to 28). In this case dof should be set to a two element vector [s1, s2], which gives the scale range that were averaged together. If, for example, the average between scales 2 and 8 was taken, then dof=[2, 8].

wct(y1, y2, dt, dj=0.08333333333333333, s0=-1, J=-1, sig=True, significance_level=0.95, wavelet='morlet', normalize=True, **kwargs)

Wavelet coherence transform (WCT).

The WCT finds regions in time frequency space where the two time series co-vary, but do not necessarily have high power.

Parameters

y1, y2 : numpy.ndarray, list Input signals. dt : float Sample spacing. dj : float, optional Spacing between discrete scales. Default value is 1/12. Smaller values will result in better scale resolution, but slower calculation and plot. s0 : float, optional Smallest scale of the wavelet. Default value is 2dt. J : float, optional Number of scales less one. Scales range from s0 up to s0 * 2(J * dj), which gives a total of (J + 1) scales. Default is J = (log2(Ndt/so))/dj. sig : bool set to compute signficance, default is True significance_level (float, optional) : Significance level to use. Default is 0.95. normalize (boolean, optional) : If set to true, normalizes CWT by the standard deviation of the signals.

Returns

WCT : magnitude of coherence aWCT : phase angle of coherence coi (array like): Cone of influence, which is a vector of N points containing the maximum Fourier period of useful information at that particular time. Periods greater than those are subject to edge effects. freq (array like): Vector of Fourier equivalent frequencies (in 1 / time units) coi : sig : Significance levels as a function of scale if sig=True when called, otherwise zero.

See also

cwt, xwt

wct_significance(al1, al2, dt, dj, s0, J, significance_level=0.95, wavelet='morlet', mc_count=300, progress=True, cache=True)

Wavelet coherence transform significance.

Calculates WCT significance using Monte Carlo simulations with 95% confidence.

Parameters

al1, al2: float Lag-1 autoregressive coeficients of both time series. dt : float Sample spacing. dj : float, optional Spacing between discrete scales. Default value is 1/12. Smaller values will result in better scale resolution, but slower calculation and plot. s0 : float, optional Smallest scale of the wavelet. Default value is 2dt. J : float, optional Number of scales less one. Scales range from s0 up to s0 * 2(J * dj), which gives a total of (J + 1) scales. Default is J = (log2(Ndt/so))/dj. significance_level : float, optional Significance level to use. Default is 0.95. wavelet : instance of a wavelet class, optional Mother wavelet class. Default is Morlet wavelet. mc_count : integer, optional Number of Monte Carlo simulations. Default is 300. progress : bool, optional If True (default), shows progress bar on screen. cache : bool, optional If True (default) saves cache to file.

Returns

TODO

xwt(y1, y2, dt, dj=0.08333333333333333, s0=-1, J=-1, significance_level=0.95, wavelet='morlet', normalize=True)

Cross wavelet transform (XWT) of two signals.

The XWT finds regions in time frequency space where the time series show high common power.

Parameters

y1, y2 : numpy.ndarray, list Input signal array to calculate cross wavelet transform. dt : float Sample spacing. dj : float, optional Spacing between discrete scales. Default value is 1/12. Smaller values will result in better scale resolution, but slower calculation and plot. s0 : float, optional Smallest scale of the wavelet. Default value is 2dt. J : float, optional Number of scales less one. Scales range from s0 up to s0 * 2(J * dj), which gives a total of (J + 1) scales. Default is J = (log2(Ndt/so))/dj. wavelet : instance of a wavelet class, optional Mother wavelet class. Default is Morlet wavelet. significance_level : float, optional Significance level to use. Default is 0.95. normalize : bool, optional If set to true, normalizes CWT by the standard deviation of the signals.

Returns

xwt (array like): Cross wavelet transform according to the selected mother wavelet. x (array like): Intersected independent variable. coi (array like): Cone of influence, which is a vector of N points containing the maximum Fourier period of useful information at that particular time. Periods greater than those are subject to edge effects. freqs (array like): Vector of Fourier equivalent frequencies (in 1 / time units) that correspond to the wavelet scales. signif (array like): Significance levels as a function of scale.

Notes

Torrence and Compo (1998) state that the percent point function (PPF) -- inverse of the cumulative distribution function -- of a chi-square distribution at 95% confidence and two degrees of freedom is Z2(95%)=3.999. However, calculating the PPF using chi2.ppf gives Z2(95%)=5.991. To ensure similar significance intervals as in Grinsted et al. (2004), one has to use confidence of 86.46%.