physio_filt_scanner_artifact

 filters scanner artifact for a given physio channel
 
   EEG = physio_filt_scanner_artifact(EEG,defs)
 
 This function uses EEGlab to epoch a given signal by TR, and then
 constructs a mean TR waveform and regresses the given signal by a design
 matrix containing a mean TR waveform regressor for each TR. It then
 removes the residuals from this regression from the signal of interest.
 
 REQUIRES
   EEG - an EEGlab data structure
   defs.physio_filt_scanner_artifact.channel - name of channel to be
       filtered, as it appears in EEG.chanlocs(:).labels
   defs.TR - length of TR, in seconds
 
 RETURNS
   EEG - an EEGlab data structure, with the given signal having been
   filtered
 
 FB 2009.04.20

0001 function EEG = physio_filt_scanner_artifact(EEG,defs)
0002 
0003 % filters scanner artifact for a given physio channel
0004 %
0005 %   EEG = physio_filt_scanner_artifact(EEG,defs)
0006 %
0007 % This function uses EEGlab to epoch a given signal by TR, and then
0008 % constructs a mean TR waveform and regresses the given signal by a design
0009 % matrix containing a mean TR waveform regressor for each TR. It then
0010 % removes the residuals from this regression from the signal of interest.
0011 %
0012 % REQUIRES
0013 %   EEG - an EEGlab data structure
0014 %   defs.physio_filt_scanner_artifact.channel - name of channel to be
0015 %       filtered, as it appears in EEG.chanlocs(:).labels
0016 %   defs.TR - length of TR, in seconds
0017 %
0018 % RETURNS
0019 %   EEG - an EEGlab data structure, with the given signal having been
0020 %   filtered
0021 %
0022 % FB 2009.04.20
0023 
0024 try TR = defs.TR; catch error('no TR defined in defs\n'); end
0025 try channel = defs.physio_filt_scanner_artifact.channel; catch
0026     error('no target channel defined in defs.physio...artifact.channel\n');
0027 end
0028 if iscell(channel), channel = channel{1}; end
0029 
0030 eeglab('initpaths');
0031 
0032 cidx = strmatch(channel,{EEG.chanlocs(:).labels});
0033 if isempty(cidx)
0034     error('target channel %s can not be found within EEG.chanlocs\n',...
0035         channel);
0036 elseif length(cidx) > 1
0037     warning('many channels match target channel %s, using first one\n',...
0038         channel);
0039     cidx = cidx(1);
0040 end
0041 
0042 % create TR event channel
0043 fprintf(1,'creating TR event channel\n');
0044 tr_idxs = 1:TR*EEG.srate:EEG.pnts; % in samples
0045 ntr = length(tr_idxs);
0046 tr_events = [ones(ntr,1) (tr_idxs/EEG.srate)']; % in seconds
0047 EEGe = pop_importevent(EEG,'append','no','event',tr_events,...
0048   'fields',{'type','latency'},'timeunit',1);
0049 
0050 % extract TR epochs
0051 fprintf(1,'extracting epochs\n');
0052 EEGe = pop_epoch(EEGe,{},[0 TR],'epochinfo','yes');
0053 
0054 % calculate mean TR waveform for each channel
0055 fprintf(1,'calculating mean TR waveform for %s channel, ',channel);
0056 W = [];
0057 for it=1:EEGe.trials
0058   W = [W; EEGe.data(cidx,:,it)];
0059 end
0060 mW = mean(W); % mean waveform across TRs
0061 
0062 % create regression design matrix
0063 dM = zeros(EEG.pnts,EEGe.trials);
0064 fprintf(1,'dM, ');
0065 for it=1:EEGe.trials
0066   start = (EEGe.pnts*(it-1)+1);
0067   stop  = EEGe.pnts*it;
0068   dM(start:stop,it) = mW';
0069 end
0070 
0071 % are there any samples not accounted for at the end of dM?
0072 if stop < size(dM,1)
0073   if (size(dM,1) - (stop - 1)) > length(mW)
0074     % left-over size is greater than the length of a TR
0075     error('there are less trials than there should be\n');
0076   else
0077     % fill out extra time at the end of the design matrix
0078     start = stop + 1;
0079     stop = size(dM,1);
0080     dM(start:stop,end+1) = mW(1:stop-start+1);
0081   end
0082 end
0083 dM(:,end+1) = 1; % must place a constant as the last regressor
0084 
0085 fprintf(1,'regressing, ');
0086 [b,bint,r,rint,stats] = regress(EEG.data(cidx,:)',dM);
0087 EEG.data(cidx,:) = EEG.data(cidx,:) - r';