Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correct factors for complex Morlet wavelet #740

Open
itamwa opened this issue May 9, 2024 · 1 comment
Open

Correct factors for complex Morlet wavelet #740

itamwa opened this issue May 9, 2024 · 1 comment
Labels

Comments

@itamwa
Copy link

itamwa commented May 9, 2024

Hi there,

I'm using the continuous complex Morlet transform to compare a local energy spectrum and qualitatively I get what I expect (from Fourier transform, for example). Quantitatively, however, I'm pretty off.

I suspect two things:

  1. In the literature I see that for the wavelet I want to use should be in this form - psi1(x)= pi^(-0.25) * e^(iwx) * e^(-x^2/2), but the cwt function only allows me to use psi2(x)= (2pi)^(-0.5) * e^(iw*x) * e^(-x^2/2). Is there a way for me to conduct the transform with psi1(x) rather than psi2(x)
  2. In order to compare the amount of energy correctly, do I need to normalize my result by the scale? I believe in the literature it's refereed to as L1 normalization.

Thank you very much,
Itamar

@eulerleibniz
Copy link

As far as I know, for a decent visualization, you need a scaling like this code, which still doesn't exactly provide the normalized values but will almost perfectly match the visualization as in MATLAB's cwt function (You will need a constant scaling on top of this I believe). You might also be interested in this https://github.com/neergaard/CWT.git which only works for cmore6.0-0.5, but comes with perfectly correct normalization out of the box :

# IMPORTS
import numpy as np
import matplotlib.pyplot as plt
# import ptwt # FOR GPU CASE WHICH I NEEDED <3
import pywt


# Parameters for the signal generation
WINDOW = 30  # window in seconds
Fs = 512  # sampling rate in Hz, I think pywt has bad implementation and usually requires higher Fs to perform well
# Three distinct frequencies for generating a nice signal, plot the results in Log space ideally
F1 = 10
F2 = 1
F3 = 100
nData = WINDOW * Fs
t = np.linspace(0, WINDOW, nData, endpoint=False)
T1 = 0.33
T2 = 0.66

# Generate the signal with 3 different sin with different time localizations and amplitudes currently set to 1
sig = (1 * np.sin(2 * np.pi * F1 * t) * (t < T1 * WINDOW) +
       1 * np.sin(2 * np.pi * F2 * t) * ((t > T1 * WINDOW) & (t < T2 * WINDOW)) +
       1 * np.sin(2 * np.pi * F3 * t) * (t > T2 * WINDOW))

# sig += 0.1 * np.random.randn(len(sig))
plt.plot(sig)

wavelet = "cmor6.0-0.5"
nv = 10 # Number of Voices
num_scale = int(nv * np.log2(nData // 2))
min_freq = 1 / WINDOW
max_freq = Fs / 2
# Choose the scales wisely. i think the following is a decent one, logarithmic frequencies 
wavescales = pywt.frequency2scale(wavelet, np.geomspace(min_freq, max_freq, num_scale, endpoint=False) / Fs)

coefficients, freqs = pywt.cwt(
    sig, wavescales,  wavelet, sampling_period=1/Fs, method='fft')
# Different visualizations exist .real .imag .... for visualizing phase and so on
coefficients = np.abs(coefficients)

# THE SCALING REQUIRED to make Sin(2 * pi * m * x) and Sin(2 * pi * n * x) appear ALMOST the same
# NOT EXACTLY SCALED
coefficients = coefficients * np.sqrt(freqs[:, np.newaxis])

plt.figure(figsize=(20, 10))
# Use seismic colormap for visualizing the case with negative values, otherwise the default is fine
print(np.max(coefficients))
print(np.min(coefficients))
plt.pcolormesh(t, freqs, coefficients, shading="auto", cmap="seismic") 
plt.yscale("log")
plt.show()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants