help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: running Person's correlation


From: Alois Schloegl
Subject: Re: running Person's correlation
Date: Mon, 02 Sep 2013 09:45:07 +0200
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:17.0) Gecko/20130827 Icedove/17.0.8

On 09/02/2013 08:50 AM, Francesco Potortì wrote:
> Hi all,
> 
> 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?
> 


filter() is the function you need to vectorize your code. However, the
big question is whether you need a causal algorithms (i.e. only samples
from the past can be used), or whether you can use a symetric sliding
window. It will make a big difference how to estimate the running mean
and standard deviation. Below, there is an outline of the code you are
looking for, the code snippets are not tested

% acausal method using symetric triangular window
N=length(template);
B=ones(1,N); A=N; %% rectangular window, because of
% mean_signal=filter(B,A,signal);   %% filter would causes a time delay
of N/2
mean_signal=filtfilt(B,A,signal);   %% thus we use filtfilt, forward
backward filtering will result in a triangular and symetric filter window
sd_signal = sqrt(filtfilt(B,A,(signal-mean_signal).^2));        %% standard
deviation
r =
filtfilt(template(end:-1:1),std(template),(signal-mean(signal))./sd_signal);

For a causal solution, I'd recommend using an exponential window to
estimate mean and sd.
UC = 0.001;     %% update coefficient, should be in the order of 1/N
B = UC; A = [1; UC-1];
mean_signal=filter(B,A,signal);   %% causes an average time delay of
sd_signal = sqrt(filter(B,A,(signal-mean_signal).^2));
r =
filter(template(end:-1:1),std(template),(signal-mean(signal))./sd_signal);

Some of these ideas are summarized in [1].


Best regards,
  Alois


[1] A. Schlögl, C. Vidaurre, K.-R. Müller
Adaptive Methods in BCI Research - An Introductory Tutorial.
(Eds.) B. Graimann, B. Allison, G. Pfurtscheller.
in "Brain Computer Interfaces - Revolutionizing Human-Computer
Interfaces" Springer, 2010, p. 331-355.
Preprint available online
http://pub.ist.ac.at/~schloegl/publications/Schloegl2010Adaptive.pdf







reply via email to

[Prev in Thread] Current Thread [Next in Thread]