I'm trying to use scipy.signal.cwt()
to do a time-frequency decomposition of neural data. However, I don't completely understand the "widths" parameter.
I already looked at this answer. On the one hand, I had a hard time understanding the question as it was written. On the other, the answers seem quite outdated (for instance, scipy.signal.cwt()
now allows complex wavelets. I also prefer to stick to Morlet wavelets than ricker
, as I understand they are more standard in the literature looking at neural data).
According to the SciPy documentation:
The second is a width parameter, defining the size of the wavelet (e.g. standard deviation of a gaussian).
The book Analyzing Neural Time Series Data by Mike X Cohen seems to define it as follows (for complex Morlet wavelets):
The number of cycles of the Gaussian taper defines its width, which in turn defines the width of the wavelet.
Unfortunately I'm still unclear on how to translate this into code, especially as they seem to be referring to different quantities.
Should the width be a vector of integers defining the number of cycles for the wavelet? Or the standard deviation of the Gaussian used to convolve the sine wave (which produces the Morlet wavelet)?
Relevant code: (Right now I'm literally passing in the frequencies I want to extract, not any "width" derivation of them.)
### Wavelet parameters
fs = 1000 # sampling frequency
min_freq = 1
max_freq = int(fs / 4) # 4 points per cycle
num_freqs = int(max_freq / 2)
freq = np.logspace(np.log10(min_freq), np.log10(max_freq), num_freqs)
sec = 5
n_samples = sec*fs
brain_signal = samples[10*n_samples:11*n_samples] # extracting just 5s from the long neural signal
cwtmatr = signal.cwt(brain_signal, signal.morlet2, freq)
cwtmatr = np.abs(cwtmatr)**2
# Extrapolating the time axis
t = np.linspace(0, sec, num=num_samples)
plt.pcolormesh(t, freq, cwtmatr, cmap='Spectral', shading='gouraud')
plt.colorbar()
I realize my y-axis is wrong, but I would save this for a separate question. Step by step. ;)