#!/usr/bin/octave

# theargs = char(argv);
arg_list = argv();
# wavename = theargs(1,:);
wavename = arg_list{1};
wavename
# pdfname = theargs(2,:);
pdfname = arg_list{2};
pdfname

# clear *;
format short g;

# profile on; # <<<



# in matlab, use flexible audioread instead of wav-only wavread
[wave, rate, bits] = wavread(wavename);

# show sampling rate and duration
rate
length(wave) / rate



# pdfname = wavename + '-stats.pdf';

# rough offset of some interesting irregularity, if any
interesting = 1;

# lowest index to consider for interesting beat
# must be at least 2 because of look back settings
# not be too low to not call early measurement artefacts interesting
mininterestingbeat = 25;

# enforce mono: if wave signal is for example stereo, use only first channel
# channels are columns: in other words, a mono wave is a vertical vector
wave = wave(1:end, 1);
# make wave a horizontal vector
wave = wave';
# size(wave)



printf('Notch filter and sampling rate conversion...\n');
# Remove 50 Hz hum with a very narrow notch filter before downsampling
# http://octave.sourceforge.net/signal/function/pei_tseng_notch.html
# frequencies are given as fraction of the Nyquist frequency (0..1)
[notch_b, notch_a] = pei_tseng_notch ( 2 * 50 / rate, 2 * 2 / rate );
wave = filter(notch_b, notch_a, wave);

# Also remove some popular higher harmonics of 50 Hz...
[notch_b, notch_a] = pei_tseng_notch ( 2 * 100 / rate, 2 * 2 / rate );
wave = filter(notch_b, notch_a, wave);
[notch_b, notch_a] = pei_tseng_notch ( 2 * 150 / rate, 2 * 2 / rate );
wave = filter(notch_b, notch_a, wave);
# Odd harmonics are usually stronger than even ones, zap some more of those:
[notch_b, notch_a] = pei_tseng_notch ( 2 * 250 / rate, 2 * 2 / rate );
wave = filter(notch_b, notch_a, wave);
[notch_b, notch_a] = pei_tseng_notch ( 2 * 350 / rate, 2 * 2 / rate );
wave = filter(notch_b, notch_a, wave);

# drop rate in ONE fast large step
if (1000 == gcd(1000,rate)) # if rate is exact multiple of 1000 samples/sec
# Default IIR 8th order Chebyshew type I, or 30 point FIR downsample:
# (fastest is "downsample": skips most samples, but does not low-pass first)
wave = decimate(wave, rate/1000); # decimate by int factor with anti-aliasing
else
# use fancy but slow polyphase resample by arbitrary quotient of integers:
wave = resample(wave, 1000/gcd(1000,rate), rate/gcd(1000,rate));
endif
rate = 1000;



# values in wave vectors have range ]-1;+1[ at most, normalize volume anyway
wave = wave .* ( 0.99 ./ max( max(wave), -min(wave) ) );
printf('Sampling rate converted to 1000/sec\n');

# Audacity exports space separated, so csvread would be unhappy
# beats = csvread('beatlist.csv'); # use row,column index!
# beats = dlmread('beatlist.txt',' ');

# only keep first column of beats
# beats = beats(:, 1);
# make that a horizontal array
# beats = beats';

# size of moving averages for heart rate and RR durations
mavgsize = 10;



# added 2016-07-16: inverse high pass to un-do soundcard distortion

# omega limit = 1 / R * C = 20 for circa 3 Hz (maybe 1 yF, 50 kOhm)
# goal is to undo the soundcard RC high pass filter with "anti high pass"!
# restored signal = 20 * integral_until_now + 1 * value_now
# use weight_int=20 and weight_now=1 for 3 Hz if dcdecay=0

