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.