ADD examples: CSSL pipeline

examples/CSSL/1 bad channels.py

@@ -0,0 +1,42 @@
+# skip test
+"""Load data from *.sqd file, look for bad channels and filter data"""
+import mne
+from eelbrain import *
+
+# load the raw data
+raw = mne.io.read_raw_kit('Data/R2290_HAYO_P3H_08.26.2015.sqd',
+                          stim_code='channel', stim=range(162, 176),
+                          preload=True)
+
+# filter the raw data
+raw.filter(1, 8)
+
+# load the events
+ds = load.fiff.events(raw)
+# look at events
+print(ds)
+# select the desired events
+ds = ds.sub("trigger == 167")
+
+# Add MEG data for trials
+ds['meg'] = load.fiff.epochs(ds, 0, 60)
+# concatenate MEG data into one long "trial"
+meg = concatenate(ds['meg'])
+# plot average correlation with neighbors
+plot.Topomap(neighbor_correlation(meg))
+
+# add bad channels
+raw.info['bads'] = ['MEG 056']
+
+# check the result of removing the channel (need to add the data with the
+# new bad channel setting)
+ds['meg'] = load.fiff.epochs(ds, 0, 60)
+plot.Topomap(neighbor_correlation(concatenate(ds['meg'])))
+
+# remove the mean from the data
+ds['meg'] -= ds['meg'].mean('time')
+# plot the first trial (ds[0]) as continuous data
+plot.TopoButterfly(ds[0, 'meg'])
+
+# save the filtered raw data
+raw.save('Data/R2290_HAYO_P3H_1-8-raw.fif')

examples/CSSL/2 dss.py

@@ -0,0 +1,20 @@
+# skip test
+"""Compute DSS on Alex's Toolbox demo data"""
+from eelbrain import *
+
+
+# If the raw file is already filtered, the events can be loaded directly with
+# the filename
+ds = load.fiff.events('DATA/R2290_HAYO_P3H_1-8-raw.fif')
+# select the desired events
+ds = ds.sub("trigger == 167")
+# add MEG trial data to the dataset
+ds['meg'] = load.fiff.epochs(ds, 0, 60)
+
+# DSS from
+todss, fromdss = dss(ds['meg'])
+
+# plot the DSS topography
+plot.Topomap(fromdss[:, :6], '.dss', h=2, ncol=6)
+# save the DSS for later use
+save.pickle((todss, fromdss), 'DATA/R2290_HAYO_P3H_DSS.pickle')

examples/CSSL/3 trf.py

@@ -0,0 +1,60 @@
+# skip test
+"""Compute TRFs with boosting"""
+from eelbrain import *
+
+
+# load the data
+ds = load.fiff.events('Data/R2290_HAYO_P3H_1-8-raw.fif')
+ds = ds.sub("trigger == 167")
+ds['meg'] = load.fiff.epochs(ds, 0, 60)
+
+# load the saved DSS
+todss, fromdss = load.unpickle('Data/R2290_HAYO_P3H_DSS.pickle')
+
+# extract the source time course for the first DSS
+meg_dss1 = todss[0].dot(ds['meg'])
+# average trials (trial number is represented by the 'case' dimension)
+y = meg_dss1.mean('case')
+# MNE-Python includes last sample; discard it for consistency with wav
+y = y[:60]
+# resample at 100 Hz
+y = resample(y, 100)
+
+# load the stimulus
+wav = load.wav('Stimuli/P3H.wav')
+# compute envelope
+env = wav.envelope()
+# filter with same settings as MEG data
+x = filter_data(env, 1, 8)
+# resample to same sampling rate as MEG data
+x = resample(x, 100)
+# plot stimulus
+plot.UTS([[wav, env, x]])
+
+# boosting
+res = boosting(y, x, -0.1, 0.5)
+
+# plot TRF
+plot.UTS(res.h)
+
+# compute TRF for multiple predictors (background and mix)
+wav_bg = load.wav('Stimuli/English_F_Norm_One_S.wav')
+x_bg = resample(filter_data(wav_bg.envelope(), 1, 8), 100)
+x_bg.name = 'background'
+
+wav_mix = wav + wav_bg
+x_mix = resample(filter_data(wav_mix.envelope(), 1, 8), 100)
+x_mix.name = 'mix'
+
+res = boosting(y, [x, x_bg, x_mix], -0.1, 0.5, error='l1')
+
+# plot kernels separately
+plot.UTS(res.h, ncol=1)
+# plot kernels in one plot
+plot.UTS([res.h])
+
+# convolve the kernel with the predictors to predict the MEG response
+y_pred = convolve(res.h_scaled, [x, x_bg, x_mix])
+y_pred.name = "y_pred"
+# compare actual and predicted response
+plot.UTS([[y, y_pred]])