Scaling mismatch between PyWavelets CWT and my implementation

[grigorishat]

I would like to understand the scaling in the CWT implementation in Python. I know that there are topics about it in Stackoverflow, although I cannot understand why my implementation does not match with the implementation in PyWavelets.

My approach is just an convolution between my signal and the wavelet. I know the implementation here is based on an old version of a matlab implementation described here.

I’ve attached an example script and resulting plots at the end. The code snippet below shows that when I do the multiplication by 1 / np.sqrt(fs) (fs being the sampling rate) , the amplitude of my implementation overlays compared to PyWavelets. Without that factor, the signals are off.

This extra factor seems unintuitive, and I’m concerned I’m missing a detail in either the wavelet definition or the normalisation.

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
import pywt

# create a test signal containing 40Hz and 10Hz 

# --- Parameters ---
fs = 3000       # Sampling rate in Hz
duration = 10   # Total duration in seconds
N = int(fs * duration)
t = np.linspace(0, duration, N, endpoint=False)

# --- 10 Hz signal (continuous) ---
f1 = 10
sig10 = np.sin(2 * np.pi * f1 * t)

# --- Create raw on/off mask for 40 Hz signal ---
# intervals: [0–1s], [2–4s], [5–8s]
mask_raw = ((t >= 0) & (t < 1)) | ((t >= 2) & (t < 4)) | ((t >= 5) & (t < 8))
mask_raw = mask_raw.astype(float)  # convert True/False to 1.0/0.0

# --- Smooth the mask with a large Gaussian filter to get ~0.8s transitions ---
sigma = 667  # ~ (0.8s * fs) / 6
mask_smooth = gaussian_filter1d(mask_raw, sigma=sigma)

# --- Create 40 Hz signal and multiply by the smoothed mask ---
f2 = 40
sig40 = np.sin(2 * np.pi * f2 * t) * mask_smooth

# --- Combine signals ---
signal = sig10 + sig40

# --- Plot the final combined signal ---
plt.figure()
plt.plot(t, signal)
plt.title("Test Signal (10 Hz + Intermittent 40 Hz)")
plt.xlabel("Time in s")
plt.ylabel("Amplitude")
plt.xlim(0, 10)

B = 14
C = 2

wavelet = f'cmor{B}-{C}'

# I aim to analyse 40 Hz so I calculate the scale
# freq_analyse is the frequency to analyse diviced by the sampling rate fs
freq_analyse = 40 / fs

scale = pywt.frequency2scale(wavelet, freq_analyse)

# calculate coefficients
coefs, freqs = pywt.cwt(signal, scale, wavelet, sampling_period=1/fs)

# now my implementation follows

a=0 # no shifting, this is done by the convolution below

b = scale/fs # b is the scale in the wavelet function, it can be also computed as: b = C/f with f =40Hz the frequency of interest

# create the wavelet with the same "sampling rate" as my test signal above

# wavelet should go from -2 to 2
wl_start = -2
wl_end = 2
t_wl = np.arange(wl_start, wl_end+1/fs,1/fs)

psi_morlet = 1/np.sqrt(B*np.pi) * np.exp(-((t_wl-a)/b)**2/B) * np.exp(1j*2*np.pi*C*((t_wl-a)/b))

# calculate the convolution 
conv_result =  np.convolve(np.conj(psi_morlet), signal, mode='full')

# scale the result by 1/sqrt(b)

conv_result = 1/np.sqrt(b) * conv_result

# plot the coefficients of the cwt over the time

fig, ax = plt.subplots(1, 1, figsize=(10, 10))

ax.plot(t, abs(coefs[0]), label='pywt implementation')

# create time vector for the convolution
t_conv = np.arange(wl_start,duration+wl_end+1/fs, 1/fs)
ax.plot(t_conv,abs(conv_result), label='my implementation') 

ax.legend()

ax.set_xlabel('Time in s')
ax.set_ylabel('CWT')


# plot the coefficients of the cwt over the time

fig, ax = plt.subplots(1, 1, figsize=(10, 10))

ax.plot(t, abs(coefs[0]), label='pywt implementation')

# create time vector for the convolution
t_conv = np.arange(wl_start,duration+wl_end+1/fs, 1/fs)
# IF IS SCALE MY IMPLEMENTATION BY 1/sqrt(fs) THAN THE RESULTS MATCH
ax.plot(t_conv,abs(conv_result)*1/np.sqrt(fs), label='my implementation') 

ax.legend()
ax.set_xlabel('Time in s')
ax.set_ylabel('CWT')