Home > database > physio > ensemble_physio.m

ensemble_physio

PURPOSE ^

processes physiological data, stores .mat file on disk, 1 per run

SYNOPSIS ^

function outdata = ensemble_physio(indata,defs)

DESCRIPTION ^

 processes physiological data, stores .mat file on disk, 1 per run
 
 If you are linking physio data to presentation file data, make sure to
 process the presentation file data first, so that it is available to this
 script.
 
 This script does not preprocess or normalize, but it does cut up and
 limit per run.
 
 NOTE: this script also exports data to EEGLAB .set format
 NOTE: this script assumes a slightly different format for the 'ri' physio
 data format that has been used in the past, in read_keithley and
 read_mate ... whereas in those scripts, channels were saved (such as
 'cardiac' and 'respir') as fields of the main 'ri' struct, but in the new
 'modified' format, the data fields are saved in ri.signal
 (ri.signal.respir, ri.signal.cardiac), so that meta data and header
 information can be saved as fields of 'ri'
 
 REQUIRES
   defs.sinfo(:).sessinfo(:).physio
   defs.sinfo(:).sessinfo(:).physio.source
   defs.sinfo(:).sessinfo(:).physio.resp_titles
   defs.physio.SIG_CHECK
   defs.physio.LINK2FMRI
 
 2008/10/06 FB - started coding
 2009/06/18 FB - moved from ensemble_fmri_physio to ensemble_physio, added
   LINK2FMRI flag, to allow for generalization outside of fmri analyses

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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=

Generated on Wed 22-May-2019 04:00:51 by m2html © 2003