# empirically, 8 to 12 seems to match soundcard filter better if dcdecay=1:
# Both 10 and 20 have good white noise frequency response but 20 has better
# square & saw reconstruction; 10 has better sine sweep frequency response,
# but soundcard OUTPUT for generating test signals may have high-pass, too!
# Without inverse pass: 12 dB/oct instead of 6 dB/oct filter, at << 3 Hz?

weight_int = 20;
weight_now = 1;
# disadvantage of high-pass removal: DC drift becomes visible again
# DC build-up prevention: use for example 1 for typical recordings
# use at most 3 to 5 to not bypass most of the high-pass removal.
dcdecay = 0;
# (adjust to make T to next P gap already flat and T not biphasic yet)
# if dcdecay 3 or 1, weight * 2 or 1.25 compensates a bit (check: square)
weight_int = weight_int * 1;

filterstate = 0;
invrate = 1.0 / rate;
filterstart = floor(0.01 * rate); # start applying filter at 0.01 sec
filterkeep = 1 - (dcdecay * invrate);
printf('Inverse low pass filter...\n');

if (dcdecay == 0)

# more efficient loop with only two additions and one multiplication inside
wave = wave * weight_now;
weight_all = weight_int * invrate / weight_now;

# for i = filterstart:length(wave)
#   # integrate whole signal over time
#   filterstate = filterstate + wave(i);
#   # now calculate the inverse of the high pass filter effect :-)
#   wave(i) = wave(i) + (weight_all * filterstate);
# endfor # filter

# Thanks to Gernot for this idea! cumsum(wave)(i) = sum of wave(1:i) ...
wave = wave + (weight_all * cumsum(wave));

else

# more flexible loop with fade out of the integrated value over time
for i = filterstart:length(wave)

  # integrate whole signal over time
  filterstate = filterstate + (wave(i) * invrate);
  # extra step against DC build-up: side effect is a high pass again...
  filterstate = filterstate * filterkeep;
  # now calculate the inverse of the high pass filter effect :-)
  wave(i) = weight_int * filterstate + weight_now * wave(i);

endfor # filter

endif

printf('Inverse low pass done!\n');
# end of code added 2016-07-16

# rotate to same orientation as wave, show amplitude:
peaksize = max( max(wave), -min(wave) )

# do NOT normalize again, it would shrink signals with DC drift!
# peaksize = 1.0;
# wave = wave .* ( 0.99 ./ peaksize );



# Savitzky Golay Filter: remove squiggly noise while preserving shape :-)
# Interesting paper: "What IS a Savitzky Golay Filter?" (zero phase, flat
# passband gain, moderate stopband attenuation, regular stopband zeros...)
# http://www-inst.eecs.berkeley.edu/~ee123/fa11/docs/SGFilter.pdf

# Defaults: 1st number order (3), 2nd number length (order+2 or 3, odd)
# For 40+ Hz, 1 ksps: 2,27 3,35 4,43 5,51 6,59 etc. (see ee123)
# For 150 Hz: 1,7 2,9 3,11 (4,13) 5,15 etc.

wave = sgolayfilt(wave, 2, 9); # postpone stronger filters to WAV output...
# wave = sgolayfilt(wave, 4, 11); # (as long as QRS finder can handle noise)



# *** set to 0 to bypass all ECG specific processing, or to 1 to enable it
isECG = 1;

# start of ECG specific processing
if (isECG)

# IMPORTANT: Settings are for Einthoven I channel, but II and III might
# still work with those without readjustments IF low pass filtered data
# For Einthoven III, qrsthreshold has to be a lot lower.

# characteristic R S peak distance on the time axis, for example 25 msec
findqrs = ceil(0.025 * rate);

# TODO different weight for positive and negative deflections, to reduce
# false positives coming from strong T waves?

# something between 0.3 and 1.3 should be appropriate here, for example 0.5
# *** 0.6 for good recordings to avoid T, lower to 0.5 to miss fewer QRS
qrsthreshold = 0.5

