Re: running Person's correlation

[Francesco Potortì]

I had written:

>>I have a signal (long) and a template (short, fixed).  I have to compute
>>the Pearson's correlation of the short signal with a sliding window of
>>the long signal.  This is a convolution where each sample is divided by
>>the (fixed) standard deviation of the short signal and the running
>>standard deviation of the long signal.
>>
>>The only loopless way I can think of is to compute a running sum, a
>>running sum of squares, and use them to compute a running standard
>>deviation to be multiplied with the convolution.  Any more
>>straightforward methods?

And then:

>Thanks again to all who have tried an answer.  Here is my solution to
>the problem.  If it turns out that there is already a canned solution
>analogous to mine, I'll adopt that one, else I'm open to suggestions and
>possibly willing to create a better solution to incorporate into the
>signal package (this one is subject to numerical cancellation).

However, I see that the code was mangled, and that a very simple change
improves code readability and reduces cancellation problems, so here is
a second version.

Also, the example provided does not work with most recent Octave.  To
test the function, try this:

octave> t=interp1([1 7 15],[1 0 0],1:15);
octave> runningcor(t,t)
ans =  1.00000
octave> plot(runningcor(t,[sin(1:15).^2 t sqrt(1:7)]))


===File ~/math/octavelib/signal/runningcor.m================
## © 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:
============================================================

-- 
Francesco Potortì (ricercatore)        Voice:  +39.050.621.3058
ISTI - Area della ricerca CNR          Mobile: +39.348.8283.107
via G. Moruzzi 1, I-56124 Pisa         Skype:  wnlabisti
(entrance 20, 1st floor, room C71)     Web:    http://fly.isti.cnr.it