function ecg2hr(stopt, hrdt)
%ECG2HR(stopt, hrdt)
%
%  stopt (300): stop execution after (stopt) seconds of input have passed
%  hrsf (20): heart rate sampling frequency; can be chosen arbitrarily
%
%Examples:
%  ecg2hr
%  ecg2hr(5)  %parse first 5 seconds of data

%From here on, a number in parenthesis will denote a time dependent number
%These numbers would likely change with a Sampling Frequency [sf] < 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) stopt = (300); end %if
if (nargin<=1) hrsf = (20); end %if
if (hrsf==0)
  error('Parameters hrsf, hbwt, and nswt cannot equal 0; check your parameters.'); end %if

%used to convert between seconds and milliseconds
TIMEBASIS = .001; INVTIMEBASIS = 1000;
%used to define heartbeat boundaries in ms
MINRRINTERVAL = (66); MAXRRINTERVAL = (500);
%used to set up running averages of the given length
AVGLENGTH = (20); SLOPERANGE = (2); RRAVGLENGTH = 5;
%used for heartbeat and noise running averages
HBTLENGTH = 5; LTTLENGTH = 5; NAVGLENGTH = 200;
%constants/thresholds used to determine heartbeats
HBTM = .3; HBTRTM = .8; NSM = 1.1; SDT = (20);
%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; waveq = 1;
%used to find the sampling frequency of the ecg
t0 = 0; tcount = 0; edt = 0; edtims = 0; sf = 1000;
%used to analyze the slope of the ecg
avgslope = 0; curslope = 0; aaslope = 0;
%used to analyze the change in slope of the ecg
slopediff = 0; prevecg = 0; prevt = 0; dt = (.001);
%used to define and utilize heartbeat threshold band
phbt = 0; nhbt = 0; ptrtims = 0; lttcount = 0; pltttims = 0;
%used to find rr intervals
rt0 = 0; recentrt = 0; prevema = 0; ishb = 0;
%used for rr calculations and output
avgrrint = 0; i = 0; rrstate = -1; rrcount = 0;
%used to analyze the rr intervals
prevrr = 0; bestrt = 0; bestrr = 0; bestema = 0; pbema = 0;
%used to find best hb match
hbmatchdiff = 0; minhbmatchdiff = inf;
%used to determine the amount of noise
navg = 0; bestnavg = 0; aflag = 0;

%file ids
ecgfid = fopen(ecgfilename,'r');
cdfid = fopen(cdfilename,'w');
rrfid = fopen(rrfilename,'w');
hrfid = fopen(hrfilename,'w');

while (t < stopt & ~feof(ecgfid))
  %input a line form the ecg file
  ecgline = fgetl(ecgfid);
  if (tcount == 0)
    if (isempty(ecgline))
      while (isempty(ecgline))
        ecgline = fgetl(ecgfid);
      end %while
    end %if
    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.');
    end %if
  end %if
  %read and analyze the ecg
  ecgline = sscanf(ecgline,'%f %f',[1 2]);
  prevt = t;
  prevecg = ecgread;
  t = ecgline(1) - t0; ecgread = ecgline(2);
  if (t <= prevt) error('ECG time column must increase continuously.'); end %if
  if (tcount >= 0)
   %start the process of figuring out the ECG sampling frequency
    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
  %figure out parameters used to determine heartbeats
  tims = (round(t*sf)/sf)*INVTIMEBASIS;
  prevema = ema;
  ema = ecgread - avg;
  if (ema > phbt*HBTRTM/HBTM)
    phbt = phbt + HBTM*ema/HBTLENGTH;
    ptrtims = tims;
  elseif (ema > 0)
    phbt = phbt*.999;
  end %if
  if (ema < nhbt*HBTRTM/HBTM)
    nhbt = nhbt + HBTM*ema/HBTLENGTH;
    ptrtims = tims;
  elseif (ema < 0)
    nhbt = nhbt*.999;
  end %if
  if (ema > phbt | ema < nhbt)
    lttcount = lttcount + 1;
    pltttims = tims;
  end %if