# minimal distance between two QRS beginnings, for example at least 0.2 sec
# *** using 0.25 sec here is a bit fishy, but it helps with noisy recordings
qrsmindistance = 0.25
# maximal distance between two QRS beginnings, for example 1.5 sec
qrsmaxdistance = 1.5

# assume that in each QRS complex, there is a time span where R comes
# roughly findqrs seconds before S: the amplitude will change a lot then
evidence = sign(max( (wave(1:end+1-findqrs) - wave(findqrs:end)) - qrsthreshold, 0));

# low pass filter for evidence vector: only believe stuff if persistent
evidenceLPFsize = ceil(0.008 * rate);
evidenceLPFweights = ones(1,evidenceLPFsize) / evidenceLPFsize;
evidence = filter(evidenceLPFweights, 1, evidence);

nbeats = 0;
inbeat = false;
invrate = 1.0 / rate;
prevbeat = -1;
prevrr = 0;
printf('QRS detector running...\n');

for i = 1:length(evidence)
  if (inbeat)
    if (evidence(i) < 0.1)
      inbeat = false;
    # else no action, inside beat and beat goes on
    endif
  elseif (evidence(i) > 0.9)
    currentrr = (i * invrate) - prevbeat;

    if (currentrr < qrsmindistance)
      # skip if too early
      warning("QRS too early absolute");
      if (interesting < 2 && nbeats > mininterestingbeat)
        interesting = round(beats(nbeats-2) * rate);
      endif
      currentrr
      i * invrate
      # skip over the whole current peak
      inbeat = true;
    else

      # could skip early beat here
      if (currentrr < (3 * prevrr / 5) && prevrr > 0)
        warning("QRS suspiciously early relative");
        if (interesting < 2 && nbeats > mininterestingbeat)
          interesting = round(beats(nbeats-2) * rate);
        endif
        currentrr
        prevrr
        i * invrate
      endif

      # could try to fill gap if too late
      if (currentrr > qrsmaxdistance && prevrr > 0)
        warning("QRS suspiciously late absolute");
        if (interesting < 2 && nbeats > mininterestingbeat)
          interesting = round(beats(nbeats-2) * rate);
        endif
        currentrr
        i * invrate
	if (currentrr > 2 * qrsmaxdistance)
          warning('ERROR: QRS twice as late as plausible')
          currentrr
          i * invrate
          # error('QRS twice as late as plausible, giving up')
        endif
      endif

      # could try to fill gap if too late
      if (currentrr > (5 * prevrr / 3) && prevrr > 0)
        warning("QRS suspiciously late relative");
        if (interesting < 2 && nbeats > mininterestingbeat)
          interesting = round(beats(nbeats-2) * rate);
        endif
        currentrr
        prevrr
        i * invrate
      endif

      inbeat = true;
      elapsed = 0;
      nbeats = nbeats + 1;
      # use some correction offset based on evidenceLPFsize and findqrs here?
      beats(nbeats) = i * invrate;
      prevbeat = beats(nbeats);

      # average over 4 previous RR periods, if so many are available
      if (nbeats > 4)
        prevrr = (beats(nbeats) - beats(nbeats - 4)) / 4;
      else
        prevrr = 0;
      endif

    endif
  endif
  # else outside beat and no new beat yet
endfor # evidence beat finder state machine loop



# for lower pulse rates, look back 0.35 sec from qRs, else 0.3 sec or less
# maybe automatically pick min(RR) - 0.2 sec or something like that here?
# *** 0.25 is good for normal pulse rates but 0.2 is better for fast pulse
theshift = 0.22

# must be able to look back by theshift sec from QRS
# ... otherwise skip first beat(s) in analysis ...
while ((beats(1) - (theshift + 0.02)) * rate < 2)
  skippedbeat = beats(1);
  # show value
  skippedbeat
  beats = beats(2:length(beats));
endwhile

# longest RR period
maxbig = 1;

