Power spectral density of EEG signal -- code sanity check

💬 Support & Discussions
1

MNE-Python version: 0.22.1
OS: Windows 10

Hi MNE folks,

I’m working on incorporating PSD calculations into my analysis pipeline, and I wanted to make sure that I’m using psd_welch function correctly. As the only one on my team working on this, a simple “looks good” would go a long way! I’ve used the ‘Sleep stage classification from polysomnography (PSG) data’ example from the docs (Sleep stage classification from polysomnography (PSG) data — MNE 0.23.4 documentation) as a guide. This is what I have to make my calculation:

freq_bands = {'alpha':[8.5,11.5],'theta':[4.5,8.5]}
psds,freqs = mne.time_frequency.psd_welch(data,picks='eeg',fmin=0.5,fmax=30.,tmin=45,tmax=345)
psds /= np.sum(psds,axis=-1,keepdims=True) # Normalize PSDs
bandpower = []
for fmin,fmax in freq_bands.values(): # Calculate bandpower for bands of interest
    psds_band = psds[:,(freqs>=fmin) & (freqs < fmax)].mean(axis=-1)
    bandpower.append(psds_band.reshape(len(psds),-1))

# Get global means
bandpower_global_means = [np.mean(band) for band in bandpower] # average across electrodes for global mean
bandpower_global_means.append(float(bandpower_global_means[0]) / float(bandpower_global_means[1])) # theta/alpha ratio
# Get frontal means
bandpower_frontal_means = [np.mean(band[3:9]) for band in bandpower]
bandpower_frontal_means.append(float(bandpower_frontal_means[0]) / float(bandpower_frontal_means[1]))

The goal is to get global (across all electrodes) PSD in alpha, theta, and alpha/theta ratio as well as the same metrics for the most frontal electrodes (the third through eighth channels in this case). Just wanted to make sure that my pipeline is doing what I think it’s doing!

As a final question, do I need to do a separate psd_welch on the subset of frontal electrodes, or is the slice that I have on the second-to-last line sufficient?

Thank you for your help!

Best,
Neurovella