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

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.