PyWavelets / pywt

PyWavelets - Wavelet Transforms in Python
http://pywavelets.readthedocs.org
MIT License
1.97k stars 460 forks source link

Ripples in output although low wavelet center frequency and high scale #727

Open faering opened 3 months ago

faering commented 3 months ago

Hi PyWavelets developers and community,

I am using a complex morlet wavelet to compute the time-frequency response and there are some ripples in the output. I am hoping someone can help me understand why this is occurring.

Figure 1) Signal and FFT To demonstrate I've created a chirp signal: image

Figure 2) Time-frequency Below you see the time-frequency of the wavelet transform and a zoom to show the ripples in question. image

This occurs even if I change the interpolation to 'none' in matplotlib's imshow function (I am using 'hanning' in the above plots).

Figure 3) CWT output for 2 Hz Below you see a plot of the output for one frequency (2 Hz) from the continuous wavelet transform with a scale of 256 for the wavelet: image

Question(s)

  1. My question is how can the output of the wavelet transform oscillate as shown in Figure 3, when the scale is 256 (consider time 0.5 s to 1.0 s)?
  2. Wouldn't a large scale i.e. 256 mean that the wavelet is overlapping with almost 0.5 seconds (fs = 512 Hz) of the data and thus not be sensitive to sudden phase changes?

Code to Reproduce


from scipy.signal import chirp
import pywt
import matplotlib.pyplot as plt

# parameters
fs                          = 512 # Hz
dt                          = 1/fs
time                        = np.arange(-1, 1+1/fs, 1/fs)
frequencies                 = np.arange(2, 32, 2)
bandwidth                   = 1.5
center_frequency            = 2.0
wavelet                     = pywt.ContinuousWavelet(f'cmor{bandwidth}-{center_frequency}') # complex morlet
cmap                        = 'RdBu_r'

# create signal
signal = chirp(t=time, f0=6, f1=1, t1=0.5, method='linear')

# fft
y = fft(signal)
xf = fftfreq(signal.size, 1/fs)
xf = xf[xf >= 0]
y = np.abs(np.ravel(y))[:len(xf)]

# plot chirp
fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(15, 3))
ax = ax.flatten()
ax[0].plot(time, signal)
ax[0].set_title('Chirp Signal')
ax[0].set_xlabel('Time (s)')
ax[0].set_ylabel('Amplitude')
ax[1].plot(xf, y)
ax[1].set_title('Mangitude Spectrum')
ax[1].set_xlabel('Frequency (Hz)')
ax[1].set_ylabel('Magnitude')
plt.show()

# scales
frequencies = np.array(frequencies) / fs  # normalise frequencies
scales = pywt.scale2frequency(wavelet, frequencies)

# compute time-frequency representation
[coeffs, freqs] = pywt.cwt(data=signal, scales=scales, wavelet=wavelet, sampling_period=dt)
magnitude = abs(coeffs)

# time-frequency plot
fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(15, 5))
ax = ax.flatten()
im = ax[0].imshow(
    magnitude, 
    aspect='auto',
    extent=[time[0], time[-1], (freqs)[0], (freqs)[-1]],
    origin='lower', 
    cmap=cmap,
    interpolation='hanning')
ax[0].set_title('Time-frequency Amplitude', fontsize=16)
ax[0].set_ylabel('Frequency (Hz)', fontsize=10)
ax[0].set_xlabel('Time (s)', fontsize=10)
ax[0].set_yscale('log')
ax[0].set_yticks([2, 10, 20, 30], [2, 10, 20, 30])
# zoom
im = ax[1].imshow(
    magnitude, 
    aspect='auto',
    extent=[time[0], time[-1], (freqs)[0], (freqs)[-1]],
    origin='lower', 
    cmap=cmap,
    interpolation='hanning')
ax[1].set_title('zoom', fontsize=16)
ax[1].set_xlabel('Time (s)', fontsize=10)
ax[1].set_yscale('log')
ax[1].set_yticks([2, 10, 20, 30], [2, 10, 20, 30])
cbar = fig.colorbar(im, orientation="vertical")
cbar.ax.set_ylabel("Magnitude", rotation=270, labelpad=15, fontsize=12)
ax[1].set_xlim(0.5, 1)
ax[1].set_ylim(2, 3)
plt.show()

# plot wavelet transform for f=2 Hz
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(15, 3))
ax.plot(time, magnitude[0])
ax.set_title('CWT Power for 2 Hz')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Magnitude')
plt.show()```