# check from QRS - theshift sec to next QRS - theshift sec for each beat
prevt = 1;
nextt = (int32)((beats(1) - theshift) * rate);
for i = 2:length(beats)
  prevt = nextt;
  nextt = (int32)((beats(i) - theshift) * rate);
  big = 1 + nextt - prevt;
  if (big > maxbig)
    maxbig = big;
  endif
endfor # find longest possible beat



nbeats = length(beats);

printf('QRS detector done!\n');
# SHOW number of beats and max duration of a beat
nbeats
maxbig

# 2016-08-15 interpolate baseline via 40 msec snippets AFTER theshift & at edges
fortymsec = (int32)(rate * 0.02);
baselinetimes(1) = -0.001;
baselimepoints(1) = mean(wave(1:fortymsec));
for i = 1:nbeats
  # TODO: reduce shift if previous beat is too nearby due to high pulse rate
  prevt = (int32)((beats(i) - theshift) * rate); # marks R peaks
  nextt = (int32)((beats(i) - theshift) * rate) + fortymsec; # marks R peaks
  if prevt < 2
    prevt = 2; # slightly after point 1
  endif
  if nextt < prevt
    nextt = prevt + fortymsec;
  endif
  baselinepoints(i+1) = mean(wave(prevt:nextt));
  # baselinetimes must not be int32, otherwise interpolation goes wrong...
  baselinetimes(i+1) = (double)(prevt + nextt) * 0.5; # time unit: samples!
endfor # collect baseline samples  
baselinetimes(nbeats+1+1) = length(wave) + 0.001;
baselinepoints(nbeats+1+1) = mean(wave(end - fortymsec : end));

bigtimes = linspace(0, length(wave), length(wave)); # time unit: samples!
# Piecewise Cubic Hermite InterPolation: smooth but follows curve closely
baselinewave = interp1(baselinetimes, baselinepoints, bigtimes, 'pchip');
# Linear interpolation between area directly after QRS-theshift of each beat:
# baselinewave = interp1(baselinetimes, baselinepoints, bigtimes, 'linear');
newwave = wave .- baselinewave;

# for i=1:length(newwave)
#   if isnan(newwave(i)) # failed to interpolate, outside range?
#     newwave(i) = wave(i);
#     error('Interpolation failure');
#   endif
# endfor # find missing samples
# plot(baselinetimes,baselinepoints, '.', bigtimes,wave, 'k', bigtimes,newwave, 'g', bigtimes,baselinewave, 'c');
# waitforbuttonpress;

clear bigtimes, baselinewave;

# create average etc. of all beats as horizontal vectors
thesum = zeros(1, maxbig);
themin = zeros(1, maxbig) + 1.01;
themax = zeros(1, maxbig) - 1.01;

prevt = 1;
nextt = (int32)((beats(1) - theshift) * rate); # marks R peaks
for i = 2:nbeats
  prevt = nextt;
  nextt = (int32)((beats(i) - theshift) * rate);
  # empirical baseline measurement point: a bit unstable
  # beatbaseline(i) = wave((int32)((beats(i) - 0.05) * rate)); # empirical
  singlepulse = wave(prevt:nextt);
  beatmean(i) = mean(singlepulse); # or use median
  # for average, use mean after baseline correction
  singlepulse = newwave(prevt:nextt);
  singlemean = mean(singlepulse);
  singlepulse = singlepulse - singlemean; # or use median
  singlesize = length(singlepulse);
  thesum(1:singlesize) = thesum(1:singlesize) .+ singlepulse;
  # TODO create heat map instead of only showing absolute min and max
  themin(1:singlesize) = min(themin(1:singlesize), singlepulse);
  themax(1:singlesize) = max(themax(1:singlesize), singlepulse);
endfor # calculate min, max and average pulse, gather baseline curve

# only use baseline corrected waveform after this point
# end of 2016-08-15 baseline correction interpolation changes
wave = newwave;
clear newwave;

