Ripples in output although low wavelet center frequency and high scale
faering opened this issue · 0 comments
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:
Figure 2) Time-frequency
Below you see the time-frequency of the wavelet transform and a zoom to show the ripples in question.
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:
Question(s)
- 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)?
- 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()```