function ecg2hr(ihbv, sdt, waveq, stopt, ihbt, hrdt, hbwt, nswt)
%ECG2HR(ihbv, sdt, waveq, stopt, ihbt, hrdt, hbwt, nswt)
%
%  ihbv (.6): initial heart beat values; an average initial-heartbeat y-value; a filter
%  sdt (75): slope differential threshold; larger means more filter; a major heartbeat filter
%  waveq (3): wave quality of the ecg: use a higher number with a more wavy ecg
%  stopt (300) : stop execution after (stopt) seconds of input have passed
%  ihbt ([]): initial heart beat times: array of initial-heartbeat x-values, e.g., [.085 .143]
%  hrdt (50): output a heart rate measurement every (hrdt) milliseconds
%  hbwt (60): heart beat warning threshold; closer to zero means more warnings
%  nswt (10): noisy signal warning threshold; closer to zero means more warnings
%
%Examples:
%  ecg2hr
%  ecg2hr(20, 25, 18);

%From here on, a number in parenthesis will denote a time dependent number
%These numbers would likely change with a Sampling Frequency [fs] < 1000
%or if this ECG program were used with, say, humans rather than mice

% Tests on the input arguments
[ecgfilename, cdfilename, rrfilename, hrfilename, ecgfilenameonly] = chooser('useold');
if (nargin==0) ihbv = (.6); end %if
if (nargin<=1) sdt = (75); end %if
if (nargin<=2) waveq = (3); end %if
if (nargin<=3) stopt = (300); end %if
if (nargin<=4) ihbt = []; end %if
if (nargin<=5) hrdt = (50); end %if
if (nargin<=6) hbwt = 60; end %if
if (nargin<=7) nswt = 10; end %if
if (ihbv==0 | sdt==0 | waveq==0 | hrdt==0 | hbwt==0 | nswt == 0)
  error('Parameters ihbv, sdt, waveq, hrdt, hbwt, and nswt cannot equal 0; check your parameters.'); end %if

%used to convert between seconds and milliseconds
TIMEBASIS = .001; INVTIMEBASIS = 1000;
%heartbeat boundaries in ms
MINRINTERVAL = (66); MAXRINTERVAL = (400);
%set up running averages of the given length
AVGLENGTH = (200); SLOPERANGE = (2); RRAVGLENGTH = 5;
%used to read and store the ecg
t = (-inf); prevecg = 0; ecgread = 0; tims = 0;
%used for basic calculations on the ecg
avg = 0; avgdiv = 0; ema = 0; ecgmult = 1/ihbv;
%used to find the fs of the ecg
t0 = 0; tcount = 0; edt = 0; edtims = 0; fs = 1000;
%used to analyze the slope of the ecg
avgslope = 0; slopediv = 0; curslope = 0;
%used to analyze the change in slope of the ecg
slopediff = 0; prevecg = 0; prevt = 0; dt = (1);
%used to find rr intervals
recentrt = -inf; rrcount = -1; prevema = 0; ishb = 0;
%used for some rr calculations and output
avgrrint = 0; minrravg = inf; i = 0; rtoffset = 0;
%used to analyze the rr intervals
prevrr = 0; bestrt = 0; bestrr = 0; bestema = 0; pbema = 0;
%used to find best hb match
fdde = 0; hbmatchdiff = 0; minfdde = inf; minhbmatchdiff = inf;
%used to warn user about possibly bad output
prevnscount = 0; nscount = 0; wrt = 3.5; iwt = 0; pwt = 0;
%file ids
ecgfid = fopen(ecgfilename,'r');
cdfid = fopen(cdfilename,'w');
rrfid = fopen(rrfilename,'w');
hrfid = fopen(hrfilename,'w');

