How to segment frequency bands in the source reconstruction process using LCMV beamformer?

💬 Support & Discussions
1
  • MNE version:1.3.0
  • operating system: Ubuntu 18.04
    Hi all,after studying the content on the page Source reconstruction using an LCMV beamformer — MNE 1.4.0 documentation, I have rewritten a similar code using the LCMV method to perform source reconstruction and obtain stc data. My experimental objective is to investigate the differences in the source reconstruction results at different frequency bands. Therefore, I tried to segment the original data into frequency bands, but I am not sure if I did it correctly. Here is my code:
data_path = sample.data_path()
subjects_dir = data_path / "subjects"
meg_path = data_path / "MEG" / "sample"
raw_fname = meg_path / "sample_audvis_filt-0-40_raw.fif"

# Read the raw data
raw = mne.io.read_raw_fif(raw_fname)
raw.info["bads"] = ["MEG 2443"]  # bad MEG channel

# Set up the epoching
event_id = 1  # those are the trials with left-ear auditory stimuli
tmin, tmax = -0.2, 0.5
events = mne.find_events(raw)

# pick relevant channels
raw.pick(["meg", "eog"])  # pick channels of interest

# Create epochs
proj = False  # already applied
epochs = mne.Epochs(raw,events,event_id,
    tmin,tmax,baseline=(None, 0),
    preload=True,proj=proj,
    reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))

fwd_fname = meg_path / "sample_audvis-meg-oct-6-fwd.fif"
fwd = mne.read_forward_solution(fwd_fname)

src_fs_sur = mne.read_source_spaces('/usr/local/freesurfer/subjects/fsaverage/bem/fsaverage-ico-5-src.fif')
morph = mne.compute_source_morph(
    fwd['src'], subject_from='sample', subjects_dir='/home/user/mne_data/MNE-sample-data/subjects',
    src_to=src_fs_sur, verbose=True)
    
subfre=[[1,4],[4,8],[8,15],[15,35]]
subfrename=['delta','theta','alpha','beta']
for fre_i, frename in zip(subfre, subfrename):
    print('-------in freband:', frename)
    epochs_filter = epochs.copy().filter(fre_i[0], fre_i[1])
    evoked = epochs_filter.average()
    data_cov = mne.compute_covariance(epochs_filter, tmin=-0.2, tmax=0.5, method="auto")
    noise_cov = mne.compute_covariance(epochs_filter, tmin=-0.2, tmax=0, method="auto")
    rank = mne.compute_rank(epochs, tol=1e-6, tol_kind='relative') 
    filters_lcmv = make_lcmv(
        epochs.info,
        fwd,
        data_cov,
        reg=0.05,
        noise_cov=noise_cov,
        pick_ori="max-power",
        weight_norm="unit-noise-gain",
        rank=None)
    lcmv_stc = apply_lcmv(evoked, filters_lcmv)
    lcmv_stc.apply_baseline((-0.5, 0))
    lcmv_stc_fsaverage = morph.apply(lcmv_stc)       
print('process completed!')

I used the following code to visualize the results again, which seems a bit out of place. (The following shows the results of the beta frequency band [15Hz to 35Hz])

brain = lcmv_stc_fsaverage.plot(
    subjects_dir='/home/user/mne_data/MNE-sample-data/subjects',
    hemi='both',
    views='dorsal',
    initial_time=0.1,
    brain_kwargs=dict(show=False))

image
image511×681 80.1 KB

The main problem is that the STC obtained has strange fluctuations, which are not consistent with the results given in the tutorial.
I also used the following code to view the evoked data obtained using my method and found that the evoked data also has the same strange fluctuations.

evoked.plot_joint()

image

I tested the same code on the somato dataset provided by mne again, and the results were even more strange. Below are the graphs of stc and the corresponding evoked data.(also in beta)

image
image519×693 89.2 KB

image
It seems that the activation regions displayed using my method are almost correct, but the source time course will fluctuate greatly, and the results obtained without dividing frequency bands will not have such large fluctuations. What is the problem? Did I make a mistake in dividing the frequency band like this? How to correct it?
I hope someone can help me, I have been troubled for a very long time, thank you very much！

2

@jshlyz - everything seems to be fine, except I assume you want an average envelope of the beta frequency versus the average of the beta-filtered timeseries.

You can calculate the hilbert transform of your epoch data epo.apply_hilbert(envelope=False). Then project your epochs to the brain to create the stcs. Then at the stcs level, iterative over all epochs (the stc files in the stcs list) and generate the envelope (I think its just np.abs – but confirm this). Then average the stcs and you will have the average envelope of the beta frequency, which will be all positive.

I think that is what you are asking at least. You can also try it at the channel level first to confirm (and at the channel level, you can use the envelope=True to make it easier - just don’t project to the brain as an envelope for LCMV beamformer).

–Jeff