Is there a similar function to the spectopo function from MATLAB?

💬 Support & Discussions
,
1

Hello everyone,

I am trying to calculate the powers of all frequencies (delta, theta, alpha, sigma, beta) during a 3 min. EEG recording with 16 (biosemi16 montage) channels. I want to calculate the average powers of each frequency over the whole recording. MATLAB has a straight-forward function for this:

% EEGLAB spectopo function
[spectra,freqs] = spectopo(EEG.data(2,:), 0, EEG.srate, 'nfft', 1024); %
For channel 2 (C4)

thetaIdx   = find(freqs>=4 & freqs<=8);                    % theta=4-8
thetaPower = 10^(mean(spectra(thetaIdx))/10);              % mean theta
power
alphaIdx   = find(freqs>=8 & freqs<=12);                   % alpha=8-12
alphaPower = 10^(mean(spectra(alphaIdx))/10);              % mean alpha
power
eeglabThetaAlphaRatio = thetaPower/alphaPower;

figure; subplot(1,2,1); plot(freqs,10.^(spectra/10))

% Matlab pwelch function
x  = EEG.data(2, :);      % For channel 2 (C4)
Fs = EEG.srate;              % Sampling frequency
L  = EEG.pnts;               % Length of signal
NFFT = 2^nextpow2(L);        % Next power of 2 from length of x
% X    = fft(x,NFFT)/L;
% freqs   = Fs/2*linspace(0,1,NFFT/2+1);
% spectra = abs(X(1:NFFT/2+1));

NOVERLAP = 0;
WINDOW = 512;
[spectra,freqs] = pwelch(x,WINDOW,NOVERLAP,NFFT,EEG.srate);

thetaIdx   = find(freqs>=4 & freqs<=8);            % theta=4-8
thetaPower = mean(spectra(thetaIdx));              % mean theta power
alphaIdx   = find(freqs>=8 & freqs<=12);           % alpha=8-12
alphaPower = mean(spectra(alphaIdx));              % mean alpha power
matlabThetaAlphaRatio = thetaPower/alphaPower;

But I want to achieve a similar output but then with MNE. I tried this function to do so:

def eeg_power_band(epochs):

    # specific frequency bands
    FREQ_BANDS = {"delta": [0.5, 4.5],
                  "theta": [4.5, 8.5],
                  "alpha": [8.5, 11.5],
                  "sigma": [11.5, 15.5],
                  "beta": [15.5, 30]}

    psds, freqs = psd_welch(epochs, picks='eeg', fmin=0.5, fmax=30.)
    # Normalize the PSDs
    psds /= np.sum(psds, axis=-1, keepdims=True)

    X = []
    for fmin, fmax in FREQ_BANDS.values():
        psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
        X.append(psds_band.reshape(len(psds), -1))

    return np.concatenate(X, axis=1)

But it’s not working because I am not working with epochs and I want to calculate the power of all frequencies over all channels and not just one signal. Is that possible to do in MNE?

Thank you for your help!

Kind regards,

Laura

2

Hi @lrajos,
could you be a bit more specific in what you mean by:

But it’s not working because I am not working with epochs (…)

mne.time_frequency.psd_welch can be used on both Epochs and Raw objects.
You also don’t have to average frequency ranges if you want to see power for each frequency - it is stored in the psds variable. For Raw data psds will be a n_channels x n_frequencies array.

3

Hello Mikolaj! Thank you a lot for your response. When I run this function on my raw data, the values seem very odd. They are extremely low and don’t seem to respresent the correct powers of each frequency. I stated as similar question here: Calculate spectral density including the values. Can you see the problem?