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

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?

Kind regards,

Laura

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.

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?