Correct factors for complex Morlet wavelet #740

Activity

1. agriyakhetarpal

   added

   question

   on May 9, 2024
2. eulerleibniz commented on Nov 3, 2024

   

   eulerleibniz

   on Nov 3, 2024

   More actions

   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()