generate_physio_evs

 phys = generate_physio_evs(phys);

 Generates explanatory variables (EVs/regressors) with physiological data.
 These regressors can then be loaded into regression models for fMRI data.

 phys - a structure containing all of the relevant data and metadata

 See also: populate_phys_struct, proc_physio

0001 function phys = generate_physio_evs(phys)
0002 % phys = generate_physio_evs(phys);
0003 %
0004 % Generates explanatory variables (EVs/regressors) with physiological data.
0005 % These regressors can then be loaded into regression models for fMRI data.
0006 %
0007 % phys - a structure containing all of the relevant data and metadata
0008 %
0009 % See also: populate_phys_struct, proc_physio
0010 
0011 % 11/01/05 Petr Janata
0012 
0013 % copy some often used values from the structure
0014 logfid = phys.ps.logfid;
0015 
0016 subid = phys.ps.sinfo.id;
0017 ri = phys.ri;
0018 pp = phys.ps.pp;
0019 
0020 nruns = length(ri);
0021 nslice_per_vol = phys.ps.pp.nslice_per_vol;
0022 
0023 % Determine how many sets of regressors we are generating
0024 nsets = length(phys.out_types);
0025 
0026 % Loop over all of the runs
0027 for irun = 1:nruns
0028   % We need to refer to runs in terms of their original indices, i.e. order in
0029   % which they were acquired.
0030   runid = phys.ps.sinfo.use_runs(irun); 
0031   
0032   nslices = length(phys.ri(irun).slice_onsets);
0033   nvol = phys.ps.sinfo.nvol(runid);
0034 
0035   % Perform a sanity check
0036   if nslices ~= nvol*nslice_per_vol
0037     fprintf(logfid,'ERROR:generate_physio_evs:Number of slices (%d) does not match expected number of slices (%d)', nslices, nvol*nslice_per_vol);
0038     error
0039   end
0040                
0041   ev_info = [];
0042   
0043   % Loop over all of the desired regressor sets
0044   for iset = 1:nsets
0045     set_type = phys.out_types{iset};
0046     outstub = sprintf('%s_run%d_%s', subid, runid, set_type);
0047 
0048     % Keep track of EV info in the ri (run info) structure
0049     ev_info(iset).set_type = set_type;
0050 
0051     % Do any necessary preprocessing and setting up of variables
0052     switch set_type
0053       case 'cardiac'
0054     curr_regress = zeros(nslices,6);
0055     ev_info(iset).ev_names = {'Sin1','Sin2','Sin3','Cos1','Cos2','Cos3'};
0056       case {'respir', 'respiration'}
0057     curr_regress = zeros(nslices,1);
0058     ev_info(iset).ev_names = {'Respiration'};
0059       otherwise
0060     end
0061     
0062     % Loop over all of the slices in the volume
0063     for islice = 1:nslices
0064       % Determine the sample offset of this slice on the original physiological
0065       % data timescale
0066       samp_idx = fix(ri(irun).slice_onsets(islice)*pp.Fs)+1;
0067 
0068       switch set_type
0069     case 'cardiac'
0070       card1_idx = max(find(ri(irun).cardiac <= ri(irun).slice_onsets(islice)));
0071       if isempty(card1_idx)
0072         card1 = 0;
0073       else
0074         card1 = ri(irun).cardiac(card1_idx);
0075       end
0076       card2_idx = min(find(ri(irun).cardiac > ...
0077           ri(irun).slice_onsets(islice)));
0078 
0079       % If physio collecting is turned off part way through last cardiac
0080       % cycle, we need to estimate when in the last cardiac cycle the slice occurred.
0081       if isempty(card2_idx)
0082         next_card = ri(irun).cardiac(end) + ...
0083         mean(diff(ri(irun).cardiac(end-9:end)));
0084         % Make sure nothing is really wrong
0085         if ri(irun).slice_onsets(end) > next_card
0086           error('Something is wrong with cardiac data')
0087         else
0088           card2 = next_card;
0089         end
0090       else
0091         card2 = ri(irun).cardiac(card2_idx);
0092       end
0093       % As described in Dagli et al. 1999, NeuroImage 9, 407-415
0094       theta = (ri(irun).slice_onsets(islice) - card1)/(card2-card1);
0095       curr_regress(islice,:) = [sin(2*pi*theta*(1:3)) cos(2*pi*theta*(1:3))];
0096       
0097     case {'respir', 'respiration'}
0098       % Figure out the range of samples in the respiration waveform (at the
0099       % sampling rate of the physio file) that we need to consider with
0100       % regard to the acquisition of this particular slice
0101       start_samp = samp_idx-phys.ps.nsamp_resp_peri;
0102       if start_samp < 1
0103         start_samp = 1;
0104       end
0105       stop_samp = samp_idx+phys.ps.nsamp_resp_peri;
0106       if stop_samp > length(ri(irun).respir)
0107         stop_samp = length(ri(irun).respir);
0108       end
0109       
0110       % Calculate the value for this slice
0111       curr_regress(islice) = mean(ri(irun).respir(start_samp:stop_samp));
0112       
0113     otherwise
0114       fprintf(logfid, 'Do not know how to handle regressor set: %s\n', set_type);
0115       end % switch phys.out_types{iset}
0116     end % for islice=
0117 
0118     % Post-process the regressors
0119     switch set_type
0120       case {'respir', 'respiration'}
0121     curr_regress = detrend(curr_regress);
0122       otherwise
0123     end
0124   
0125     % Rearrange into a 3D design matrix which is of the form
0126     % nvol X nev X nslice
0127     regset = reshape(curr_regress, nslice_per_vol, nvol, size(curr_regress,2));
0128     
0129     regset = permute(regset,[2 3 1]);
0130     
0131     % Correct for slice acquisition order
0132     regset(:,:,phys.ps.slice_order) = regset;
0133     nev = size(regset,2);
0134     ev_info(iset).nev = nev;
0135   
0136     for islice = 1:nslice_per_vol
0137       for iev = 1:nev
0138     outfname = fullfile(phys.outdir, sprintf('%s_ev%d_slice%02d.txt', outstub, iev, islice));
0139     fprintf(logfid,'Saving EV to file: %s\n', outfname);
0140     
0141     ev_info(iset).fnames{islice,iev} = outfname;
0142     
0143     outdata = regset(:,iev,islice);
0144     switch phys.file_format
0145       case 'fsl'
0146         save(outfname,'outdata','-ascii')
0147       otherwise
0148         fprintf(logfid,'Do not know how to handle desired EV output format: %s\n', phys.file_format);
0149     end % switch phys.file_format
0150       end % for iev
0151     end % for islice
0152   end % for iset=
0153   
0154   phys.ri(irun).ev_info = ev_info;
0155 end % for irun=