pythonnumpyscipysignal-processingwavelet

Creating wavelet transform of brain signals using scipy

I am trying to create time-frequency representation of my brain signal. I would like to create the data for frequencies from 0hz - 120hz (so it can cover, delta, theta, alpha, beta, low-gamma and high-gamma frequency bands).

Here is my code:

First, I visualize my brain signal:

ax1.plot(list(range(stc_broca_tmp2.shape[1])), brain_signal, linewidth=3, color='blue')

Results:

enter image description here

brain_signal.shape

I am trying to create time-frequency representation of my brain signal. I would like to create the data for frequencies from 0hz - 120hz (so it can cover, delta, theta, alpha, beta, low-gamma and high-gamma frequency bands).

Here is my code:

First, I visualize my brain signal:

ax1.plot(list(range(stc_broca_tmp2.shape1)), brain_signal, linewidth=3, color='blue') Results:

enter image description here

brain_signal.shape (353191,)

Then I use the following code to compute its time-frequency representation:

t = np.linspace(0, signal_length, signal_length, endpoint=False)
#brain_signal = stc_broca[3, :]
widths = np.arange(0, 130)
cwtmatr = signal.cwt(brain_signal, signal.ricker, widths)
plt.imshow(cwtmatr, extent=[0, signal_length, 0, 130], cmap='PRGn', aspect='auto',
           vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.show()

However, I received the following error:

--------------------------------------------------------------------------- ValueError Traceback (most recent call last) in 2 #brain_signal = stc_broca[3, :] 3 widths = np.arange(0, 130) ----> 4 cwtmatr = signal.cwt(brain_signal, signal.ricker, widths) 5 plt.imshow(cwtmatr, extent=[0, signal_length, 0, 130], cmap='PRGn', aspect='auto', 6 vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())

/usr/local/lib/python3.6/dist-packages/scipy/signal/wavelets.py in cwt(data, wavelet, widths, dtype, **kwargs) 478 N = int(N) 479 wavelet_data = np.conj(wavelet(N, width, **kwargs)[::-1]) --> 480 output[ind] = convolve(data, wavelet_data, mode='same') 481 return output

/usr/local/lib/python3.6/dist-packages/scipy/signal/signaltools.py in convolve(in1, in2, mode, method) 1276 # fastpath to faster numpy.convolve for 1d inputs when possible 1277 if _np_conv_ok(volume, kernel, mode): -> 1278 return np.convolve(volume, kernel, mode) 1279 1280 return correlate(volume, _reverse_and_conj(kernel), mode, 'direct')

<array_function internals> in convolve(*args, **kwargs)

~/.local/lib/python3.6/site-packages/numpy/core/numeric.py in convolve(a, v, mode) 813 raise ValueError('a cannot be empty') 814 if len(v) == 0: --> 815 raise ValueError('v cannot be empty') 816 mode = _mode_from_name(mode) 817 return multiarray.correlate(a, v[::-1], mode)

ValueError: v cannot be empty

My brain signal is not empty. Why do I receive such an error?

Thanks in advance

Solution

  • Your widths has to start from 1 not 0. Here's an illustration using random data.

    import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal


N = 5000
rnd = np.random.RandomState(12345)

brain_signal = np.sin(np.linspace(0, 1000, N)) + rnd.uniform(0, 1, N)
widths = np.arange(1, N//8)
cwtmatr = signal.cwt(brain_signal, signal.ricker, widths)

fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(10, 6))
axes = ax.flatten()

sns.lineplot(np.linspace(0, 1000, N), brain_signal, ax=axes[0], lw=2)
sns.heatmap(cwtmatr, cmap='Spectral', ax=axes[1]);

axes[0].set_title('Brain signal')
axes[1].set_title('CWT of brain signal')

plt.tight_layout()

    CWT brain signal