Extracting different frequency bands (alpha, beta, gamma) from MEG source estimated data

Question

I have a dataset that consists of source estimated data of MEG brain signals.

I need to extract the different frequency band features (i.e., alpha, beta, gamma, etc) from them.

What steps/procedure should I take?

My data consists of 1D numpy array for each brain region.

From a little search, it seems I have to take fourier transform first. But I couldn't understand why I should use fourier transform and what to do after that.

Information about source estimated data can be found here.

I want to use python to calculate them.

Thanks in advance...

Answer 1

MEG data is non-stationary and FFT can provide a crude approximation on clean data at best. Instead I'd recommend time-frequency analysis, like STFT or CWT.

CWT is favored as it adjusts its resolution to better discriminate higher frequencies in time, and lower frequencies in frequency, on logscale (which is appropriate for brain waves). Afterwards one can apply ridge extraction to extract independent components as 1D frequency vs time arrays. CWT and STFT can both be enhanced by synchrosqueezing, which studied have applied to EEG.

Examples

Ridge extraction example

EEG example:

References

1. ssqueezepy - CWT, STFT, synchrosqueezing
2. EEG synchrosqueezing
3. I'd recommend learning some SP basics, as it'll save much time and hassle in the long run; I've made a list here.

Alternative

FFT offers a simpler but limited approach: algorithm below will estimate M most dominant frequencies. Improve accuracy by splitting up input and feeding separately.

import numpy as np
def est_freqs(x, M=3, T=1, smoothing=None):
    xf = np.abs(np.fft.rfft(x))
# to avoid picking from around the same bin
if smoothing not in (0, 1):
    if smoothing is None:
        smoothing = max(len(x) // 500, 2)
    # ensure we can reshape
    xf = np.pad(xf, [0, smoothing - len(xf) % smoothing])
    xfs = xf.reshape(-1, smoothing).mean(axis=1)
else:
    xfs = xf

# sort descending
idxs = list(range(len(xfs)))
xf_idxs_sort = sorted(zip(list(xfs), idxs), key=lambda x: x[0], reverse=True)
idxs_sort = [a[1] for a in xf_idxs_sort]    
# map back to un-smoothed and account for `T`
idxs_sort = np.round(len(x) / len(xfs) / 2 * np.array(idxs_sort) / T
                     ).astype(int)

return idxs_sort[:M]


fs = 400         # sampling frequency [Hz]
T  = 3           # duration [sec]
noise_std = 1    # WGN standard deviation
freqs = [4, 28]  # frequencies to simulate
amps  = [3,  2]  # their amplitudes
N = int(fs * T)
t = np.linspace(0, T, N, endpoint=False)
x = np.zeros(N)
for f, a in zip(freqs, amps):
    x += a * np.cos(2np.pi  f * t)
x += noise_std * np.random.randn(N)
print(est_freqs(x, T=T, M=len(freqs)))

>>> [4, 28]