# normalize average pulse shape
thesum = thesum * (1/nbeats);
invrate = 1.0 / rate;

prevt = 1;
nextt = (int32)((beats(1) - theshift) * rate); # marks R peaks
for i = 2:nbeats
  prevt = nextt;
  nextt = (int32)((beats(i) - theshift) * rate);
  singlepulse = wave(prevt:nextt);
  singlenoise = singlepulse .- thesum(1:length(singlepulse));
  beatnoise(i) = sqrt(meansq(singlenoise));
endfor # track noise amounts over time



# show 700 to 1400 msec to keep time scale within "40 ms per grid div"
# alternatively use xlim to achieve the same effect...
if (length(thesum) > (1.4 * rate))
  thesum = thesum(1:round(1.4 * rate));
  themin = themin(1:round(1.4 * rate));
  themax = themax(1:round(1.4 * rate));
endif
if (length(thesum) < (0.7 * rate))
  thepad = ones(1,round(0.71 * rate) - length(thesum));
  thesum = horzcat(thesum, thepad * 0);
  themin = horzcat(themin, thepad * min(themin));
  themax = horzcat(themax, thepad * max(themax));
endif

timeaxis = 1:length(thesum);
timeaxis = 1000 .* ((invrate .* timeaxis) .- theshift);

myfig = figure;
orient(myfig, 'tall'); # orient(myfig, 'portrait');
# set(myfig, 'PaperOrientation', 'portrait', 'PaperUnits', 'centimeters');

subplot(4,2,1);
plot(timeaxis, thesum, ".k", timeaxis, themin, "-r", timeaxis, themax, "-r");
line([timeaxis(1) timeaxis(end)], [0 0]); # draw horizontal line: baseline
line([0 0], [min(themin) max(themax)]); # draw vertical line: R peak time
title('Gemittelte Pulsform');
xlabel('Zeit relativ zu R-Zacke [ms]');
ylabel('Spannung [unkalibriert]');
axis('tight');
# xlim([min(timeaxis) max(timeaxis)]);
# ylim([-1.6 1.6])
# legend('Average','Minimum amplitude','Maximum amplitude');
grid on;
grid minor on;

durations = beats(2:nbeats) - beats(1:nbeats-1);

subplot(4,2,2);
plot( 1000 .* sort(durations), ".k");
title('RR Intervalle nach Dauer sortiert');
xlabel(wavename);
ylabel('RR Abstand [ms]');
xlim([1 length(durations)]);
# ylim([ min(600, 1000 * min(durations) - 100) max(1000, 1000 * max(durations) + 100) ]);
ylim([ (1000 * min(durations) - 100) (1000 * max(durations) + 100) ]);
grid on;
grid minor on;

avgweights = ones(1,mavgsize) / mavgsize;
avgdurations = filter(avgweights, 1, durations);

# early values are too low because filter assumes 0s before data
if (length(avgdurations) >= mavgsize)
  avgdurations(1:mavgsize) = ones(1,mavgsize) * avgdurations(mavgsize);
endif

subplot(4,1,2);
plot( beats(2:nbeats), 1000 .* durations, "+k", beats(2:nbeats), 1000 .* avgdurations, "-r" )
title('RR Abstand und gleitender Durchschnitt');
xlabel('Zeit [s]');
ylabel('RR Abstand [ms]');
xlim([min(beats) max(beats)]);
# ylim([ min(600, 1000 * min(durations) - 100) max(1000, 1000 * max(durations) + 100) ]);
ylim([ (1000 * min(durations) - 100) (1000 * max(durations) + 100) ]);
grid on;
grid minor on;

heartrates = 60 ./ durations;
avgweights = ones(1,mavgsize) / mavgsize;
avgrates = filter(avgweights, 1, heartrates);

# early values are too low because filter assumes 0s before data
if (length(avgrates) >= mavgsize)
  avgrates(1:mavgsize) = ones(1,mavgsize) * avgrates(mavgsize);
endif

