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

Question

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.