Contents

Add paths

%NOTE: BE SURE TO ADJUST PATHS TO YOUR MACHINE
clc
clear variables;
addpath('C:\MATLAB\mTRF_1.6');
addpath C:\MATLAB\fieldtrip-20170801\fieldtrip-20170801
ft_defaults

Initialize some variables

%initialize variables
fs = 250; %sample rate
nCh = 64; %number of eeg channels
sCh = 65; %location of stimulus channel
bpF = [4 8]; %band pass filter of EEG
bpF_o = 2; %band pass filter order of EEG

%stuff that we will eventually cross validate
trilen = 60; %"trial" length in seconds
lambda = 256; % ridge regression penalty coefficient, in powers of 2
intWin = [-100 350]; %integration window, in seconds


%matlab default colors for plotting.
C1 = [0.8500, 0.3250, 0.0980];
C2 = [0, 0.4470, 0.7410];

Import data and visualize

%load in a subject
load S04_raw;

%separate EEG and Audio
cfg = [];
cfg.channel = 1:nCh; %eeg channels
eeg = ft_selectdata(cfg,data);

cfg.channel = sCh;
stim = ft_selectdata(cfg,data);

%let's look at the EEG
cfg = [];
cfg.ylim = [-10 10];
cfg.blocksize = 20;
cfg.continuous = 'yes';
cfg.viewmode = 'vertical';
ft_databrowser(cfg,eeg);

figure;
%let's look at the stimulus
cfg = [];
cfg.xlim = [15 30]; %15 to 30 second chunk
cfg.graphcolor = C2;
ft_singleplotER(cfg,stim);
ylabel('Amplitude'); xlabel('Time')

%low pass filter EEG and reference
cfg = [];
cfg.bpfilter = 'yes';
cfg.bpfreq = bpF;
cfg.bpfiltord = bpF_o;
cfg.reref = 'yes';
cfg.refchannel = 'all';
eeg = ft_preprocessing(cfg,eeg);

%look at the low-pass filtered data
cfg = [];
cfg.ylim = [-10 10];
cfg.blocksize = 20;
cfg.continuous = 'yes';
cfg.viewmode = 'vertical';
ft_databrowser(cfg,eeg);

the call to "ft_selectdata" took 0 seconds
the call to "ft_selectdata" took 0 seconds
the input is raw data with 64 channels and 1 trials
detected   0 visual artifacts
the call to "ft_singleplotER" took 0 seconds
the call to "ft_selectdata" took 0 seconds
preprocessing
preprocessing trial 1 from 1

the call to "ft_preprocessing" took 3 seconds
the input is raw data with 64 channels and 1 trials
detected   0 visual artifacts

TRFs. Step 1: Encoding Models

%extract eeg and audio so they are simply matrices
data = eeg.trial{1};
stim = stim.trial{1};
data(sCh:end,:) = [];

%we're going to reshape the data so it is in arbitrary "trials" of 60
%seconds long.

nTri = floor((length(data)/fs) / trilen); % get N trials
data  = data(:,1:nTri * trilen * fs); %extract possible samps
data = reshape(data,[nCh, trilen * fs, nTri]); %reshape data
stim = stim(:,1:nTri * trilen * fs); %do the same for stimulus
stim = reshape(stim,[1,trilen * fs, nTri]);

%we are left with a matrix that is channels by time by "trials"

%put into mTRF-friendly data structure (cells) and z-score
dat = []; sig = []; sig2 = [];