subplot(4,1,3);
# Pulsfrequenz korreliert mit Atmung und Anspannung
plot( beats(2:nbeats), heartrates, ".k", beats(2:nbeats), avgrates, "-r" );
title('Pulsfrequenz und gleitender Durchschnitt');
xlabel('Zeit [s]');
ylabel('Pulsfrequenz [pro Minute]');
xlim([min(beats) max(beats)]);
# ylim([ min(50, min(heartrates) - 10) max(100, max(heartrates) + 10) ]);
ylim([ (min(heartrates) - 10) (max(heartrates) + 10) ]);
grid on;
grid minor on;

subplot(4,1,4);
# Baseline Drift korreliert mit Atmung und wandernden Elektrodenpotentialen
[ax h1 h2] = plotyy( beats(2:nbeats), beatmean(2:nbeats), beats(2:nbeats), beatnoise(2:nbeats) );
title('Baseline Drift: Durchschnitts-Spannung Intervall (Linie) und RMS Abweichung vom Mittel (Punkte)'); # jeweils R - theshift
xlabel('Zeit [s]');
line(ax(1), [beats(2) beats(nbeats)], [0 0]); # draw horizontal line
ylabel(ax(1), 'Spannung [unkalibriert]');
ylabel(ax(2), 'RMS Abweichung Pulsform');
xlim(ax(1), [min(beats) max(beats)]);
ylim(ax(1), [min(beatmean) max(beatmean)]);
# set(h1, 'Marker', '+');
set(h1, 'Color', 'k');
xlim(ax(2), [min(beats) max(beats)]);
ylim(ax(2), [min(beatnoise) max(beatnoise)]);
set(h2, 'Marker', '.');
set(h2, 'LineStyle', ':');
set(h2, 'Color', 'k');
grid on;
grid minor on;

print(pdfname);

# end of ECG specific (isECG) processing
endif



# Filtfilt: filter, time-reverse, filter: cancels out phase distortion :-)
# Note that this means that filtfilt doubles the effective filter order...

# Apply a strong low pass filter pair without phase distortion
# [butb buta] =  butter(3, 150 / (rate/2)); # 3rd order low pass
# wave = filtfilt(butb, buta, wave); # apply & ylppa

# Remove DC baseline offset drift by applying a high pass filter
# Hardware filters would use 0.05 Hz to keep phase distortion low,
# but we mostly have to care about 2/3 Hz being the lowest ECG part
# *** [butb buta] =  butter(3, 0.5 / (rate/2), "high"); # 3rd order high pass
# *** wave = filtfilt(butb, buta, wave); # apply & ylppa

# Rate is known to be 1000/sec at this point... Low pass to print 0 to 40 Hz
# https://www.wavemetrics.com/products/igorpro/dataanalysis/signalprocessing/smoothing.htm
# For 40+ Hz, 1 ksps: 2,27 3,35 4,43 5,51 6,59 etc. (see ee123)
wave = sgolayfilt(wave, 4, 43);

# FAST median filter: to omit spectral components above 40 Hz in PDF output
# http://octave.sourceforge.net/signal/function/medfilt1.html
# wave = medfilt1(wave,round(rate/(2*40))); # in our case: 13



# set to 1 to enable output of a WAV file with the processed data
outputRawData = 1;

# output raw data as WAV file
if (outputRawData)
# normalize amplitudes
peaksize = max( max(wave), -min(wave) );
bigwave = wave .* ( 0.99 ./ peaksize );
# bigwave = resample(bigwave, 8, 1); # polyphase resample by any quotient
# rate is known to be 1000 here, which makes the resample a lot faster:
bigwave = interp(bigwave, 8, 4); # upsample 8x using 2*8*4+1 FIR filter
wavwrite(bigwave, 8000, 16, 'rohdaten_ekg.wav');
clear bigwave;
endif



# simply set this to (0 or) 1 to (not) output a PDF with curve excerpts
showRawData = 1;

