# finestra 128 e a=12, con round(p+0.09)
# soglia 0.25

function [rcor l] = arduino (x, y)
  if (length(x) > length(y))
    signal = x; proto = y;
  else
    signal = y; proto = x;
  endif
  pvar = var(proto,1);
  plen = length(proto);		# length of prototype is equal to
       	 			# length of running correlation window
  slen = length(signal);	# length of signal
  ps = prepad(signal,slen+1);	# prepend signal with a single 0
  scsum = cumsum(ps);		# signal cumulative sum
  scsqr = cumsum(ps.^2);	# signal cumulative sum of squares
  we = 1+(plen:slen);		# moving window end
  ws = we-plen;			# moving window start
  srsum = scsum(we)-scsum(ws);	# signal running sum
  srsqr = scsqr(we)-scsqr(ws);	# signal running sum of squares
  cov = fftconv(signal,rot90(proto,2));
  cw = plen:slen;		# central window of convolution discarding tails
  rcor = max(0,cov(cw)).^2;	# discard negative covariances
  rvar = (plen*srsqr-srsum.^2);	# running variance
  lowdev = 5;			# minimum deviation of a blink
  bw = 8;			# blink width
  lowv = lowdev^2*bw/plen;	# low variance
  lowvar = lowv*plen^2;		# low variance threshold
  rcor(rvar<lowvar) = 0;	# discard low variance samples
  rcor .*= plen / sumsq(proto) ./ rvar;
  l = [signal(plen/2+(0:slen-plen)) srsum srsqr lowvar*ones(size(rvar)) rvar cov(cw) rcor];
  l = [l(65:end,1) l(1:end-64,2:5) l(4:end-61,6:7)]; # align with logging from arduino
endfunction

## Local Variables:
## fill-column: 100
## End:
