Re: running Person's correlation

[Francesco Potortì]

>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?

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).

As a usage example, try:

octave> t=triangle_lw(15,0.1);
octave> plot(runningcor(t',[sin(1:15).^2 t' sqrt(1:7)]))

You'll see that the graph has a maximum at 16, that is, where the second
argument (the signal) equals the first (the prototype).


# © Francesco Potortì 2013
# 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
  plen = length(proto);         # length of prototype, same as
                                #  length of running correlation window
  pmean = mean(proto);             # mean of prototype
  pstd = std(proto,1);               # population standard deviation of
  prototype
  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-pmean*srmean; # running covariance of signal and
  prototype
  rcor = rcov / pstd ./ 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