if (showRawData)

# offset of interesting events or 1: round down to previous full second
interesting = 1 + rate * floor(interesting/rate);
# if possible, show two extra seconds before the event
if (interesting > 2*rate)
interesting = interesting - 2*rate;
endif

# *** manual selection of interesting time offset
# interesting = 496 * rate;
# interesting = 1;

# number of raw data rows per page
stripes = 4;
# stripe size in seconds: 3 to 7 for 1 second as 1 major & 5 minor grid units
secperstripe = 7;

# figure; # landscape
myfig2 = figure;
orient(myfig2, 'tall'); # orient(myfig2, 'portrait');
# set(myfig2, 'PaperOrientation', 'portrait', 'PaperUnits', 'centimeters');

for stripe = 1:stripes

subplot(stripes,1,stripe); # several non-tall stripes on top of each other
stripeStart = interesting + round(secperstripe*rate*(stripe-1));
stripeEnd = interesting + round(secperstripe*rate*stripe) - 1;

if (stripeEnd > length(wave))
padwave = ones(1, stripeEnd - length(wave)) * wave(end);
wave = horzcat(wave, padwave);
clear padwave;
endif

timeaxis2 = stripeStart:stripeEnd;
timeaxis2 = timeaxis2 * (1000/rate);

# only plot if at least 1 second left to show
if (stripeEnd - stripeStart) > 1
plot(timeaxis2, wave(stripeStart:stripeEnd), "k", "linewidth", 2);
if (max(wave(stripeStart:stripeEnd)) < 1.1 && min(wave(stripeStart:stripeEnd)) > -1.1)
  ylim([-1.1 1.1]); # try to stick to a normalized default vertical scale
else
  axis('tight'); # allow bigger Y amplitude but stay reasonable
endif
if (stripe < 2)
  title(wavename);
endif
xlabel('Zeit [ms]');
ylabel('Spannung [unkalibriert]');
grid on;
grid minor on;
endif

endfor # loop over stripes

print('rohdaten_ekg.pdf');
endif



# simply set this to (0 or) 1 to (not) output a PDF with curve averages at
# QRS-theshift, size 0.8s, which represent 1/20 of the complete recording
showAvgData = 1;

if (showAvgData && isECG) # and maybe nbeats > nexcerpts, calculated below?

# number of average data rows per page
stripes = 4;
# number of seconds per average data row (3 to 7 for 1 sec per major grid)
secperstripe = 6.4;

# figure; # landscape
myfig2 = figure;
orient(myfig2, 'tall'); # orient(myfig2, 'portrait');
# set(myfig2, 'PaperOrientation', 'portrait', 'PaperUnits', 'centimeters');

for stripe = 1:stripes

subplot(stripes,1,stripe); # several non-tall stripes on top of each other

timeaxis2 = 1 : 1 + (secperstripe*rate); # should be a bit more than N * 800 msec
timeaxis2 = timeaxis2 * (1000/rate); # convert to msec units
waveexcerpts = zeros(1, length(timeaxis2));
minexcerpts = zeros(1, length(timeaxis2)) + 0.5;
maxexcerpts = zeros(1, length(timeaxis2)) - 0.5;
excerptsize = (int32)(0.8 * rate); # at most 800 msec per excerpt
excerptcolumns = floor(secperstripe/0.8);
nexcerpts = excerptcolumns * stripes;

for excerptcolumn = 1:excerptcolumns
excerptoffset = (int32)(excerptsize * (excerptcolumn-1)) + 1;
avgblock = (excerptcolumns * (stripe-1)) + excerptcolumn - 1;
blockstart = round((nbeats-1) * avgblock/nexcerpts) + 1;
blockend = round((nbeats-1) * (avgblock+1)/nexcerpts);
blockrangestart(excerptcolumn) = beats(blockstart) - theshift;
blockrangeend(excerptcolumn) = beats(blockend) - theshift;

