Get a PSD per atlas label

Hi all,

I am trieng to get a PSD for each atlas label in source space and wanted to make sure that my approach is correct:

I was recommended to bring the epochs and not the evoked data into source space ( I am using the LCMV beamformer). So I have a source estimate for each epoch. As a next stept I extract the label time courses of each epoch and compute a psd for each epoch for each label. Then I average over the psds within one label, following the suggestion made here.
So now I have a PSD for each atlas label. Does my approach make sense, so first averaging over space (by extracting the label time courses), then computer the psd and then average over the epochs?
Here my code snippet:

labels = mne.read_labels_from_annot(f'sub-{sub:02d}', parc="aparc", subjects_dir=subjects_dir,hemi = "both")
psd_all={}

# create dictionary with labels
for label in range(len(labels)):
    psd_all[labels[label]]=[]

# add epoch psd to each label
for ep in range(len(stc)):
    stc_ep=stc[ep]
    label_ts = mne.extract_label_time_course(stc_ep, labels, src, mode = "mean_flip")
    for label in range(len(labels)):
       ts=label_ts[label,]
      psd, freqs = mne.time_frequency.psd_array_welch(x=ts, sfreq=1000, fmin=0.5, fmax=120, n_fft=10000, n_overlap=5000, average="median", window="hamming")
      psd_all[labels[label]].append(psd)
# average over epochs
for label in range(len(labels)):
    psd_all[labels[label]] = np.average(psd_all[labels[label]], axis = 0)

I am using MNE version 1.2.3

Thank you for your help!

Best,
Antonia

Hi,

Sounds ok to me.

Is there a particular reason/question that you have in mind, that would potentially require another approach?

We ended up discussing this with @antonia at last week’s office hours. Here’s a summary of what I remember (@antonia feel free to chime in to expand on this!)

  1. general approach (extract label time course, then compute PSD on that) is fine
  2. There was some discussion of specific parameters (esp. mode="pca_flip")
  3. computing the PSD on each vertex and then using extract_label_time_course should also work, if there were strong (e.g. hypothesis-driven) reasons to do it that way.
1 Like

Hi,
Sorry for my late response. Yes, those were the three points we discussed. It was recommended that I change the label extraction mode to “pca_flip” since it better represents the dominant power spectra in that label. Also, as I am looking at resting state data, I was not really sure how to quality check my analysis steps since I do not know the ground truth of the data (I cannot check my data by i.e. looking for a motor signal when participants press a button as one would do in task-related data). Hence, it was recommended that for quality checks I should check for alpha-oscillation in parietal and lower occipital regions.
All suggestions were really helpful. Thank you!

1 Like