for tr = 1:size(stim,3)
    dat{tr} = zscore(data(:,:,tr)');
    sig{tr} = zscore(squeeze(stim(:,:,tr))'); %forward speech


    %to understand if our TRFs represent real signals or spurious
    %noise, we will create a control condition in which the speech is
    %reversed in time. This preserves the statistics of the speech
    %signal but the brain response should be absent.

    sig2{tr} = flipud(zscore(squeeze(stim(:,:,tr))')); %reverse speech
end


% trial loop for "forward" speech (i.e., aligned speech).
% the output of this loop is a TRF for each of the trials. We also
% obtain regression constants used for predicting EEG from the model. A
% time vector is also produced.

%Pay attention to the 4th input of mTRF train below. 1 is an
% "encoding" model, while -1 is a backward or decoding model.

trfsF = []; cF = []; %predefine variable for "forward" TRFs
trfsR = []; cR = []; %predfine variable for reversed TRFs
for tr = 1:length(sig)
    [trfsF(tr,:,:),t,cF(tr,:)] = mTRFtrain(sig{tr},dat{tr},fs,1,intWin(1),intWin(2),lambda);
    [trfsR(tr,:,:),t,cR(tr,:)] = mTRFtrain(sig2{tr},dat{tr},fs,1,intWin(1),intWin(2),lambda);
end

%by averaging over trials (dim 1) we get an average model
modlF = mean(trfsF,1);  %forward model
modlR = mean(trfsR,1); %reverse model
cF = mean(cF,1); %regression constants for forward model
cR = mean(cR,1); %regression constants for reverse model

%let's have a look at the TRFs for each channel

figure;
subplot(2,1,1)
plot(t,squeeze(modlF),'color',C1) %plotting the forward aligned model
axis([-75 325 -.0075 .0075])
xlabel('Time Lag (ms)'); ylabel('TRF Weight (a.u.)');
title('Forward speech')

subplot(2,1,2)
plot(t,squeeze(modlR),'color',C2) %plotting the reversed aligned model
xlabel('Time Lag (ms)'); ylabel('TRF Weight (a.u.)');
title('Reversed speech')
axis([-75 325 -.0075 .0075])

set(gcf,'color','w');

% now we are going to use the mTRF predict function to predict EEG
% from the audio envelope, and the quality of this prediction is
% assessed by simple Pearson's R.

predF = []; predR = []; %predefine variables for pearson R
rF = []; rR = [];
for tr = 1:length(sig)
    [predF(:,:,tr),rF(tr,:)] = mTRFpredict(sig{tr},dat{tr},modlF,fs,1,intWin(1),intWin(2),cF);
    [predR(:,:,tr),rR(tr,:)] = mTRFpredict(sig2{tr},dat{tr},modlR,fs,1,intWin(1),intWin(2),cR);
end

%let's look at an example of how the predicted data looks with the real
%data

figure; plot(zscore(squeeze(predF(2e3:7e3,1,1))),'color',C1); hold on;
plot(zscore(squeeze(dat{1}(2e3:7e3,1))),'color',C2)
xlabel('Time (samples)'); ylabel('Amplitude (uV');
title(['Accuracy: r = ' num2str(rF(1,1))]);
xlim([0 5e3])

avAcc = [mean(mean(abs(rF))) mean(mean(abs(rR)))];
disp(['Average Encoding Accuracy for Forward Speech: r = ' num2str(avAcc(1))])
disp(['Average Encoding Accuracy for Reverse Speech: r = ' num2str(avAcc(2))])

Average Encoding Accuracy for Forward Speech: r = 0.045845
Average Encoding Accuracy for Reverse Speech: r = 0.034466

Step 2: Decoding Models

%The steps above are the same, except pay attention to the fact that the
%dimensions of the input strutures are a bit different. See mTRFpredict
%documentation.

%Note that the 4th input in mTRFtrain is now -1 for decoding.

%NOTE ALSO THAT DECODING MODELS TAKE LONGER TO COMPUTE
disp('Running decoding models. This process is slower than for encoding.')

%trial loop for "forward" speech (i.e., aligned speech)
trfsF = []; trfsR = [];
for tr = 1:length(sig)
    %Because this takes a while, it is sometimes good to let yourself know
    %where you are in the loop.
    disp(['Decoding Model. Trial: ' num2str(tr) ' of ' num2str(length(sig))])
    [trfsF(:,:,tr),t,cF] = mTRFtrain(sig{tr},dat{tr},fs,-1,intWin(1),intWin(2),lambda);
    [trfsR(:,:,tr),t,cR] = mTRFtrain(sig2{tr},dat{tr},fs,-1,intWin(1),intWin(2),lambda);
end


modlF = mean(trfsF,3);
modlR = mean(trfsR,3);
cF = mean(cF,2);
cR = mean(cR,2);

predF = []; predR = []; %predefine variables for pearson R
rF = []; rR = [];
for tr = 1:length(sig)
    [predF(:,tr),rF(tr,:)] = mTRFpredict(sig{tr},dat{tr},modlF,fs,-1,intWin(1),intWin(2),cF);
    [predF(:,tr),rR(tr,:)] = mTRFpredict(sig2{tr},dat{tr},modlR,fs,-1,intWin(1),intWin(2),cR);
end

figure; plot(zscore(squeeze(predF(5e3:10e3,2))),'color',C1); hold on;
plot(zscore(sig{2}(5e3:10e3)),'color',C2)
xlabel('Time (samples)'); ylabel('Amplitude (uV');
title(['Accuracy: r = ' num2str(rF(1))]);
xlim([0 5e3])

avAcc = [mean(abs(rF)) mean(abs(rR))];
disp(['Average Encoding Accuracy for Forward Speech: r = ' num2str(avAcc(1))])
disp(['Average Encoding Accuracy for Reverse Speech: r = ' num2str(avAcc(2))])

Running decoding models. This process is slower than for encoding.
Decoding Model. Trial: 1 of 16
Decoding Model. Trial: 2 of 16
Decoding Model. Trial: 3 of 16
Decoding Model. Trial: 4 of 16
Decoding Model. Trial: 5 of 16
Decoding Model. Trial: 6 of 16
Decoding Model. Trial: 7 of 16
Decoding Model. Trial: 8 of 16
Decoding Model. Trial: 9 of 16
Decoding Model. Trial: 10 of 16
Decoding Model. Trial: 11 of 16
Decoding Model. Trial: 12 of 16
Decoding Model. Trial: 13 of 16
Decoding Model. Trial: 14 of 16
Decoding Model. Trial: 15 of 16
Decoding Model. Trial: 16 of 16
Average Encoding Accuracy for Forward Speech: r = 0.11755
Average Encoding Accuracy for Reverse Speech: r = 0.060606

Published with MATLAB® R2018a