for beat = blockstart:blockend
  prevt = (int32)((beats(beat) - theshift) * rate);
  nextt = (int32)((beats(beat+1) - theshift) * rate);
  if (nextt > prevt + excerptsize)
    nextt = prevt + excerptsize; # only use first 800 msec
  endif
  realexcerptsize = nextt - prevt; # at most excerptsize
  wavesnip = wave(prevt:nextt);
  minexcerpts(excerptoffset:excerptoffset+realexcerptsize) = min(minexcerpts(excerptoffset:excerptoffset+realexcerptsize), wavesnip);
  maxexcerpts(excerptoffset:excerptoffset+realexcerptsize) = max(maxexcerpts(excerptoffset:excerptoffset+realexcerptsize), wavesnip);
  # wavesnip = wavesnip - mean(wavesnip); # would amplify drift edge artifacts
  # note: block sizes can vary per column due to rounding!
  wavesnip = wavesnip * (1 / (1+blockend-blockstart));
  waveexcerpts(excerptoffset:excerptoffset+realexcerptsize) = waveexcerpts(excerptoffset:excerptoffset+realexcerptsize) .+ wavesnip;
endfor # average over all beats in one block

endfor # loop over columns in one stripe

plot(timeaxis2, minexcerpts, 'r', maxexcerpts, 'r', timeaxis2, waveexcerpts, "k", "linewidth", 2);
stripemin = min(min(waveexcerpts), -0.90);
stripemax = max(max(waveexcerpts), 1.1);
for i = 0:excerptcolumns-1
  line([i*0.8*rate i*0.8*rate], [stripemin stripemax]); # draw vertical line: excerpt boundary
  line([(i*0.8+theshift)*rate (i*0.8+theshift)*rate], [stripemin stripemax], 'color', 'r'); # vertical line: QRS
  text((i*0.8+0.25)*rate, -0.7, sprintf('%i-%i', round(blockrangestart(i+1)), round(blockrangeend(i+1))));
endfor
xlim([0 max(timeaxis2)]);
ylim([stripemin stripemax]); # may clip the min/max plot lines, intentionally
# if (stripemax < 1.1 && stripemin > -1.1)
#   ylim([-1.1 1.1]); # try to stick to a normalized default vertical scale
# else
#  axis('tight'); # allow bigger Y amplitude but stay reasonable
# endif
if (stripe < 2)
  title(wavename);
else
  title('Gemittelte EKG-Pulsformen pro Abschnitt');
endif
xlabel('Zeit virtuell [ms]');
ylabel('Spannung [unkalibriert]');
grid on;
grid minor on;

endfor # loop over stripes

print('gemitteltes_ekg.pdf');
endif



# profile off; # <<<
# profshow(profile('info')); # <<<

# Profile: handles 10 min of ECG recording in 30/16 sec, with/without dcdecay
# (most REAL time is spent in the for loop of the dcdecy anti-high pass step)
# If dcdecay is not used, use cumsum: Saves huge FOR loop with +, +, * :-)

# Time spent to get the whole recorded data down to 1000 samples/second:
# decimate filter 0.6 sec (2 calls) versus resample upfirdn 5 sec versus
# interp fft+ifft+fftfilt = 3.2 sec (2+2+1 calls)

# Function  Time  Calls (2016-07-28 version)
# drawnow   3.190   9 (plot)
# fft       1.569   2 (interp)
# ifft      1.384   2 (interp)
# filter    1.209  16 (notches, sgolayfilt, low pass, moving avg, decimate)
# binary +  0.777 [2 per sample after downsample, anti-highpass with dcdecay]
# binary *  0.467 [1 per sample after downsample, anti-highpass with dcdecay]
# binary >  0.424 [1 per sample after downsample, QRS finder]
# fftfilt   0.422   1 (interp)
# binary .* 0.351 186 (various filters and direct invocations)
# set       0.291 262