/**************************************************************************

              	    C O R R E L A T E . C   

***************************************************************************/
/*****************************************************************************
FILENAME   : correlate.c
PURPOSE    : This performs the normalized correlation of Isthmii.
*****************************************************************************/
#include <math.h>
#define M 5

float correlate();
float sum_sqr_2();
float sum_both_2();
float correl_control();


float correl_control(f,g, gs, ge, fs, fe, correl_offset)
 float f[], g[]; 
 int gs, ge, fs, fe;
 int *correl_offset;
{
  int n;
  int sum, m;
  float num; 
  float g_bar, g_sqrs;
  float biggest = -3.0;
  if ( (ge-gs) < (fe-fs) )
     n = ge-gs;
  else
     n = fe-fs;
  
  sum   = sum_funct (g, gs, n); 
  g_bar = (float)sum / n;
  g_sqrs = sum_sqr_2 (g, gs, n, g_bar);
  for (m = -M; m <= M; m++)
  {
    num = correlate (g_bar, g_sqrs, f,g, gs, fs+m, n);
    if (num > biggest)
    {
      biggest = num;
      *correl_offset = m;
    }
  }
  return (biggest);
}

sum_funct (array, s,n)
  float array[]; int s,n;
{
  int x;
  int sum = 0;
  for (x=s;x<= s+n;x++)
  {
    sum = sum + (int)array[x];
  }
  return (sum);
}

float sum_both_2 (f,g, gs, n, fs, f_bar, g_bar)
  float f[], g[]; int gs, n, fs;
  float f_bar, g_bar;
{
  int x;
  float numer = 0.0;
  int ind;
  ind = fs; 
  for (x=gs;x<=gs+n;x++)
  {
    numer = numer + ( (f[ind] - f_bar) * (g[x] - g_bar) );
    ind++; 
  }
  return (numer);
}

float sum_sqr_2 (array, s,n, ave)
  float array[]; int s,n;
  float ave;
{
  int x;
  float num;
  float sum = 0.0;
  for (x=s;x<=s+n;x++)
  {
    num = array[x] - ave; 
    sum = sum + (num * num);
  }
  return (sum);
}

float correlate(g_bar, g_sqrs, f,g, gs, fs, n)
  float g_bar, g_sqrs, f[],g[];
  int gs, fs, n;
{
  int one, two, mt; 
  float numer, denom, first, second;
  float correl_value;
  float f_bar, f_sqrs;

  f_bar = sum_funct (f, fs, n); 
  f_bar = (float)f_bar / n; 
  numer = sum_both_2 (f,g, gs, n, fs, f_bar, g_bar);
  
  f_sqrs = sum_sqr_2 (f, fs, n, f_bar);
  denom = sqrt (f_sqrs * g_sqrs); 
  correl_value = (float)numer/denom; 
  return (correl_value);
}

#ifdef ROGER
  printf ("g_bar: %f", g_bar);
  printf ("  g_sqrs: %f\n", g_sqrs);
  printf (" f_bar: %f", f_bar);
  printf ("  f_sqrs: %f\n", f_sqrs);
  printf ("numer: %f  denom: %f corr: %f\n", 
           numer, denom, correl_value);
    printf ("=========================================\n"); 
    printf ("m: %d ----> correl_value: %f\n", m, num);
    printf ("=========================================\n"); 
#endif