Computation of psd on dipole data

💬 Support & Discussions
1

Hi everyone,

I am running source reconstruction with MNE, and am running into some issues when trying to use compute_psd on the source reconstructed data (where channel type is dipole).

This is my setup:

  • MNE version: 1.5.0
  • operating system: Manjaro Linux

I have the following workflow:

  1. Create a forward solution using the fsaverage template. I load the Desikan-Killiany atlas parcelation:
trans = "fsaverage"
src = os.path.join(fs_dir, trans, "bem", "fsaverage-ico-5-src.fif")
bem = os.path.join(fs_dir, trans, "bem", "fsaverage-5120-5120-5120-bem-sol.fif")

sourcespace = mne.read_source_spaces(src, verbose=0)

# make forward solution based on the first file
# this assumes that all files have the same channel layout!!
signal_container = SignalContainer(to_source_reconstruct_files[0])
signal_container.load()
fwd = mne.make_forward_solution(signal_container.raw.info, trans=trans, src=src,
                                    bem=bem, eeg=True, mindist=5.0)
mne.write_forward_solution(fname_fwd, fwd, overwrite=True)

# Load parcellations
# Parcellation
# Desikan-Killiany atlas (34 ROI from both hemispheres = 68 ROIs)
# Named aparc.annot in MNE python fsaverage folder
labels = mne.read_labels_from_annot("fsaverage", parc="aparc",
                                    subjects_dir=fs_dir, verbose=0)
labels = labels[:-1]  # remove 2 "unknown" vertices
  1. Given a loaded EEG file, run source reconstruction, and create a new file were the signals corresponding to the 68 DK-atlas parcels are stored as dipoles. I exlicilty set the unit to Ampere.
n_ch, n_timepoints = raw.get_data().shape
# Use different forward solutions depending on number of channels

# Using assumption about equal variance and no correlations I make a diagonal matrix
# Using the default option for 0.2µV std for EEG data
noise_cov = mne.make_ad_hoc_cov(raw.info, None)

# Make inverse operator
# Using default depth parameter = 0.8 and free orientation (loose = 1)
inverse_operator = mne.minimum_norm.make_inverse_operator(raw.info,
                                                        fwd, noise_cov,
                                                        loose = 1, depth = 0.8,
                                                        verbose = 0)
src_inv = inverse_operator["src"]
# Compute inverse solution and retrieve time series for each label
# Preallocate memory
label_ts = np.full((len(labels),n_timepoints),np.nan)
# Define regularization
snr = 3 # Default setting
# A for loop is used for each label due to memory issues when doing all labels at the same time
for l in range(len(labels)):
    stc = mne.minimum_norm.apply_inverse_raw(raw,inverse_operator,
                                                lambda2 = 1/(snr**2),
                                                label = labels[l],
                                                pick_ori = "vector",
                                                method = "MNE", verbose = 0)
    # Use PCA to reduce the 3 orthogonal directions to 1 principal direction with max power
    # There can be ambiguity about the orientation, thus the one that
    # is pointing most "normal", i.e. closest to 90 degrees to the skull is used

    stc_pca, pca_dir = stc.project(directions="pca", src=src_inv)
    # Get mean time series for the whole label
    temp_label_ts = mne.extract_label_time_course(stc_pca, labels[l], src_inv, mode="mean_flip",
                                    return_generator=False, verbose=0)
    # Save to array
    label_ts[l,:] = np.squeeze(np.array(temp_label_ts))
    del stc # Free up memory
    # print("Finished estimating STC for {} out of {} ROIs".format(l+1,len(labels)))

# Prepare variables
sfreq=raw.info["sfreq"]

source_info = mne.create_info(label_names, sfreq, ch_types="dipole")
#set units to ampere
for s in range(num_sources):
    source_info['chs'][s]['unit']= mne.io.constants.FIFF.FIFF_UNIT_A
new_raw = mne.io.RawArray(label_ts,source_info)
  1. I compute the power spectrum for the source reconstructed signals, :
spectrum = new_raw.compute_psd(method='welch', fmin=l_freq_full, fmax=h_freq_full, n_fft=n_fft, n_overlap=n_overlap, n_per_seg=n_fft, picks='all')
    freqs = spectrum.freqs
    psds = spectrum.get_data(picks='all')

However, I get the following error:

ValueError: picks (NoneNone, treated as "data_or_ica") yielded no channels, consider passing picks explicitly

It looks like the dipole type channels are not considered data channels, and it is not possible to run the compute_psd function on them.

Is there a reason why it is not possible to run compute_psd on dipole data ? Is there an appropriate solution for this (other than inadvertently setting the source channels as eeg type instead of dipole type) ?

I have found another function: https://mne.tools/stable/generated/mne.minimum_norm.compute_source_psd.html , which seems to be intended for source reconstructed data. What is the difference between running this function vs running compute_psd on the source reconstructed data?