How to implement scale-dependent Gaussian averaging using Morlet wavelet envelope in Python?

51 Views Asked by At

I'm trying to reproduce the scale-dependent Gaussian averaging of a time series as described in this paper: https://arxiv.org/pdf/1706.01126.pdf

The process involves performing a continuous wavelet transform and then using a Gaussian window to estimate a scale-dependent average of the timeseries using the envelope of the wavelet at each scale

The equations from the paper are as follows, and I am looking to implement them in Python:

Wavelet Transform:

$W_i(t, \sigma) = \frac{1}{\sqrt{\sigma}} \int B_i(t')\psi\left(\frac{t' - t}{\sigma}\right) dt'$

Morlet Wavelet: $\psi(\eta) = \pi^{-1/4}e^{-i\omega_0\eta}e^{-\eta^2/2}$

Frequency to Scale Conversion: $f_e = \frac{(\omega_0 + \sqrt{2 + \omega_0^2})}{(4\pi\sigma)}$

Local Scale-Dependent Background Field: $b_i(t, s_b) = \int B_i(t') \exp\left(-\frac{(t' - t)^2}{2\sigma_b^2}\right) dt'$

Here's the code I've written so far using the pycwt library, but I'm unsure if my implementation of the Gaussian window and its convolution with the time series is correct:

import numpy as np
import pycwt as wavelet
from scipy.signal import fftconvolve

def cwt_and_local_average(signal, dt, dj=1/12, s0=-1, J=-1, alpha=1, omega0=6):
    # Perform CWT
    mother = wavelet.Morlet(omega0)
    W, sj, freqs, coi, signal_ft, ftfreqs = wavelet.cwt(signal, dt, dj, s0, J, wavelet=mother)
    
    # Calculate local average using Gaussian windows
    local_averages = np.zeros_like(W)
    for n, scale in enumerate(sj):
        # Convert scale to frequency
        fe = (omega0 + np.sqrt(2 + omega0**2)) / (4 * np.pi * scale)
        # Compute the standard deviation of the Gaussian
        sigma = alpha / fe  # Since alpha is scale factor for the Gaussian width
        # Create Gaussian window
        window_size = int(np.round(3 * sigma / dt))
        t = np.linspace(-window_size, window_size, 2 * window_size + 1)
        gaussian_window = np.exp(-(t * dt)**2 / (2 * sigma**2))
        gaussian_window /= np.sum(gaussian_window)  # Normalize
        
        # Convolve signal with Gaussian window
        local_averages[n, :] = fftconvolve(signal, gaussian_window, mode='same')
    
    # Return everything including local averages
    return W, sj, freqs, coi, signal_ft, ftfreqs, local_averages

# Example usage
dt = 1/1000  # Sampling period
signal = np.random.randn(1000)  # Your signal
W, sj, freqs, coi, signal_ft, ftfreqs, local_averages = cwt_and_local_average(signal, dt)

In particular I am not sure whether the way define the quantity in the exponential is in line with what the paper suggests.

0

There are 0 best solutions below