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:
- 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
- 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)
- 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?