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, Alexey Brazhe, Christopher Torrence, Gilbert P. Compo 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 pycwt/mothers.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
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 / np.sqrt(gamma(self.m + 0.5)) * f ** self.m *
                np.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) * np.polyval(p, t) *
                np.exp(-t ** 2 / 2) / np.sqrt(gamma(self.m + 0.5)))

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

    def coi(self):
        """e-Folding Time as of Torrence and Compo (1998)."""
        return 1 / np.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 pycwt/mothers.py 
197
198
199
def coi(self):
    """e-Folding Time as of Torrence and Compo (1998)."""
    return 1 / np.sqrt(2)

flambda()

Fourier wavelength as of Torrence and Compo (1998).

Source code in pycwt/mothers.py 
193
194
195
def flambda(self):
    """Fourier wavelength as of Torrence and Compo (1998)."""
    return (2 * np.pi / np.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 pycwt/mothers.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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) * np.polyval(p, t) *
            np.exp(-t ** 2 / 2) / np.sqrt(gamma(self.m + 0.5)))

psi_ft(f)

Fourier transform of the DOG wavelet.

Source code in pycwt/mothers.py 
170
171
172
173
def psi_ft(self, f):
    """Fourier transform of the DOG wavelet."""
    return (- 1j ** self.m / np.sqrt(gamma(self.m + 0.5)) * f ** self.m *
            np.exp(- 0.5 * f ** 2))

sup()

Wavelet support defined by the e-Folding time.

