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;
   for (i=0;i<n;i++) {
      mx += x[i];
      my += y[i];
   }
   mx /= n;//lasketaan keskiarvot
   my /= n;

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

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

         // Or should it be (?)
  /*       if (j < 0 || j >= n)
            sxy += (x[i] - mx) * (-my);
         else
            sxy += (x[i] - mx) * (y[j] - my);
*/
      }
      r = sxy / denom;
      if(r>prev_r)
      {
          prev_r=r;
          *calcdelay=delay;
      }
      /* r is the correlation coefficient at "delay" */
   }
   return prev_r;
}