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
- ssqueezepy - CWT, STFT, synchrosqueezing
- EEG synchrosqueezing
- 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(2*np.pi * f * t)
x += noise_std * np.random.randn(N)
print(est_freqs(x, T=T, M=len(freqs)))
>>> [4, 28]