Source code in pycwt/mothers.py 
201
202
203
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 pycwt/mothers.py 
225
226
227
228
229
230
231
232
233
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 pycwt/mothers.py 
 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
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 (np.pi ** -0.25) * np.exp(-0.5 * (f - self.f0) ** 2)

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

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

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

    def sup(self):
        """Wavelet support defined by the e-Folding time."""
        return 1. / 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 * np.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 = np.exp(-0.5 * (snorm[:, np.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 np.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(np.round(wsize)), normalize=True)
        T = convolve2d(T, win[:, np.newaxis], 'same')  # Scales are "vertical"

        return T

coi()

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

Source code in pycwt/mothers.py 
38
39
40
def coi(self):
    """e-Folding Time as of Torrence and Compo (1998)."""
    return 1. / np.sqrt(2)

flambda()

Fourier wavelength as of Torrence and Compo (1998).

Source code in pycwt/mothers.py 
34
35
36
def flambda(self):
    """Fourier wavelength as of Torrence and Compo (1998)."""
    return (4 * np.pi) / (self.f0 + np.sqrt(2 + self.f0 ** 2))

psi(t)

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

Source code in pycwt/mothers.py 
30
31
32
def psi(self, t):
    """Morlet wavelet as described in Torrence and Compo (1998)."""
    return (np.pi ** -0.25) * np.exp(1j * self.f0 * t - t ** 2 / 2)

psi_ft(f)

Fourier transform of the approximate Morlet wavelet.

Source code in pycwt/mothers.py 
26
27
28
def psi_ft(self, f):
    """Fourier transform of the approximate Morlet wavelet."""
    return (np.pi ** -0.25) * np.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 pycwt/mothers.py 
 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
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 * np.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 = np.exp(-0.5 * (snorm[:, np.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 np.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(np.round(wsize)), normalize=True)
    T = convolve2d(T, win[:, np.newaxis], 'same')  # Scales are "vertical"

    return T

sup()

Wavelet support defined by the e-Folding time.

Source code in pycwt/mothers.py 
42
43
44
def sup(self):
    """Wavelet support defined by the e-Folding time."""
    return 1. / 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 pycwt/mothers.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
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 /
                np.sqrt(self.m * np.prod(range(2, 2 * self.m))) *
                f ** self.m * np.exp(-f) * (f > 0))

    def psi(self, t):
        """Paul wavelet as described in Torrence and Compo (1998)."""
        return (2 ** self.m * 1j ** self.m * np.prod(range(2, self.m - 1)) /
                np.sqrt(np.pi * np.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 * np.pi / (2 * self.m + 1)

    def coi(self):
        """e-Folding Time as of Torrence and Compo (1998)."""
        return np.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 pycwt/mothers.py 
134
135
136
def coi(self):
    """e-Folding Time as of Torrence and Compo (1998)."""
    return np.sqrt(2)

flambda()

Fourier wavelength as of Torrence and Compo (1998).

Source code in pycwt/mothers.py 
130
131
132
def flambda(self):
    """Fourier wavelength as of Torrence and Compo (1998)."""
    return 4 * np.pi / (2 * self.m + 1)

psi(t)

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

Source code in pycwt/mothers.py 
124
125
126
127
128
def psi(self, t):
    """Paul wavelet as described in Torrence and Compo (1998)."""
    return (2 ** self.m * 1j ** self.m * np.prod(range(2, self.m - 1)) /
            np.sqrt(np.pi * np.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 pycwt/mothers.py 
118
119
120
121
122
def psi_ft(self, f):
    """Fourier transform of the Paul wavelet."""
    return (2 ** self.m /
            np.sqrt(self.m * np.prod(range(2, 2 * self.m))) *
            f ** self.m * np.exp(-f) * (f > 0))

sup()

Wavelet support defined by the e-Folding time.

Source code in pycwt/mothers.py 
138
139
140
def sup(self):
    """Wavelet support defined by the e-Folding time."""
    return 1 / self.coi

cwt(signal, dt, dj=1 / 12, 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 / so)) / 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)

Source code in pycwt/wavelet.py 
 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
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def cwt(signal, dt, dj=1/12, 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 2*dt.
    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 / so)) / 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)

    """
    wavelet = _check_parameter_wavelet(wavelet)

    # Original signal length
    n0 = len(signal)
    # If no custom frequencies are set, then set default frequencies
    # according to input parameters `dj`, `s0` and `J`. Otherwise, set wavelet
    # scales according to Fourier equivalent frequencies.
    if freqs is None:
        # Smallest resolvable scale
        if s0 == -1:
            s0 = 2 * dt / wavelet.flambda()
        # Number of scales
        if J == -1:
            J = int(np.round(np.log2(n0 * dt / s0) / dj))
        # The scales as of Mallat 1999
        sj = s0 * 2 ** (np.arange(0, J + 1) * dj)
        # Fourier equivalent frequencies
        freqs = 1 / (wavelet.flambda() * sj)
    else:
        # The wavelet scales using custom frequencies.
        sj = 1 / (wavelet.flambda() * freqs)

    # Signal Fourier transform
    signal_ft = fft.fft(signal, **fft_kwargs(signal))
    N = len(signal_ft)
    # Fourier angular frequencies
    ftfreqs = 2 * np.pi * fft.fftfreq(N, dt)

    # Creates wavelet transform matrix as outer product of scaled transformed
    # wavelets and transformed signal according to the convolution theorem.
    # (i)   Transform scales to column vector for outer product;
    # (ii)  Calculate 2D matrix [s, f] for each scale s and Fourier angular
    #       frequency f;
    # (iii) Calculate wavelet transform;
    sj_col = sj[:, np.newaxis]
    psi_ft_bar = ((sj_col * ftfreqs[1] * N) ** .5 *
                  np.conjugate(wavelet.psi_ft(sj_col * ftfreqs)))
    W = fft.ifft(signal_ft * psi_ft_bar, axis=1,
                 **fft_kwargs(signal_ft, overwrite_x=True))

    # Checks for NaN in transform results and removes them from the scales if
    # needed, frequencies and wavelet transform. Trims wavelet transform at
    # length `n0`.
    sel = np.invert(np.isnan(W).all(axis=1))
    if np.any(sel):
        sj = sj[sel]
        freqs = freqs[sel]
        W = W[sel, :]

    # Determines the cone-of-influence. Note that it is returned as a function
    # of time in Fourier periods. Uses triangualr Bartlett window with
    # non-zero end-points.
    coi = (n0 / 2 - np.abs(np.arange(0, n0) - (n0 - 1) / 2))
    coi = wavelet.flambda() * wavelet.coi() * dt * coi

    return (W[:, :n0], sj, freqs, coi, signal_ft[1:N//2] / N ** 0.5,
            ftfreqs[1:N//2] / (2 * np.pi))

icwt(W, sj, dt, dj=1 / 12, 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)

Source code in pycwt/wavelet.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def icwt(W, sj, dt, dj=1/12, 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)

    """
    wavelet = _check_parameter_wavelet(wavelet)

    a, b = W.shape
    c = sj.size
    if a == c:
        sj = (np.ones([b, 1]) * sj).transpose()
    elif b == c:
        sj = np.ones([a, 1]) * sj
    else:
        raise Warning('Input array dimensions do not match.')

    # As of Torrence and Compo (1998), eq. (11)
    iW = (dj * np.sqrt(dt) / (wavelet.cdelta * wavelet.psi(0)) *
          (np.real(W) / np.sqrt(sj)).sum(axis=0))
    return iW

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

Source code in pycwt/wavelet.py 
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def 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].

    """
    wavelet = _check_parameter_wavelet(wavelet)

    try:
        n0 = len(signal)
    except TypeError:
        n0 = 1
    J = len(scales) - 1
    dj = np.log2(scales[1] / scales[0])

    if n0 == 1:
        variance = signal
    else:
        variance = signal.std() ** 2

    if alpha is None:
        alpha, _, _ = ar1(signal)

    period = scales * wavelet.flambda()  # Fourier equivalent periods
    freq = dt / period                   # Normalized frequency
    dofmin = wavelet.dofmin              # Degrees of freedom with no smoothing
    Cdelta = wavelet.cdelta              # Reconstruction factor
    gamma_fac = wavelet.gamma            # Time-decorrelation factor
    dj0 = wavelet.deltaj0                # Scale-decorrelation factor

    # Theoretical discrete Fourier power spectrum of the noise signal
    # following Gilman et al. (1963) and Torrence and Compo (1998),
    # equation 16.
    def pk(k, a, N):
        return (1 - a ** 2) / (1 + a ** 2 - 2 * a * np.cos(2 * np.pi * k / N))
    fft_theor = pk(freq, alpha, n0)
    fft_theor = variance * fft_theor     # Including time-series variance
    signif = fft_theor

    try:
        if dof == -1:
            dof = dofmin
    except ValueError:
        pass

    if sigma_test == 0:  # No smoothing, dof=dofmin, TC98 sec. 4
        dof = dofmin
        # As in Torrence and Compo (1998), equation 18.
        chisquare = chi2.ppf(significance_level, dof) / dof
        signif = fft_theor * chisquare
    elif sigma_test == 1:  # Time-averaged significance
        if len(dof) == 1:
            dof = np.zeros(1, J+1) + dof
        sel = find(dof < 1)
        dof[sel] = 1
        # As in Torrence and Compo (1998), equation 23:
        dof = dofmin * (1 + (dof * dt / gamma_fac / scales) ** 2) ** 0.5
        sel = find(dof < dofmin)
        dof[sel] = dofmin  # Minimum dof is dofmin
        for n, d in enumerate(dof):
            chisquare = chi2.ppf(significance_level, d) / d
            signif[n] = fft_theor[n] * chisquare
    elif sigma_test == 2:  # Time-averaged significance
        if len(dof) != 2:
            raise Exception('DOF must be set to [s1, s2], '
                            'the range of scale-averages')
        if Cdelta == -1:
            raise ValueError('Cdelta and dj0 not defined '
                             'for {} with f0={}'.format(wavelet.name,
                                                        wavelet.f0))
        s1, s2 = dof
        sel = find((scales >= s1) & (scales <= s2))
        navg = sel.size
        if navg == 0:
            raise ValueError('No valid scales between {} and {}.'.format(s1,
                                                                         s2))
        # As in Torrence and Compo (1998), equation 25.
        Savg = 1 / sum(1. / scales[sel])
        # Power-of-two mid point:
        Smid = np.exp((np.log(s1) + np.log(s2)) / 2.)
        # As in Torrence and Compo (1998), equation 28.
        dof = (dofmin * navg * Savg / Smid) * \
              ((1 + (navg * dj / dj0) ** 2) ** 0.5)
        # As in Torrence and Compo (1998), equation 27.
        fft_theor = Savg * sum(fft_theor[sel] / scales[sel])
        chisquare = chi2.ppf(significance_level, dof) / dof
        # As in Torrence and Compo (1998), equation 26.
        signif = (dj * dt / Cdelta / Savg) * fft_theor * chisquare
    else:
        raise ValueError('sigma_test must be either 0, 1, or 2.')

    return signif, fft_theor

wct(y1, y2, dt, dj=1 / 12, 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

Source code in pycwt/wavelet.py 
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
def wct(y1, y2, dt, dj=1/12, 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 2*dt.
    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/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

    """
    wavelet = _check_parameter_wavelet(wavelet)

    # Checking some input parameters
    if s0 == -1:
        # Number of scales
        s0 = 2 * dt / wavelet.flambda()
    if J == -1:
        # Number of scales
        J = int(np.round(np.log2(y1.size * dt / s0) / dj))

    # Makes sure input signals are numpy arrays.
    y1 = np.asarray(y1)
    y2 = np.asarray(y2)
    # Calculates the standard deviation of both input signals.
    std1 = y1.std()
    std2 = y2.std()
    # Normalizes both signals, if appropriate.
    if normalize:
        y1_normal = (y1 - y1.mean()) / std1
        y2_normal = (y2 - y2.mean()) / std2
    else:
        y1_normal = y1
        y2_normal = y2

    # Calculates the CWT of the time-series making sure the same parameters
    # are used in both calculations.
    _kwargs = dict(dj=dj, s0=s0, J=J, wavelet=wavelet)
    W1, sj, freq, coi, _, _ = cwt(y1_normal, dt, **_kwargs)
    W2, sj, freq, coi, _, _ = cwt(y2_normal, dt, **_kwargs)

    scales1 = np.ones([1, y1.size]) * sj[:, None]
    scales2 = np.ones([1, y2.size]) * sj[:, None]

    # Smooth the wavelet spectra before truncating.
    S1 = wavelet.smooth(np.abs(W1) ** 2 / scales1, dt, dj, sj)
    S2 = wavelet.smooth(np.abs(W2) ** 2 / scales2, dt, dj, sj)

    # Now the wavelet transform coherence
    W12 = W1 * W2.conj()
    scales = np.ones([1, y1.size]) * sj[:, None]
    S12 = wavelet.smooth(W12 / scales, dt, dj, sj)
    WCT = np.abs(S12) ** 2 / (S1 * S2)
    aWCT = np.angle(W12)

    # Calculates the significance using Monte Carlo simulations with 95%
    # confidence as a function of scale.
    if sig:
        a1, b1, c1 = ar1(y1)
        a2, b2, c2 = ar1(y2)

        sig = wct_significance(a1, a2, dt=dt, dj=dj, s0=s0, J=J,
                               significance_level=significance_level,
                               wavelet=wavelet, **kwargs)
    else:
        sig = np.asarray([0])

    return WCT, aWCT, coi, freq, sig

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

Source code in pycwt/wavelet.py 
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
def 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 2*dt.
    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/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

    """

    if cache:
        # Load cache if previously calculated. It is assumed that wavelet
        # analysis is performed using the wavelet's default parameters.
        aa = np.round(np.arctanh(np.array([al1, al2]) * 4))
        aa = np.abs(aa) + 0.5 * (aa < 0)
        cache_file = 'wct_sig_{:0.5f}_{:0.5f}_{:0.5f}_{:0.5f}_{:d}_{}'\
            .format(aa[0], aa[1], dj, s0 / dt, J, wavelet.name)
        cache_dir = get_cache_dir()
        try:
            dat = np.loadtxt('{}/{}.gz'.format(cache_dir, cache_file),
                             unpack=True)
            print('NOTE: WCT significance loaded from cache.\n')
            return dat
        except IOError:
            pass

    # Some output to the screen
    print('Calculating wavelet coherence significance')

    # Choose N so that largest scale has at least some part outside the COI
    ms = s0 * (2 ** (J * dj)) / dt
    N = int(np.ceil(ms * 6))
    noise1 = rednoise(N, al1, 1)
    nW1, sj, freq, coi, _, _ = cwt(noise1, dt=dt, dj=dj, s0=s0, J=J,
                                   wavelet=wavelet)

    period = np.ones([1, N]) / freq[:, None]
    coi = np.ones([J + 1, 1]) * coi[None, :]
    outsidecoi = (period <= coi)
    scales = np.ones([1, N]) * sj[:, None]
    sig95 = np.zeros(J + 1)
    maxscale = find(outsidecoi.any(axis=1))[-1]
    sig95[outsidecoi.any(axis=1)] = np.nan

    nbins = 1000
    wlc = np.ma.zeros([J + 1, nbins])
    # Displays progress bar with tqdm
    for _ in tqdm(range(mc_count), disable=not progress):
        # Generates two red-noise signals with lag-1 autoregressive
        # coefficients given by al1 and al2
        noise1 = rednoise(N, al1, 1)
        noise2 = rednoise(N, al2, 1)
        # Calculate the cross wavelet transform of both red-noise signals
        kwargs = dict(dt=dt, dj=dj, s0=s0, J=J, wavelet=wavelet)
        nW1, sj, freq, coi, _, _ = cwt(noise1, **kwargs)
        nW2, sj, freq, coi, _, _ = cwt(noise2, **kwargs)
        nW12 = nW1 * nW2.conj()
        # Smooth wavelet wavelet transforms and calculate wavelet coherence
        # between both signals.
        S1 = wavelet.smooth(np.abs(nW1) ** 2 / scales, dt, dj, sj)
        S2 = wavelet.smooth(np.abs(nW2) ** 2 / scales, dt, dj, sj)
        S12 = wavelet.smooth(nW12 / scales, dt, dj, sj)
        R2 = np.ma.array(np.abs(S12) ** 2 / (S1 * S2), mask=~outsidecoi)
        # Walks through each scale outside the cone of influence and builds a
        # coherence coefficient counter.
        for s in range(maxscale):
            cd = np.floor(R2[s, :] * nbins)
            for j, t in enumerate(cd[~cd.mask]):
                wlc[s, int(t)] += 1

    # After many, many, many Monte Carlo simulations, determine the
    # significance using the coherence coefficient counter percentile.
    wlc.mask = (wlc.data == 0.)
    R2y = (np.arange(nbins) + 0.5) / nbins
    for s in range(maxscale):
        sel = ~wlc[s, :].mask
        P = wlc[s, sel].data.cumsum()
        P = (P - 0.5) / P[-1]
        sig95[s] = np.interp(significance_level, P, R2y[sel])

    if cache:
        # Save the results on cache to avoid to many computations in the future
        np.savetxt('{}/{}.gz'.format(cache_dir, cache_file), sig95)

    # And returns the results
    return sig95

xwt(y1, y2, dt, dj=1 / 12, 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%.

Source code in pycwt/wavelet.py 
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
def xwt(y1, y2, dt, dj=1/12, 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 2*dt.
    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/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%.

    """
    wavelet = _check_parameter_wavelet(wavelet)

    # Makes sure input signal are numpy arrays.
    y1 = np.asarray(y1)
    y2 = np.asarray(y2)
    # Calculates the standard deviation of both input signals.
    std1 = y1.std()
    std2 = y2.std()
    # Normalizes both signals, if appropriate.
    if normalize:
        y1_normal = (y1 - y1.mean()) / std1
        y2_normal = (y2 - y2.mean()) / std2
    else:
        y1_normal = y1
        y2_normal = y2

    # Calculates the CWT of the time-series making sure the same parameters
    # are used in both calculations.
    _kwargs = dict(dj=dj, s0=s0, J=J, wavelet=wavelet)
    W1, sj, freq, coi, _, _ = cwt(y1_normal, dt, **_kwargs)
    W2, sj, freq, coi, _, _ = cwt(y2_normal, dt, **_kwargs)

    # Calculates the cross CWT of y1 and y2.
    W12 = W1 * W2.conj()

    # And the significance tests. Note that the confidence level is calculated
    # using the percent point function (PPF) of the chi-squared cumulative
    # distribution function (CDF) instead of using Z1(95%) = 2.182 and
    # Z2(95%)=3.999 as suggested by Torrence & Compo (1998) and Grinsted et
    # al. (2004). If the CWT has been normalized, then std1 and std2 should
    # be reset to unity, otherwise the standard deviation of both series have
    # to be calculated.
    if normalize:
        std1 = std2 = 1.
    a1, _, _ = ar1(y1)
    a2, _, _ = ar1(y2)
    Pk1 = ar1_spectrum(freq * dt, a1)
    Pk2 = ar1_spectrum(freq * dt, a2)
    dof = wavelet.dofmin
    PPF = chi2.ppf(significance_level, dof)
    signif = (std1 * std2 * (Pk1 * Pk2) ** 0.5 * PPF / dof)

    # The resuts:
    return W12, coi, freq, signif