CODEWITHSUNDEEP
pythonnumpyscipysignal-processingwavelet-transform
Kadaj13
Kadaj13

Reputation: 1551

Creating time-frequency representation of brain-waves using scipy

I want to create time-frequency representation of brain waves (neural activity data obtained by MEG neuroimaging).

My question is followed up by my previous question

What I need, is the time-frequency representation (amplitude values or power values) for different frequency bands of brain waves.

More specifically, frequency band of brain waves are defined as below:

delta: 0.5-4 hz theta: 4-8 hz alpha: 8-12 hz beta: 12-30 hz low_gamma 30-80 hz high_gamma = 80 - 120 hz

My (fixed) code from the previous thread:

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

My question is, do the values on y-axis correctly show frequency information? If no, how should I fix it?

Update: Results of the above code:

enter image description here

Thanks in advance

Upvotes: 0

Views: 581

Answers (1)

cmbfast
cmbfast

Reputation: 489

Like we discussed on comments to your older question, the Y-axis of CWT is just a scaled version of the input wave based on a certain wavelet, where as X-axis is the time point on the original wave. Wavelet transform is just a convolution of wavelet and the original wave. As mentioned in comments by Tim, you might probably be looking for FFT. Here's an example:

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq

# Number of data points    
N = 5000
rnd = np.random.RandomState(12345)
# Sampling interval
T = 1/1000

I don't have your data, so I just used your frequencies above for different brain signal to make a fictional brain signal. I also added a bit of noise to make it random and more realistic.

## Frequencies
delta = 2
theta = 5
alpha = 10
beta = 20
low_gamma = 75
high_gamma = 100


x = np.linspace(0, N*T, N)
brain_signal = np.sin(delta * 2.0*np.pi*x) + \
               np.sin(theta * 2.0*np.pi*x) + \
               4 *np.sin(alpha * 2.0*np.pi*x) + \
               np.sin(beta * 2.0*np.pi*x) + \
               2 * np.sin(low_gamma * 2.0*np.pi*x) + \
               10 * np.sin(high_gamma * 2.0*np.pi*x) + rnd.uniform(-10, 10, N)

Notice that in this brain_signal, high_gamma has the highest amplitude (10). Now let's see if FFT can recover any of these different type of waves form one signal.

yf = fft(brain_signal)
xf = fftfreq(N, T)[:N//2]

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

sns.lineplot(x, brain_signal, ax=axes[0], lw=2)
sns.lineplot(xf[:1000], 2.0/N * np.abs(yf[0:N//2])[:1000], ax=axes[1]);


axes[0].set_title('Brain signal \nTime domain')
axes[0].set_xlabel('Time')
axes[0].set_ylabel('Amplitude')

axes[1].set_title('FFT of brain signal \nFrequency domain')
axes[1].set_xlabel('Frequency')
axes[1].set_ylabel('Amplitude')

plt.tight_layout()

FFT brain waves Indeed, it did. See the different peaks, and as expected, since high_gamma was the strongest in our signal, we get a sharp peak at 100Hz. We also see different peaks for different type of brain waves. However, note that your real world data may not give a clear FFT spectrum as this one because of noise. You can test that by making this data more noisy.

Upvotes: 2

Related Questions