MNE处理脑电数据

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mne
import os.path as op

import mne
from mne.time_frequency import tfr_morlet, psd_multitaper, psd_welch
from mne.datasets import somato

from mne.channels import compute_native_head_t, read_custom_montage
from mne.datasets import sample
from mne import read_source_estimate


'''读取数据，这是刺激下采集的数据，所以也会有事件'''

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_raw.fif')

raw = mne.io.read_raw_fif(sample_data_raw_file, preload=False)


raw.pick(['eeg']).load_data()  # 挑选EEG通道
raw.load_data()  # 挑选EEG通道

raw.plot(duration=4, show_scrollbars=False)


''' 预处理 '''

''' 滤波 '''

raw.filter(l_freq=0.1, h_freq=70)

''' 重采样 '''

raw.resample(sfreq=256)


''' ICA '''
ica_raw = raw.copy()
ica = mne.preprocessing.ICA(n_components=20, random_state=97, max_iter=800)
ica.fit(ica_raw)
ica.exclude = [0, 18, 19]

# ica.plot_sources(ica_raw, show_scrollbars=False)

ica.apply(ica_raw)

ica_raw.plot(duration=4, show_scrollbars=False)

''' power spectral density '''

raw.plot_psd(fmax=70)

raw.plot_psd(fmin=0, fmax=4)
raw.plot_psd(fmin=4, fmax=8)
raw.plot_psd(fmin=8, fmax=12)
raw.plot_psd(fmin=12, fmax=30)
raw.plot_psd(fmin=30, fmax=45)

# raw.plot_psd_topo(proj=True)




''' 等分数据，转为epoch '''

new_events = mne.make_fixed_length_events(raw, start=20, stop=70, duration=10.)

epochs = mne.Epochs(raw, new_events)
epochs = epochs.set_eeg_reference(projection=True)
# print(epochs)
# input()

''' 频率分析 '''

# epochs.plot_psd_topomap(normalize=False)


''' 时频分析 '''

''' 要研究frequencies的作用'''

frequencies = np.arange(7, 30, 3)
power = mne.time_frequency.tfr_morlet(epochs, n_cycles=2, return_itc=False,
                                      freqs=frequencies, decim=3)
power.plot('EEG 035')


''' Inspect power '''

# power.plot_topo(baseline=(-0.5, 0), mode='logratio', title='Average power')
# power.plot([35], baseline=(-0.5, 0), mode='logratio', title=power.ch_names[35])

fig, axis = plt.subplots(1, 2, figsize=(7, 4))
power.plot_topomap(ch_type='eeg', fmin=8, fmax=12,
                   baseline=(-0.5, 0), mode='logratio', axes=axis[0],
                   title='Alpha', show=False)
power.plot_topomap(ch_type='eeg', fmin=12, fmax=30,
                   baseline=(-0.5, 0), mode='logratio', axes=axis[1],
                   title='Beta', show=False)
power.plot_topomap(ch_type='eeg', fmin=4, fmax=8,
                   baseline=(-0.5, 0), mode='logratio', axes=axis[2],
                   title='Theta', show=False)
power.plot_topomap(ch_type='eeg', fmin=1, fmax=4,
                   baseline=(-0.5, 0), mode='logratio', axes=axis[3],
                   title='Delta', show=False)

mne.viz.tight_layout()
plt.show()


''' inverse modeling '''

noise_cov = mne.compute_covariance(
    epochs, tmax=0., method=['shrunk', 'empirical'], rank=None, verbose=True)

fig_cov, fig_spectra = mne.viz.plot_cov(noise_cov, raw.info)

evoked = epochs.average().pick('meg')
evoked.plot(time_unit='s')
evoked.plot_topomap(times=np.linspace(0.05, 0.15, 5), ch_type='mag',
                    time_unit='s')