while (t < stopt & ~feof(ecgfid))
  %reading and analyzing the ecg
  prevt = t;
  prevecg = ecgread;
  ecgline = fgetl(ecgfid);
  if (isempty(ecgline))
    if (tcount == 0)
      if ~isempty(sscanf(ecgline,'%[# b-df-mo-zB-DF-MO-Z]'))
        error('Data from ECG input file is not valid.  Please run Splice if appropriate.');
        return;
      end %if
      while (isempty(ecgline))
        ecgline = fscanf(ecgfid,'%f %f',[1 2]);
      end %while
    else
      break;
    end %if
  end %if
  ecgline = sscanf(ecgline,'%f %f',[1 2]);
  t = ecgline(1) - t0; ecgread = ecgline(2)*ecgmult;
  if (t <= prevt) error('ECG time column must increase continuously.'); end %if
  if (tcount >= 0)
    if (tcount == 0)
      t0 = t;
      t = 0;
      prevt = 0;
      prevecg = ecgread;
    else
      dt = t - prevt;
      edt = edt + (dt - edt)/tcount;
    end %if
    tcount = tcount+1;
  end %if
  tims = (round(t*fs)/fs)*INVTIMEBASIS;
  prevema = ema;
  ema = ecgread - avg;
  avgdiv = avgdiv + (AVGLENGTH-avgdiv)/AVGLENGTH;
  avg = avg + ema/avgdiv;
  curslope = (ecgread - prevecg)/dt;
  slopediff = curslope - avgslope;
  if (slopediv < SLOPERANGE) slopediv = slopediv + 1; end %if
  avgslope = avgslope + slopediff/slopediv;
  temprrint = tims - recentrt;

  %is this ecg value a heartbeat?
  if (size(ihbt, 2) >= 2 & rrcount <= size(ihbt, 2) - 2)
    ishb = sum(t == ihbt);
  else
    ishb = rrcount > 0 | abs(slopediff*prevema) > ihbv*(125);
    ishb = abs(slopediff) > sdt & temprrint > MINRINTERVAL & (ishb);
  end %if
  if (ishb)
   %a new possible heartbeat has been found
    if (rrcount == -1)
     %first heartbeat, need one more to calculate heart rate
      rrcount = 0;
      recentrt = tims; %should be tims-edtims, but we don't know edtims yet (see below)
    elseif (rrcount == 0)
     %second heartbeat, now heart rate (and edtims) is known
      rrcount = 1;
      %calculate final fs, edt
      fs = round(.1/edt)*10;
      edt = 1/fs;
      edtims = INVTIMEBASIS/fs;
      tcount = -1;
      %output fs
      fprintf(cdfid, '%5g\n', fs);
      fprintf(rrfid, '%5g\n', fs);
      fprintf(hrfid, '%5g\n', fs);
      %set rr variables
      bestrt = tims-edtims;
      recentrt = recentrt-edtims; %-edtims to account for above discrepancy
      bestrr = temprrint; %(1-1)*edtims to account for recentrt discrepancy
      avgrrint = bestrr;
      if (bestrr > MAXRINTERVAL)
        fclose('all');
        error(['Heartbeat interval=' num2str(bestrr*TIMEBASIS) ' is too long.  Try with ihbv=' ...
            num2str(round(ihbv*.75)) '.']);
      end %if
      %output this rr interval      
      fprintf(rrfid, '%10g %10g\n', bestrt, bestrr);
      for i=(0:edtims:recentrt-edtims) fprintf(cdfid, '%11g %10g\n', i*TIMEBASIS, 0); end %for
      fprintf(cdfid, '%11g %10g\n', recentrt*TIMEBASIS, 1);
      for i=(recentrt+edtims:edtims:bestrt-edtims) fprintf(cdfid, '%11g %10g\n', i*TIMEBASIS, 0); end %for
      fprintf(cdfid, '%11g %10g\n', bestrt*TIMEBASIS, 1);
      for i=(recentrt:hrdt:bestrt)
        fprintf(hrfid, '%10g %10g\n', i, INVTIMEBASIS/bestrr);
      end %for
      rtoffset = mod(bestrt-recentrt+rtoffset-hrdt,hrdt);
      prevrr = bestrr;
      recentrt = bestrt;
    else
     %all other heartbeats execute these statements
      nscount = nscount + 1;
      if (rrcount <= size(ihbt, 2) - 2 | (temprrint > avgrrint * 1.5) | ...
          (temprrint > avgrrint * 1.4125 & nscount > 1))
       % past likely RR interval boundary, so calculate heart rate
        if (rrcount < RRAVGLENGTH) rrcount = rrcount + 1; end %if
        avgrrint = avgrrint + (bestrr-avgrrint)/rrcount;
        if (avgrrint < minrravg)
          minrravg = avgrrint;
        elseif (minrravg == 0)
        elseif (avgrrint > minrravg*2.125)
          disp('Please check both the output as well as your data.');
          disp(['Around or before t=' num2str(t) ', this program may']);
          disp('have skipped and continued to skip heartbeats.');
          disp('Quick fix: Fudge data to make heartbeats obvious.');
          minrravg = 0;
        end %if
        if (bestrr > MAXRINTERVAL)
          fclose('all');
          error(['Heartbeat interval=' num2str(bestrr*TIMEBASIS) ...
              ' is too long.  Adjust parameters ihbv and sdt (type ''help ecg2hr'' for more info).']);
        end %if
        if (minfdde > hbwt | (prevnscount > nswt & minfdde > hbwt/(1.1)))
         %list all possibly bad output points
          if (iwt == 0)
            disp('Ecg2hr generated questionable output during the following time ranges:');
            iwt = floor((recentrt*TIMEBASIS-.3)*2)/2;
            pwt = iwt;
          elseif (t-pwt < wrt)
            pwt = ceil((recentrt*TIMEBASIS)*2+.2)/2;
          else
            if (pwt == iwt)
              disp(num2str(iwt));
            else
              disp([num2str(iwt) '-' num2str(pwt)]);
            end %if
            iwt = floor((recentrt*TIMEBASIS-.3)*2)/2;
            pwt = iwt;
          end %if
        end %if
        %output rr interval
        fprintf(rrfid, '%10g %10g\n', bestrt, bestrr);
        for i=(recentrt+edtims:edtims:bestrt-edtims) fprintf(cdfid, '%11g %10g\n', i*TIMEBASIS, 0); end %for
        fprintf(cdfid, '%11g %10g\n', bestrt*TIMEBASIS, .5*bestema);
        for i=(recentrt-rtoffset+hrdt:hrdt:bestrt)
          fprintf(hrfid, '%10g %10g\n', i, INVTIMEBASIS/bestrr+...
            (i<=recentrt+hrdt)*(INVTIMEBASIS/prevrr-INVTIMEBASIS/bestrr)/2);
        end %for
        rtoffset = mod(bestrt-recentrt+rtoffset-hrdt,hrdt);
        prevrr = bestrr;
        recentrt = bestrt;
        %reset rr variables
        temprrint = tims - recentrt;
        pbema = bestema;
        bestema = 0;
        minhbmatchdiff = inf;
        prevnscount = nscount;
        nscount = 0;
        AVGLENGTH = avgrrint/waveq;
      end %if
      if (prevema~=0 & temprrint > MINRINTERVAL)
       %test heartbeat candidate to see if it's a good match
        fdde = abs((abs(INVTIMEBASIS/(temprrint-edtims)-INVTIMEBASIS/prevrr)+10)/prevema);
        if (prevnscount == 0) prevnscount = 1; end %if
        pnscdw = ceil(prevnscount/nswt);
        hbmatchdiff = (abs(INVTIMEBASIS/(temprrint-edtims)-INVTIMEBASIS/prevrr)+10)*...
          (-atan(pnscdw*nswt*abs(prevema)-.5)/pi+.5+.28/pnscdw)*...
          (-2*atan(prevnscount/nswt)/pi*(sign(prevema) == sign(pbema))+1);
        if (hbmatchdiff <= minhbmatchdiff)
          minfdde = fdde;
          minhbmatchdiff = hbmatchdiff;
          bestema = prevema;
          bestrt = tims-edtims;
          bestrr = temprrint-edtims;
        end %if
      end %if
    end %if
  end %if

end %while

if (iwt ~= 0)
 %output final questionable points
  if (pwt <= iwt)
    disp(num2str(iwt));
  else
    disp([num2str(iwt) '-' num2str(pwt)]);
  end %if
end %if

disp(['Finished at t = ' num2str(t)]); 

fclose(ecgfid);
fclose(hrfid);
fclose(rrfid);
fclose(cdfid);