fprintf(hrfid, '%10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n', t, phbt+avg, nhbt+avg, abs(ema)+.2, navg+.2, (ishb + aflag)*.1+.6);
  if (tims - pltttims > MAXRRINTERVAL & rrstate ~= -1)
    rrstate = -1;
    disp(sprintf('Analysis stopped at t=%g since the value at t=%g was too large.', t, ptrtims*TIMEBASIS));
    rt0 = recentrt;
  end %if
  avgdiv = avgdiv + (AVGLENGTH-avgdiv)/AVGLENGTH;
  avg = avg + ema/avgdiv;
  curslope = (ecgread - prevecg)/dt;
  slopediff = curslope - avgslope;
  avgslope = avgslope + slopediff/SLOPERANGE;
  aaslope = aaslope + abs(slopediff)/(SLOPERANGE);
  navg = navg + (abs(ema) - navg)/NAVGLENGTH;
  temprrint = tims - recentrt;
  aflag = 0;

  %is this ecg value a heartbeat?
  ishb = (rrstate >= 0 | temprrint > MAXRRINTERVAL) & (lttcount);
  ishb = abs(slopediff) > SDT & temprrint > MINRRINTERVAL & (ishb);
  if (ishb)
   %a new possible heartbeat has been found
    if (rrstate == -1)
     %first heartbeat, need one more to calculate heart rate
      rrstate = 0;
      lttcount = 0;
      recentrt = tims; %should be tims-edtims, but possibly we don't know edtims yet
    elseif (rrstate == 0)
     %second heartbeat, now heart rate (and edtims) is known
      rrstate = 1;
      if (tcount >= 0)
       %calculate final sf, edt
        sf = round(.1/edt)*10;
        edt = 1/sf;
        edtims = INVTIMEBASIS/sf;
        %output sf
        fprintf(cdfid, '%5g\n', sf);
        fprintf(rrfid, '%5g\n', sf);
        tcount = -1;
      end %if
      %set rr variables
      bestrt = tims-edtims;
      recentrt = recentrt-edtims; %-edtims to account   above discrepancy
      bestrr = temprrint; %+(1-1)*edtims to account for recentrt discrepancy
      avgrrint = bestrr;
      lttcount = 0;
      if (bestrr > MAXRRINTERVAL)
        rrstate = -1;
        disp(sprintf('Analysis stopped at t=%g since the next heartbeat was too far away.', recentrt*TIMEBASIS));
        rt0 = tims;
      else
        %output this rr interval      
        disp(sprintf('Heartbeat found at t=%g.', (recentrt-edtims)*TIMEBASIS));
        rrcount = rrcount + 2;
        fprintf(rrfid, '%10g %10g\n', bestrt, bestrr);
        for i=(rt0:edtims:recentrt-edtims) fprintf(cdfid, '%11.3f %10.3f\n', i*TIMEBASIS, 0); end %for
        fprintf(cdfid, '%11.3f %10.3f\n', recentrt*TIMEBASIS, 1);
        for i=(recentrt+edtims:edtims:bestrt-edtims) fprintf(cdfid, '%11.3f %10.3f\n', i*TIMEBASIS, 0); end %for
        fprintf(cdfid, '%11.3f %10.3f\n', bestrt*TIMEBASIS, 1);
        prevrr = bestrr;
        recentrt = bestrt;
      end %if
    else
     %all other heartbeats execute these statements
      if (bestrr > MAXRRINTERVAL)
        rrstate = -1;
        disp(sprintf('Analysis stopped at t=%g since the next heartbeat was too far away.', recentrt*TIMEBASIS));
        rt0 = tims;
      end %if
      if ((temprrint > avgrrint * 1.5) & rrstate > 0)
       % past likely RR interval boundary, so calculate heart rate
        if (abs(bestema)*NSM < bestnavg)
         %best heartbeat was not discernable
            rrstate = -1;
            disp(sprintf('Analysis stopped at t=%g since there was too much intra-beat noise.', (tims-edtims)*TIMEBASIS));
            rt0 = recentrt;
        else
         %best heartbeat was discernable
          rrcount = rrcount + 1;
          if (rrstate < RRAVGLENGTH) rrstate = rrstate + 1; end %if
          avgrrint = avgrrint + (bestrr-avgrrint)/rrstate;
          %output rr interval
          fprintf(rrfid, '%10g %10g\n', bestrt, bestrr);
          for i=(recentrt+edtims:edtims:bestrt-edtims) fprintf(cdfid, '%11.3f %10.3f\n', i*TIMEBASIS, 0); end %for
          fprintf(cdfid, '%11.3f %10.3f\n', bestrt*TIMEBASIS, .8); %.5*bestema);
          prevrr = bestrr;
          recentrt = bestrt;
          pbema = bestema;
          %reset rr variables
          temprrint = tims - recentrt;
          bestema = 0;
          minhbmatchdiff = inf;
          waveq = 4; %waveq + .01*(lttcount-waveq)/LTTLENGTH;
          AVGLENGTH = avgrrint/waveq;
          lttcount = 0;
        end %if
      end %if
      if (prevema~=0 & temprrint > MINRRINTERVAL & rrstate > 0)
       %test heartbeat candidate to see if it's a good match
        hbmatchdiff = abs((abs(INVTIMEBASIS/(temprrint-edtims)-INVTIMEBASIS/prevrr)+10)/prevema)*...
          ((sign(prevema) ~= sign(pbema))+.1);
        if (hbmatchdiff <= minhbmatchdiff)
          minhbmatchdiff = hbmatchdiff;
          bestema = prevema;
          bestrt = tims-edtims;
          bestrr = temprrint-edtims;
          bestnavg = navg;
        end %if
      end %if
    end %if
  end %if

end %while

disp(sprintf('Finished at t=%g with %g heartbeats', t, rrcount)); 

fclose(ecgfid);
fclose(cdfid);
fclose(rrfid);
%fprintf(hrfid, '%5g %10g\n', sf, hrsf);
%hrt = 0:t
%[rt rr] = textread(rrfilename, '%f %f', -1, 'headerlines', 1);
%[
fclose(hrfid);