%	to open physiological data file and do some plotting of
%	of the ecg and respiratory data.  
%	
%	rev 0 	11/3/99		original from motplot
%	rev 1 	12/7/99		does the resp phase from peaks
%	rev 2 	12/14/99	calculates running HR, resp rate, and
%				breath volume
%				makes an output for each acquisition
%	rev 3 	1/25/00		makes an output every desired TR
%	rev 11 	2/21/04		to read file from excite physio
%	rev 12 	2/25/04		file has trig time pts and then resp waveform
%				trig table is terminated by -9999 entry
%		2/14/05		changed the resp threshold and made the
%				output in resp rate rather then interval
%		7/26/05		fix odd factor of two in the 'ignore time'
%		9/3/07		add rvt plot, a la R. Birn.
%	rev 20  5/1/08		20x samples PG at 10 ms, resp at 40ms
%		5/25/08		use the initial 30s to set up the convolver
%		8/4/08          offset time start as Tscan - nframes*TR
%                               because start of record not tied to scan start

clear; 
kern2 = [19:-1:-19];

PRESC = 30.0;		% seconds of prepended data
dtr = 0.040;		% in seconds
dth = 0.010;		% in seconds
Twin = 6;		% integration window

%  open the input file

fname = input('gimme file = ', 's');
fid = fopen(fname, 'r');
[dat nf] = fscanf(fid, '%g\n');
fclose(fid);

foo = input('gimme [nframes TR(s)] = ');
nframes = foo(1);  TR = foo(2);
Tscan = nframes*TR;

% parse the triggers and resp waveform

ntrig = find(dat == -9999)-1;
etrig = dat(1:ntrig)*dth;	% trigger times
respr = dat(ntrig+2:end);	% resp waveform

nresp = length(respr);
ne = length(etrig);
aveecg = (etrig(ne) - etrig(1))/ne;
avehr = 60/aveecg;		% BPM
fprintf('average HR = %d BPM\n', fix(avehr));

ts = nresp*dtr;	% sampled time
time0 = ts - Tscan;		% how much to ignore
fprintf('start time = %.3f\n', time0 - PRESC);
if(time0 < 0)
  fprintf('Eek!  Not enough data- I give up!\n');
  return
end

Tout = input('output sample interval, s [2s]) = ');
if(isempty(Tout))
  Tout = 2;
end
nout = fix(Tscan/Tout);
time = (0:nout-1)*Tout;
fprintf('num output samples = %d\n', nout);

% find ecg intvl for each cardiac cycle and thence hrate

clear hrate;
for j=1:nout
  t = time0 + time(j);
  i1 = 1; i2 = 1;
  for i1=1:ntrig
    if(etrig(i1)>=t-Twin*.5)
      break;
    end
  end
  for i2=i1:ntrig
    if(etrig(i2)>t+Twin*.5)
      break;
    end
  end
  i2 = i2 - 1;
  if(i2 == i1)              % end of trace
    i1= i1 - 1;
  end
  hrate(j) = (i2-i1)*60/(etrig(i2) - etrig(i1));       % bpm 
end

subplot(2,2,1); 
plot(time, hrate);grid
ylabel('hrate, BPM');
xlabel('time, s');
title(fname);

%  now do the resp.  find the peaks

respx = max(respr);
respn = min(respr);
resp = 100*(respr - respn)/(respx - respn);
subplot(2,2,2); 
n1 = fix(time0/dtr);
plot((1:nresp-n1+1)*dtr, resp(n1:nresp));grid
ylabel('respiration'); 
xlabel('time, s');

drdt = conv(resp, kern2);
drdt = drdt(19:nresp+18);
d2rdt2 = conv(drdt, kern2);
d2rdt2 = d2rdt2(19:nresp+18);
rpeak = (d2rdt2 > 0.5e-6);  	% nice threshold

% find the resp trigs

nr = 0; 
for (j=2:nresp)
  if (rpeak(j)==1 & rpeak(j-1)==0)  % first only
    nr = nr + 1;
    rtrig(nr) = j*dtr;
  end
end
averesp = (rtrig(nr) - rtrig(1))/nr;
averrate = 60/averesp;		% breaths/min
fprintf('average resp = %d breaths/min\n', fix(averrate));

% find resp intvl for each breath

resptime = diff(rtrig);
nrespt = nr - 1;
tresptime = (rtrig(1:nrespt) + rtrig(2:nr))*.5;
respintvl = spline(tresptime, resptime, time+time0);
rrate = 60./respintvl;
subplot(2,2,3); 
plot(time, rrate);grid
ylabel('resp rate, BrPM');
xlabel('time, s');

% plot rv(t)

clear rv;
for j = 1:nout
  t = time0 + time(j);
  i1 = fix((t - Twin*.5)/dtr);
  i2 = min(nresp, fix(t + Twin*.5)/dtr);
  rv(j) = std(resp(i1:i2));
end
subplot(2,2,4); 
plot(time, rv);  grid
ylabel('RV')
xlabel('time, s');
  
% save the output

fnout = input('output file [cr for default, s for none]= ', 's');
if(isempty(fnout))
  fnout = sprintf('%s.out',fname);
  elseif(strcmp(fnout, 's'))
  return;
end
fout = fopen(fnout, 'w');
fprintf(fout, '%f  %f %f\n', [hrate', rrate', rv']');
fclose(fout); 
fprintf('wrote file  %s\n', fnout);