#include <stdio.h>
#include <math.h>

#define N 8000000
double data[N];
double x[N];
double y[N];

void hilbert(double *x, double *y, double *src, int n) {
   for(int i = 0; i < n; i++) {
      double t = 0;
      for(int j = -127; j < 127; j+=2) {
         if(i+j >= 0 && i+j < N)
            t += src[i+j] * (2.0/M_PI) /j;
      }

      x[i] = src[i];
      y[i] = t;
   }
}

void dump(double *x, double *y, int start, int end) {
   for(int i = start; i < end; i++) {
      double mag = sqrt(x[i]*x[i] + y[i]*y[i]);
      double phase_last = atan2(y[i-1],x[i-1]);
      double phase      = atan2(y[i],x[i]);
      double angular_v  = (phase_last-phase);
      if(angular_v < -M_PI) 
         angular_v += 2*M_PI;
      else if(angular_v > M_PI)
         angular_v -= 2*M_PI;
      printf("%i, %7.4f, %7.4f, %7.4f, %7.4f, %7.4f\n", i, x[i], y[i], phase, mag, angular_v/(2*M_PI));
   }
}

int main(void) {
   double phase = 0.0;
   for(int i = 0; i < N; i++) {
      data[i] = sin(phase) * (1.0 - 0.5*i/N);
      phase += 2*M_PI*(0.05+ 0.40*i/N);
   }

   hilbert(x,y,data,N);

   printf("i, x, y, phase, mag, angular_v\n");   
   dump(x, y, 1000,  1200);
   dump(x, y, N/2,   N/2+200);
   dump(x, y, N-400, N-200);
}