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

def calculate_local_mean(x, dt, dj, alpha=1, omega0=5.5):
    # Initialize the Morlet wavelet with the desired parameter
    mother_wavelet = wavelet.Morlet(w=omega0)
    # Perform the continuous wavelet transform using pycwt
    wavelet_coeffs, scales, freqs, coi, _, _ = wavelet.cwt(x, dt, dj, wavelet=mother_wavelet)
    # Initialize the array to store local means
    local_means = np.zeros_like(wavelet_coeffs)
    # Iterate over each scale and compute the Gaussian window and then the local mean
    for i, scale in enumerate(scales):
        # The frequency associated with this scale
        fe = freqs[i]
        # The standard deviation of the Gaussian in terms of time
        sigma_t = alpha / (2 * np.pi * fe)
        # Create the Gaussian window, ensuring it spans enough points in time
        window_size = int(np.ceil(3 * sigma_t / dt))
        t = np.arange(-window_size, window_size + 1) * dt
        gaussian_window = np.exp(-0.5 * (t / sigma_t) ** 2)
        gaussian_window /= gaussian_window.sum()
        # Convolve the input signal with the Gaussian window
        local_means[i, :] = np.convolve(x, gaussian_window, mode='same')
    return scales, freqs, local_means

# Example usage
dt = 1 / 1000  # Sampling interval in seconds
dj = 0.1       # Scale resolution
x = np.random.randn(1000)  # Replace with your actual time series data

scales, frequencies, local_means = calculate_local_mean(x, dt, dj)

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



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.