## © Francesco Potortì 2013
## $Revision: 1.2 $
## This program is licensed under the terms of the GPL v3

## Running Pearson's correlation between X and Y
##
## Useful as a detector with automatic gain control.
## Compares the two arguments vectors X and Y, giving a result of
## lenght abs(length(x)-lenght(y))+1.
## It takes the shortest of X and Y (the prototype) and computes a
## correlation of it versus a moving window of the longest (the signal).
## The absolute magnitude of the result is not greater than 1.

function rcor = runningcor (x, y)
  if (length(x) > length(y))
    signal = x; proto = y;
  else
    signal = y; proto = x;
  endif
  proto -= mean(proto);		# standardise prototype for simpler
  proto /= std(proto,1);	# code and better numerical stability
  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
  srmean = srsum/plen;		# signal running mean
  srmsqr = srsqr/plen;		# signal running mean of squares
  srvar = (srmsqr-srmean.^2);	# signal running population variance
  srstd = sqrt(srvar);		# signal running population standard deviation
  p = (rot90(proto,2));		# invert prototype
  rprod = fftconv(signal,p);	# running product of signal and prototype
  cw = plen:slen;		# central window of convolution discarding tails
  rcov = rprod(cw)/plen;	# running covariance of signal and prototype
  rcor = rcov ./ srstd;		# running Pearson's correlation
endfunction

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