Source localization from resting state data

[dokato]

Hi, in my project I'm trying to perform source recontruction from resting state MEG data. I tried to follow the steps from this tutorial: https://martinos.org/mne/stable/auto_tutorials/plot_mne_dspm_source_localization.html

with one modification. As this is resting state conitnuous recording I don't have events in my data, so I create 2 seconds artificial epochs this way:

def create_fake_epochs(raw, n_epochs):
    Fs = raw.info['sfreq']
    rdata = raw.get_data()
    events = np.array([[int(i*rdata.shape[1]/n_epochs), 0, 1] for i in range(n_epochs)])
    epochs = mne.Epochs(raw, events, 1, 0, int(0.9*(rdata.shape[1]/n_epochs)/Fs),
                        proj=True)
    return epochs

Everything goes smooth, but when I computed labelled time series after parcellation:

    labels = mne.read_labels_from_annot(subject_number, parc='aparc', subjects_dir=head_shape_dir)
    label_names = [l.name for l in labels]
    label_ts = mne.extract_label_time_course(stc, labels, inverse_operator['src'], mode='mean_flip',
                                             return_generator=False)

And plotted power spectrum in occipital region of the brain I don't get any alpha peak. Power spectrum looks rather flat, like white noise. Note the scale difference on the plot above.

I'm struggling to find a mistake in my code. I tried both sLORETA and Beamforming, but all lead to the same results. Do you have any idea, why? Am I doing something conceptually wrong?

[agramfort]

you should be using:

https://martinos.org/mne/stable/generated/mne.make_fixed_length_events.html

how do you estimate power? did you remove bad data segments?

[dokato]

Thanks for that function - it will make life easier.

For sake of this visualisation I used custom function:

def periodogram(s, window, F_samp):
    s = s*window
    N_fft = len(s)
    S = np.fft.rfft(s,N_fft)
    P = S*S.conj()/np.sum(window**2)   
    P = P.real/F_samp
    F = np.fft.rfftfreq(N_fft, 1/F_samp)
    if len(s)%2 ==0:
        P[1:-1] *=2
    else:
        P[1:] *=2
    plt.plot(F,P)
    plt.xlim([0,int(F_samp/2)])
    plt.show()

All steps of data preprocessing I'm using is to notch-filter (50Hz) the data, highpass filter from 1.5Hz, resample to 256 Hz, visual inspection and bad channels removal, ICA artefacts rejection (ECG and EOG), but then after creating my epochs I didn't remove any segments. I might try if that helps... but I had almost 200 2sec long epochs, so the change shouldn't be that obvious.

[agramfort]

use this function:

https://martinos.org/mne/dev/generated/mne.time_frequency.psd_array_welch.html

never ever use FFT to estimate a good PSD. It has too much variance.

if you have short signals you could use:

https://martinos.org/mne/dev/generated/mne.time_frequency.psd_array_multitaper.html

take home message: avoid rewriting too much stuff from scratch.

HTH

[dokato]

never ever use FFT to estimate a good PSD. It has too much variance.

Yeah, I know, it was just a quick sanity check, not that fluent in mne yet.

take home message: avoid rewriting too much stuff from scratch.

Thanks (but I guess no one likes black boxes).

So now I used psd_array_welch and applied rejection reject = dict(mag=4e-12) while creating epochs. Still getting:

RHS is reconstructed lateraloccipital region.

[dokato]

Not sure how relevant that is to this problem, but mag noise std looks fine:

[dokato]

Okay, my bad... I was using apply_inverse instead of apply_inverse_raw. Now everything works fine. Thanks for all the hints.