fourier1p.c

double fourier1p(double x_in[], double n, int length, double* phase_r, double* phase_i) { 
  double x_complex[2] = { 0, 0 }; 
  int i; 
  
  for(i = 0; i < length; i++) { 
    x_complex[0] += x_in[i] * cos(M_PI * 2 * i * n / (double) length); 
    x_complex[1] += x_in[i] * sin(M_PI * 2 * i * n / (double) length); 
  } 
  
  double norm = sqrt(x_complex[0]*x_complex[0] + x_complex[1]*x_complex[1]); 
  *phase_i = x_complex[1] / norm; 
  *phase_r = x_complex[0] / norm; 
  return norm / length; 
}