Description of ensemble

0001 function outdata = ensemble_physio(indata,defs)
0002 
0003 % processes physiological data, stores .mat file on disk, 1 per run
0004 %
0005 % If you are linking physio data to presentation file data, make sure to
0006 % process the presentation file data first, so that it is available to this
0007 % script.
0008 %
0009 % This script does not preprocess or normalize, but it does cut up and
0010 % limit per run.
0011 %
0012 % NOTE: this script also exports data to EEGLAB .set format
0013 % NOTE: this script assumes a slightly different format for the 'ri' physio
0014 % data format that has been used in the past, in read_keithley and
0015 % read_mate ... whereas in those scripts, channels were saved (such as
0016 % 'cardiac' and 'respir') as fields of the main 'ri' struct, but in the new
0017 % 'modified' format, the data fields are saved in ri.signal
0018 % (ri.signal.respir, ri.signal.cardiac), so that meta data and header
0019 % information can be saved as fields of 'ri'
0020 %
0021 % REQUIRES
0022 %   defs.sinfo(:).sessinfo(:).physio
0023 %   defs.sinfo(:).sessinfo(:).physio.source
0024 %   defs.sinfo(:).sessinfo(:).physio.resp_titles
0025 %   defs.physio.SIG_CHECK
0026 %   defs.physio.LINK2FMRI
0027 %
0028 % 2008/10/06 FB - started coding
0029 % 2009/06/18 FB - moved from ensemble_fmri_physio to ensemble_physio, added
0030 %   LINK2FMRI flag, to allow for generalization outside of fmri analyses
0031 
0032 outdata = ensemble_init_data_struct();
0033 
0034 global r
0035 
0036 r = init_results_struct;
0037 
0038 r.type = 'physio';  % Identify the type of this reporting instance
0039 r.report_on_fly = 1;
0040 
0041 %%%% FIXME: allow for a param?
0042 % load colormap
0043 load new_seismic;
0044 colormap(new_seismic);
0045 
0046 % Parse out the input data
0047 for idata = 1:length(indata)
0048   if isfield(indata{idata},'type')
0049     switch indata{idata}.type
0050       case {'paths'}
0051         pathdata = indata{idata};
0052         pcol = set_var_col_const(pathdata.vars);
0053       case {'sinfo'}
0054         sinfo = indata{idata};
0055         sinfo = sinfo.data;
0056         proc_subs = {sinfo(:).id};
0057         nsub_proc = length(proc_subs);
0058     end
0059   end
0060 end
0061 
0062 % check for required vars
0063 check_vars = {'sinfo','pathdata'};
0064 check_required_vars;
0065 
0066 if (iscell(indata) && ~isempty(indata) && isfield(indata{1},'task') && ...
0067         ~isempty(strmatch('return_outdir',indata{1}.task))) || ...
0068         (isstruct(indata) && isfield(indata,'task') && ...
0069         ~isempty(strmatch('return_outdir',indata.task)))
0070   if exist('pathdata','var') && ~isempty(pathdata.data{1})
0071     if length(nsub_proc) == 1
0072       pfilt = struct();
0073       pfilt.include.all.subject_id = proc_subs;
0074       lpathdata = ensemble_filter(pathdata,pfilt);
0075       if ~isempty(lpathdata.data{1})
0076         sfilt = pfilt;
0077         sfilt.include.all.path_type = {'physio_outdir'};
0078         spathdata = ensemble_filter(lpathdata,sfilt);
0079         if length(spathdata.data{1}) == 1
0080           % one epi outdir, save outdata = epi_outdir
0081           outdata = spathdata.data{pcol.path}{1};
0082         else
0083           sfilt = pfilt;
0084           sfilt.include.all.path_type = {'sess_outdir'};
0085           spathdata = ensemble_filter(lpathdata,sfilt);
0086           if length(spathdata.data{1}) == 1;
0087             outdata = spathdata.data{pcol.path}{1};
0088           else
0089             sfilt = pfilt;
0090             sfilt.include.all.path_type = {'sub_outdir'};
0091             spathdata = ensemble_filter(lpathdata,sfilt);
0092             if length(spathdata.data{1}) == 1;
0093               outdata = spathdata.data{pcol.path}{1};            
0094             end
0095           end
0096         end
0097       end
0098     end
0099   end
0100   if ~exist('outdata','var') || ~exist(outdata,'dir'), outdata = ''; end
0101   return
0102 end
0103 
0104 outdata.vars = [outdata.vars 'sinfo'];
0105 sinfo_idx = length(outdata.vars);
0106 outdata.data{sinfo_idx} = ensemble_init_data_struct();
0107 outdata.data{sinfo_idx}.type = 'sinfo';
0108 outdata.data{sinfo_idx}.data = sinfo;
0109 
0110 outdata.vars = [outdata.vars 'physio_paths'];
0111 ppaths_idx = length(outdata.vars);
0112 outdata.data{ppaths_idx} = ensemble_init_data_struct();
0113 outdata.data{ppaths_idx}.type='physio_paths';
0114 outdata.data{ppaths_idx}.vars = {'subject_id','session',...
0115     'ensemble_id','run','path','preprocessed'};
0116 outdata.data{ppaths_idx}.data{1} = {};
0117 outdata.data{ppaths_idx}.data{2} = {};
0118 outdata.data{ppaths_idx}.data{3} = [];
0119 outdata.data{ppaths_idx}.data{4} = [];
0120 outdata.data{ppaths_idx}.data{5} = {};
0121 outdata.data{ppaths_idx}.data{6} = [];
0122 ppcol = set_var_col_const(outdata.data{ppaths_idx}.vars);
0123 
0124 ecdparams.outDataName = 'physio_data';
0125 
0126 try SIG_CHECK = defs.SIG_CHECK; catch SIG_CHECK = 1; end
0127 try LINK2FMRI = defs.LINK2FMRI; catch LINK2FMRI = 0; end
0128 
0129 % init EEGLAB
0130 
0131 % eeglab('initpaths');
0132 % Add only the eeglab paths we need
0133 eeglab_root = fileparts(which('eeglab'));
0134 addpath(fullfile(eeglab_root,'functions/popfunc'));
0135 addpath(fullfile(eeglab_root,'functions/sigprocfunc'));
0136 
0137 %
0138 % START OF THE SUBJECT LOOP
0139 %
0140 
0141 for isub=1:nsub_proc
0142   subid = sinfo(isub).id;
0143   msg = sprintf('\t\tPROCESSING SUBJECT (%d/%d): %s\n', isub, nsub_proc,subid);
0144   r = update_report(r,msg);
0145 
0146   % Determine number of sessions for this subject
0147   nsess = length(sinfo(isub).sessinfo);
0148   
0149   %
0150   % START OF THE SESSION LOOP
0151   %
0152   
0153   for isess = 1:nsess
0154     sess = sinfo(isub).sessinfo(isess);
0155     
0156     if isfield(sess,'use_session') && ~sess.use_session
0157       msg = sprintf('\t\t\tSkipping session %s\n',sess.id);
0158       r = update_report(r,msg);
0159       continue
0160     elseif ~isfield(sess,'physio') || ~isstruct(sess.physio)
0161       msg = sprintf(['\t\t\tNo Physio sinfo for session %s, '...
0162           'SKIPPING!\n'],sess.id);
0163       r = update_report(r,msg);
0164       continue
0165     end
0166     
0167     r = update_report(r,sprintf('\t\t\t%s\n', sess.id));
0168     
0169     % init session vars
0170     physio = sinfo(isub).sessinfo(isess).physio;
0171     pparams = init_physmon_params(physio.source);
0172     pparams.channels = physio.channels;
0173     if isfield(physio,'baseline')
0174       pparams.baseline = physio.baseline;
0175     end
0176     if isfield(physio,'link2pres')
0177       pparams.link2pres = physio.link2pres;
0178     end
0179     
0180     if LINK2FMRI
0181       % Obtain information about the scanning protocol that was used in this
0182       % session
0183       protocol_idx = strmatch(sess.protocol_id, {defs.fmri.protocol.id}, ...
0184           'exact');
0185       pparams.TR = defs.fmri.protocol(protocol_idx).epi.tr;
0186       pparams.nslice_per_vol = defs.fmri.protocol(protocol_idx).epi.nslices;
0187     end
0188     
0189     
0190     % get physio_indir/outdir, behav_outdir
0191     sfilt = struct();
0192     sfilt.include.all.subject_id = {subid};
0193     sfilt.include.all.session = {sess.id};
0194     spaths = ensemble_filter(pathdata,sfilt);
0195     
0196     indx = strmatch('physio_indir',spaths.data{pcol.path_type});
0197     physio_indir = spaths.data{pcol.path}{indx};
0198 
0199     outdx = strmatch('physio_outdir',spaths.data{pcol.path_type});
0200     physio_outdir = spaths.data{pcol.path}{outdx};
0201     
0202     behdx = strmatch('behav_outdir',spaths.data{pcol.path_type});
0203     if isempty(behdx)
0204       behav_outdir = physio_outdir;
0205     else
0206       behav_outdir = spaths.data{pcol.path}{behdx};
0207     end
0208     
0209     % Determine how many runs we're dealing with
0210     if LINK2FMRI
0211       runs = sess.use_epi_runs;
0212     else
0213       if isfield(sess,'use_physio_runs')
0214         runs = sess.use_physio_runs;
0215       else
0216         runs = 1:size(physio.datafiles,1);
0217       end
0218     end
0219     nruns = length(runs);  
0220     
0221     %
0222     % START OF RUN LOOP
0223     %
0224     for irun = 1:nruns
0225 
0226       msg = sprintf('\t\tPROCESSING RUN (%d/%d)\n', irun,nruns);
0227       r = update_report(r,msg);
0228         
0229       % get physio data
0230       physiofname = physio.datafiles{runs(irun),1};
0231       if isempty(physiofname)
0232         msg = sprintf('no physiofname found, skipping run %d\n',irun);
0233         r = update_report(r,msg);
0234         continue
0235       end
0236       physiopath = fullfile(physio_indir,physiofname);
0237       if ~exist(physiopath,'file')
0238         msg = sprintf('physiopath doesn''t exist, skipping run %d\n',irun);
0239         r = update_report(r,msg);
0240         continue
0241       end
0242       lparams = pparams;
0243       rparams = physio.datafiles{runs(irun),2};
0244       if isfield(rparams,'link2pres')
0245         lparams.run=irun;
0246         lparams.link2pres = rparams.link2pres;
0247         lparams.link2pres.presfname = fullfile(behav_outdir,...
0248             sprintf('%s_%s_present.mat',subid,sess.id));
0249         lparams.link2pres.filt.include.any.RUN = irun;
0250       end
0251       if isfield(rparams,'signal_bounds')
0252         lparams.signal_bounds = rparams.signal_bounds;
0253       end
0254       
0255       % read-in physio data
0256       ri = proc_physio_data(physiopath,lparams);
0257 
0258       if isempty(fields(ri))
0259         continue
0260       end
0261       
0262       % import into EEGLAB
0263       snames = fields(ri.signal);
0264       nchans = length(snames);
0265       npts   = length(ri.signal.(snames{1}));
0266       signal = zeros(nchans,npts);
0267       chanlocs = struct();
0268       for ic=1:nchans
0269         signal(ic,:) = ri.signal.(snames{ic});
0270         chanlocs(ic).labels = snames{ic};
0271       end
0272       if isfield(ri.meta,'baseline_samps')
0273         xmin = -ri.meta.baseline_samps/ri.meta.srate;
0274       else
0275         xmin = 0;
0276       end
0277       
0278       EEG = pop_importdata('dataformat','array','data',signal,'srate',...
0279           ri.meta.srate,'subject',subid,'session',sess.ensemble_id,...
0280           'chanlocs',chanlocs,'nbchan',nchans,'pnts',npts,'xmin',xmin);
0281       
0282       set_fname = sprintf('%s_%s_run%d_physio.set',subid,sess.id,irun);
0283       pop_saveset(EEG,'filename',set_fname,'filepath',physio_outdir);
0284       set_fpn = fullfile(physio_outdir,set_fname);
0285 
0286       outdata.data{ppaths_idx} = ensemble_add_data_struct_row(...
0287           outdata.data{ppaths_idx},'subject_id',subid,'session',sess.id,...
0288           'ensemble_id',sess.ensemble_id,'run',irun,'path',set_fpn,...
0289           'preprocessed',0);
0290       
0291       if SIG_CHECK
0292         sigs = {EEG.chanlocs(:).labels};
0293         nsig = length(sigs);
0294         figure();
0295         for i=1:nsig
0296           sig = sigs{i};
0297           subplot(nsig,1,i);
0298           plot(EEG.data(i,abs(EEG.xmin*EEG.srate)+1:end));
0299           ylabel(sig);
0300         end % for i=1:nsig
0301         if isfield(defs,'figs') && isfield(defs.figs,'write2file') && ...
0302                 defs.figs.write2file
0303             figfname = fullfile(physio_outdir,sprintf(...
0304                 '%s_%s_run%d_physioSigCheck.ps',subid,sess.id,irun));
0305             print(figfname,defs.figs.printargs{:});
0306         end
0307       end % if SIG_CHECK
0308       
0309       %%%% Filter this run data?
0310       if ~isempty(strmatch('filter',fieldnames(lparams.channels)))
0311         EEGf = EEG;
0312         for ic=1:length(lparams.channels)
0313           if isfield(lparams.channels(ic),'filter')
0314             % get EEG channel index for this channel
0315             cidx = strmatch(lparams.channels(ic).name,...
0316                 {EEGf.chanlocs(:).labels});
0317             if length(cidx) > 1
0318               cidx = cidx(1);
0319             elseif isempty(cidx)
0320               warning(['no match between lparams.channel %d and EEG ',...
0321                   'struct!? skipping\n'],ic);
0322               continue
0323             end % if length(cidx
0324           
0325             % filter data?
0326             lf = lparams.channels(ic).filter;
0327             if ~isstruct(lf)
0328               warning(['filter struct for channel %d is not a struct, '...
0329                   'skipping\n'],ic);
0330             else
0331               for icf=1:length(lf)
0332                 if ~isfield(lf(icf),'fun')
0333                   warning(['filter struct %d for channel %d does not '...
0334                       'contain a function handle (.fun field), '...
0335                       'skipping\n'],icf,ic);
0336                 else
0337                   lfh = parse_fh(lf(icf).fun);
0338                   fname = func2str(lfh);
0339 
0340                   cparams = lparams;
0341                   cparams.(fname) = lf(icf);
0342                   EEGf = lfh(EEGf,cparams);
0343                 end % if ~isfield(lf(icf
0344               end % for icf=1:
0345             end % if ~isstruct(lf
0346           end % if isfield(lparams.channels(ic
0347         end % for ic=1:length(lparams.ch
0348 
0349         set_fname = sprintf('%s_%s_run%d_physio_filtered.set',subid,sess.id,irun);
0350         pop_saveset(EEGf,'filename',set_fname,'filepath',physio_outdir);
0351         set_fpn = fullfile(physio_outdir,set_fname);
0352 
0353         outdata.data{ppaths_idx} = ensemble_add_data_struct_row(...
0354             outdata.data{ppaths_idx},'subject_id',subid,'session',sess.id,...
0355             'ensemble_id',sess.ensemble_id,'run',irun,'path',set_fpn,...
0356             'preprocessed',1);
0357         
0358         if SIG_CHECK
0359           sigs = {EEGf.chanlocs(:).labels};
0360           nsig = length(sigs);
0361           figure();
0362           for i=1:nsig
0363             sig = sigs{i};
0364             subplot(nsig,1,i);
0365             plot(EEGf.data(i,abs(EEGf.xmin*EEGf.srate)+1:end));
0366             ylabel(sig);
0367           end % for i=1:nsig
0368           if isfield(defs,'figs') && isfield(defs.figs,'write2file') && ...
0369                   defs.figs.write2file
0370               figfname = fullfile(physio_outdir,sprintf(...
0371                   '%s_%s_run%d_physioFilteredSigCheck.ps',subid,sess.id,irun));
0372               print(figfname,defs.figs.printargs{:});
0373           end
0374         end % if SIG_CHECK
0375       end % if ~isempty(strmatch('filter
0376       
0377     end % for irun
0378     
0379   end % for isess
0380 end % for isub=
ensemble_physio

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
