Page History

double Statistic::Correlation(double*x,double*y,int maxdelay,int *calcdelay,int n)

{
//http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/correlate/
int i,j,delay;

     double mx,my,sx,sy,sxy,denom,r,prev_r;

     /\* Calculate the mean of the two series x\[\], y\[\] \*/
     mx = 0;
     my = 0;
&nbsp;&nbsp;   for (i=0;i<n;i++) {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      mx \+= x\[i\];
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      my \+= y\[i\];
&nbsp;&nbsp;   }
&nbsp;&nbsp;   mx /= n;//lasketaan keskiarvot
&nbsp;&nbsp;   my /= n;

&nbsp;&nbsp;   /\* Calculate the denominator \*/
&nbsp;&nbsp;   sx = 0;
&nbsp;&nbsp;   sy = 0;
&nbsp;&nbsp;   for (i=0;i<n;i++)
&nbsp;&nbsp;   {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      sx \+= (x\[i\] - mx) * (x\[i\] - mx);
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      sy \+= (y\[i\] - my) * (y\[i\] - my);
&nbsp;  }
&nbsp;&nbsp;   denom = sqrt(sx*sy);
&nbsp;&nbsp;   r=0;
&nbsp;&nbsp; \   *calcdelay=0;
&nbsp;&nbsp;   prev_r=0;
&nbsp;&nbsp;   /\* Calculate the correlation series \*/
&nbsp;&nbsp;   for (delay=-maxdelay;delay<maxdelay;delay++)
&nbsp;&nbsp;&nbsp;    {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      sxy = 0;

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      for (i=0;i<n;i++) {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;         j = i + delay;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;         if (j < 0 \|\| j >= n)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;            continue;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;         else
            sxy \+= (x\[i\] - mx) * (y\[j\] - my);

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;         // Or should it be \(?)
&nbsp;  /*&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;       if (j < 0 \|\| j >= n)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;            sxy \+= (x\[i\] - mx) * (-my);
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;         else
            sxy \+= (x\[i\] - mx) * (y\[j\] - my);
\*/
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      r = sxy / denom;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      if(r>prev_r)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      {
          prev_r=r;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; \          *calcdelay=delay;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      }



&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;      /\* r is the correlation coefficient at "delay" \*/

&nbsp;&nbsp;   }
&nbsp;&nbsp;   return prev_r;
}