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

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

Open
faering opened this issue Mar 14, 2024 · 0 comments
Open
Labels

Comments

@faering
Copy link

faering commented Mar 14, 2024

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

2 participants