1 function [Y, T] = phys_noise(varargin)
   2 %PHYS_NOISE Summary of this function goes here
   3 %   Detailed explanation goes here
   4 
   5 ppg_file   = varargin{1};
   6 dic_folder = varargin{2};
   7 TR         = varargin{4};
   8 Hz         = varargin{10};
   9 
  10 %% start
  11 P.TR        = TR; % Acquisition time (for working around dummy scans)
  12 P.EndClip   = 1000; % Always seems to take a second and change to finish
  13 P.StartClip = 0;    % Some systems have to have a pause between record onset and gradient onset
  14 P.Hz        = Hz;   % Sampling Frequency
  15 P.StartJump = 0;    % Number of shift samples to adjust at the beginning of the series
  16 P.FirstScan = 1;
  17 P.TickTime  = 1000/P.Hz;
  18 [Y,T,P] = read_ppg(ppg_file,P);
  19 
  20 % get the start and end of scanning scanning
  21 % start - first series folder and first dicom
  22 % end   - last series fold and last dicom
  23 
  24 % get series folders
  25 ser_folder = dic_folder(1:strfind(dic_folder, '/Series'));
  26 dce_f = dir(ser_folder);
  27 dc1_folder = fullfile(ser_folder,dce_f(3).name);
  28 dce_folder = fullfile(ser_folder,dce_f(end).name);
  29 
  30 dc1_files = fullfl(dc1_folder,'1.3','dcm');
  31 dce_files = fullfl(dce_folder,'1.3','dcm');
  32 dc1_file  = dc1_files{1};
  33 dce_file  = dce_files{end};
  34 dc1_hdr   = spm_dicom_headers(dc1_file, false);
  35 dce_hdr   = spm_dicom_headers(dce_file, false);
  36 mri_start = dc1_hdr{1}.AcquisitionTime;
  37 mri_stop  = dce_hdr{1}.ContentTime;
  38 mri_dur   = mri_stop - mri_start;
  39 disp(['--- Total MRI time ---']);
  40 disp(['MRI start        : ' secs2hms(mri_start)]);
  41 disp(['MRI stop         : ' secs2hms(mri_stop)]);
  42 disp(['MRI duration     : ' secs2hms(mri_dur)]);
  43 
  44 % load dicom headers
  45 dicoms = fullfl(dic_folder,'1.3','dcm');
  46 dicoms = dicoms(P.FirstScan:end);
  47 d_1st = 1;
  48 d_end = size(dicoms,1);
  49 hdr_1 = spm_dicom_headers(dicoms{d_1st}, false);
  50 hdr_e = spm_dicom_headers(dicoms{d_end}, false);
  51 
  52 % dicom start (in seconds since midnight)
  53 d_start = hdr_1{1}.AcquisitionTime;
  54 d_stop  = hdr_e{1}.ContentTime;
  55 
  56 % ppg start
  57 p_start = P.MRIStartTime/1000; % secs
  58 p_stop  = P.MRIStopTime/1000;  % secs
  59 p_dur   = p_stop - p_start;    % secs
  60 offset  = d_start - p_start;   % secs
  61 i_s = ceil(offset*P.Hz);       % secs * Hz = indices
  62 i_e = ceil((offset+(d_end*P.TR/1000))*P.Hz);
  63 %i_e = ceil((d_stop - p_start)*P.Hz);
  64 
  65 disp(['Phys log start   : ' secs2hms(p_start)]);
  66 disp(['Phys log stop    : ' secs2hms(p_stop)]);
  67 disp(['Phys log duration: ' secs2hms(p_dur)]);
  68 disp(['Freq estimate    : ' num2str(size(Y,2)/p_dur) 'hz, Freq used: ' num2str(P.Hz) 'hz'])
  69 disp(['--- Run information ---']);
  70 disp(['Number of dicoms : ' num2str(d_end)]);
  71 disp(['First dicom      : ' secs2hms(d_start)]);
  72 disp(['Last dicom       : ' secs2hms(d_stop)]);
  73 disp(['Phys-dicom       : ' secs2hms(d_start-p_start)]);
  74 disp(['Phys indices     : ' num2str(i_s) ' ... ' num2str(i_e)]);
  75 
  76 % extract ppg data pertaining to this run
  77 Y = Y(i_s:i_e);
  78 T = T(i_s:i_e);
  79 T = T - T(1);
  80 disp(['Run duration from dicoms      : ' num2str((size(Y,2).*P.TickTime)/1000/60) ' mins']);
  81 disp(['Phys. data for this run       : ' num2str((T(end)-T(1))/60) ' mins']);
  82 
  83 % check if there is any data in Y
  84 if numel(unique(Y))==1
  85     disp('-- ERROR --');
  86     disp('The physiological signal extracted for this time period is completely flat!');
  87     disp(['Check the file ' ppg_file]);
  88     disp('for the time period reported above');
  89     disp('Noise regressor cannot be estimated.');
  90     disp('This might be because your recording equipment was detached during this scanning period');
  91     error('Exiting ..');    
  92 end
  93 
  94 % saving
  95 [fpath fname fext] = fileparts(ppg_file);
  96 save(fullfile(fpath,[fname '.mat']),'Y','T','P');
  97 
  98 function [signal,time,P] = read_ppg(pulsefile,P)
  99 
 100 fclose('all');
 101 
 102 fid=fopen(pulsefile);
 103 ignore=textscan(fid,'%s',4); %Ignore first 4 values.
 104 
 105 data   = textscan(fid,'%u16'); %Read data until end of u16 data.
 106 footer = textscan(fid,'%s');   %Read in remaining data (time stamps and statistics).
 107 
 108 if ~size(footer{1},1)
 109     disp('No footer found in the log file!');
 110     disp('Hence, the start of the logging is unknown');
 111     disp('and match with dicoms cannot be calculated.');
 112     disp('This might indicate that the equipment was switched off unexpectedly.');
 113     error('Exiting..');
 114 end
 115 
 116 %Get time stamps from footer:
 117 for n=1:size(footer{1},1)
 118     if strcmp(footer{1}(n),'LogStartMDHTime:')  %log start time
 119         P.LogStartTime=str2num(footer{1}{n+1});
 120     end
 121     if strcmp(footer{1}(n),'LogStopMDHTime:')   %log stop time
 122         P.LogStopTime=str2num(footer{1}{n+1});
 123     end
 124     if strcmp(footer{1}(n),'LogStartMPCUTime:') %scan start time
 125         P.MRIStartTime=str2num(footer{1}{n+1});
 126     end
 127     if strcmp(footer{1}(n),'LogStopMPCUTime:')  %scan stop time
 128         P.MRIStopTime=str2num(footer{1}{n+1});
 129     end
 130 end
 131 
 132 % Remove the systems own evaluation of triggers.
 133 t_on  = find(data{1} == 5000);  % System uses identifier 5000 as trigger ON
 134 t_off = find(data{1} == 5003);  % System uses identifier 5003 as trigger OFF
 135 
 136 % Filter the trigger markers from the data
 137 Hz = P.Hz;
 138 data_t=transpose(1:length(data{1}));
 139 indx = setdiff(data_t(:),union(t_on,t_off)); %Note: depending on when the scan ends, the last size(t_off)~=size(t_on).
 140 signal = data{1}(indx);
 141 time = (1:length(signal))./Hz;
 142 
 143 % Clip the time series to begin and end with the scan.
 144 if P.TR < 1000
 145   start = P.StartClip*(Hz/1000);
 146   stop  = length(signal) - Hz*floor([length(signal) - EndClip*(Hz/1000)]/Hz);
 147 else
 148   start = P.StartClip*(Hz/1000)+Hz/2;
 149   stop  = P.EndClip*(Hz/1000)+mod(length(signal(1:end)),Hz)+ Hz/2;
 150   %stop  = length(signal) - Hz*floor([length(signal) - P.EndClip*(Hz/1000)]/Hz) + Hz/2;
 151 end;
 152 
 153 % Reset vectors;
 154 signal = double(signal(start+1-P.StartJump:end-stop-P.StartJump)');
 155 time = time(start+1-P.StartJump:end-stop-P.StartJump);
 156 time = time(:)-time(1);
 157 
 158 function [Y,T,P] = read_pulsef(pulsefile)
 159 
 160 [fid, emsg]  = fopen(pulsefile);
 161 
 162 if fid==-1
 163     error(emsg);
 164 end
 165 
 166 tline = fgetl(fid);
 167 C = [];
 168 
 169 while ischar(tline)
 170     C{end+1,1} = tline;
 171     tline = fgetl(fid);
 172 end
 173 fclose(fid);
 174 
 175 Y = str2num(C{1,1});
 176 p = Y(1:5); % first 5 values are parameters
 177 Y = Y(6:end); % resize Y, starting with the 6th value of the first line are the voltage values themselves that reflect the pulse ox.
 178 
 179 P.TickTime   = p(3); %time (in msec) between subsequent measurements
 180 P.StartMDH   = get_timeval(C{11,1});
 181 P.StopMDH    = get_timeval(C{12,1});
 182 P.StartMPCU  = get_timeval(C{13,1});
 183 P.StopMPCU   = get_timeval(C{14,1});
 184 
 185 function timeval = get_timeval(ccell)
 186 [~, timeval] = strtok(ccell);
 187 timeval = str2num(rmwhite(timeval));