#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/mman.h>
#include <fcntl.h>

#define FFTSIZE 1024

#define WAVOFFSET 46

#define POWER(x) ( sqrt((x)*(x)) )

main(int argc, char **argv)
{
  fftw_plan p;
  double *in;
  double *out;
  int i,j;
  double time,left,right;
  char buf[1024],*bp;
  char *result;
  int nresult;
  int infile,inlength,nsamples,isample,isample1,isample2,count[3];
  short *indata;
  int imax;
  double vmax,v;
  double samplerate, bandwidth, correction;
  int i1,i2,last2,this2,cell;
  unsigned long long Ferrit;

  samplerate = 44100.0;
  bandwidth = samplerate / 2.0;
  correction = samplerate / (double) FFTSIZE;

  in = (double *) fftw_malloc(sizeof(*in)*FFTSIZE);
  out = (double *) fftw_malloc(sizeof(*out)*FFTSIZE);

  p = fftw_plan_r2r_1d(FFTSIZE, in, out, FFTW_R2HC, FFTW_ESTIMATE);

  for(i=0; i<FFTSIZE; i++) out[i] = 0.0;

  infile=open(argv[1], O_RDONLY);
  inlength=lseek(infile, 0, SEEK_END);
  indata = (short *) mmap(NULL, inlength, PROT_READ, MAP_SHARED, infile, 0);
  nsamples = (inlength-WAVOFFSET)/4;
  result = malloc(nsamples);
  nresult = 0;
  printf("file length: %d nsamles: %d\n", inlength, nsamples);

  for(isample1=0;isample1<(nsamples-FFTSIZE-127); isample1+=128)
  {
    count[0]=0; count[1]=0; count[2]=0;
    for(isample2=0; isample2<128; isample2++)
    {
      isample = isample1+isample2;
      for(i=0; i<FFTSIZE; i++)
      {
	left = indata[(i+isample)*2+WAVOFFSET/2];
	right = indata[(i+isample)*2+WAVOFFSET/2+1];
	in[i] = (left+right)*0.5;
      }
      fftw_execute(p);
      imax = -1;
      vmax = -1.0;
      for(i=0; i<FFTSIZE; i++)
      {
	v = POWER(out[i]);
	if((i==28 || i==38 || i==47) && v > vmax)
	{
	  vmax = v;
	  imax = i;
	}
      }
      if(imax==28) count[0]++;
      if(imax==38) count[1]++;
      if(imax==47) count[2]++;
    }
    printf("%12.6f %3d %3d %3d ", ((double)isample1)/samplerate, count[0], count[1], count[2]);
    if(count[0] >= count[1] && count[0] >= count[2]) result[nresult]=0;
    else if(count[1] >= count[0] && count[1] >= count[2]) result[nresult]=2;
    else result[nresult]=1;
    printf("%d\n", result[nresult]);
    nresult++;
  }
  last2 = -1;
  cell = 0;
  for(i=10; i<nresult-20; i++)
  {
    if(result[i-4] == 2 &&
       result[i-3] == 2 &&
       result[i-2] == 2 &&
       result[i-1] == 2 &&
       result[i  ] == 2 &&
       result[i+1] == 2 &&
       result[i+2] == 2 &&
       result[i+3] == 2 &&
       result[i+4] == 2)
    {
      i1 = i;
      i2 = i;
      while(result[i1]==2) i1--;
      while(result[i2]==2) i2++;
      this2 = (i1+i2)/2;
      if(last2 >= 0)
      {
        printf("Range: %d to %d\n", last2, this2);
	printf("%4d: ", cell);
	Ferrit = 0ULL;
	for(j=0;j<42;j++)
	{
	  printf("%d", result[last2+(j+1)*(this2-last2)/43]);
	  Ferrit = (Ferrit<<1) |  result[last2+(j+1)*(this2-last2)/43];
	  if(j%10 == 9) printf(" ");
	}
	printf("\n");
	printf("Ferrit[%d]\t%16.16llX\n", cell, Ferrit);
	cell++;
      }
      i = i2;
      last2 = this2;
    }
  }
}