#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdint.h>
#include <time.h>
#include <fftw3.h>
#include <alsa/asoundlib.h>

#include "fano.h"
#include "wsprd_utils.h"
//#include "wsprsim_utils.h"
#define GNUPLOT "gnuplot -persist"
//#define GNUPLOT "gnuplot "

#define max(x,y) ((x) > (y) ? (x) : (y))
// Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT,
// FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
#define PATIENCE FFTW_ESTIMATE
fftwf_plan PLAN1,PLAN2,PLAN3;

unsigned char pr3[162]=
{1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0,
    0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1,
    0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1,
    1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1,
    0,0,1,0,1,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0,
    0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1,
    0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1,
    0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,0,1,1,0,
    0,0};

unsigned long nr;
int printdata=0;
int iterasi = 0;

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

unsigned long readwavfile(char *infile, float *idat, float *qdat )

{
  FILE *fp;
  char header[44];
  
  size_t i, j, npoints;
  int nfft1, nfft2, nh2, i0;
  double df;

  int16_t *buf2_int16, *buf2;   // buffer integer 16bits
  float *buf2_float;

  int rtime =114 ; // 114 detik estimasi jitter transmisi WSPR
  int buffer_frames_48 = rtime*48000; // alokasi buffer samping 48KHz sesuai jumlah sampel
  int buffer_frames_12 = rtime*12000; // alokasi buffer sampleing 12 KHz

  //  infile=malloc(40*sizeof(char));
  buf2 = malloc(buffer_frames_48 * 4);
  buf2_int16 = malloc(buffer_frames_48 * 2);
  buf2_float = malloc(buffer_frames_48*sizeof(float));
   

  //  printf("\n%s\n",infile);
  if(strstr(infile,".raw")==NULL && strstr(infile,".WAV")==NULL ){
    printf("\nAnda tidak memasukkan argumen dengan benar..!!!!\n\n");
    exit(1);
  } 
  if(strstr(infile,".raw") ){
    fp = fopen(infile,"rb");
    if (fp == NULL) {
      fprintf(stderr, "\nTidak ada file %s... \n",infile);
      exit(1);
    }

    nr=fread(buf2,4,buffer_frames_48,fp);       //Read raw data
    fclose(fp);


  /* ambil 1 sampel mono  */

  for (i = 0; i < buffer_frames_48; i++){
    buf2[i]=buf2[2*i];
 }
  /**********************/

  /* anti aliasing filter  */

  float pi = 3.141592654; 
  float fc = 2500;
  float frq=48000;
  float a1,b1,c1,tau1;

  tau1=1/(2*pi*fc);
  a1=(1.0-1.0/(2.0*tau1*frq))/(1.0+1.0/(2.0*tau1*frq));
  b1=c1=(1.0/(2.0*tau1*frq))/(1.0+1.0/(2.0*tau1*frq));

  buf2_int16[0]=buf2[0];
  buf2_float[0]=(float)buf2[0];
  for (i = 1; i < buffer_frames_48; i++){
    buf2_float[i]=a1*buf2_float[i-1]+b1*(float)buf2[i]+c1*(float)buf2[i-1];
    buf2_int16[i]=(int16_t)buf2_float[i];
  }
  /*************************/

  /**** decimasi menjadi 12 KHz sampling ***/

  for (i = 0; i < buffer_frames_12; i++){
    buf2_int16[i] = buf2_int16[i*4];
  }
  /*
  for (i = buffer_frames_12; i < buffer_frames_48; i++){
    buf2_int16[i] = 0;
  }
  */
  
  /*****************************************/
  }

  else  if(strstr(infile,".WAV") != 0){
    fp = fopen(infile,"rb");
    if (fp == NULL) {
      fprintf(stderr, "\nTidak ada  file %s... \n",infile);
      exit(1);
    }

    nr=fread(header,2,22,fp);            //Read and ignore header
    nr=fread(buf2_int16,2,buffer_frames_12,fp);       //Read raw data
    fclose(fp);
  }
                                                                   
  /*
  for (i = 0; i < buffer_frames_12; i++){
    if(i>100000){
      printf("buf2_int16[%d]= %d\n",i,buf2_int16[i]);
      getchar();
    }
  }
  */


  nfft2=46080; /* ini adalah jumlah sampel setelah proses downsampling */
  nh2=nfft2/2;
    

  nfft1=nfft2*32;      //need to downsample by a factor of 32
  df=12000.0/nfft1;
  i0=1500.0/df+0.5;
  npoints=114*12000;

    
  float *realin;
  fftwf_complex *fftin, *fftout;
    
    realin=(float*) fftwf_malloc(sizeof(float)*nfft1);
    fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*(nfft1/2+1));
    PLAN1 = fftwf_plan_dft_r2c_1d(nfft1, realin, fftout, PATIENCE);
    
    for (i=0; i<npoints; i++) {
      realin[i]=buf2_int16[i]/32768.0;
    }
    

    for (i=npoints; i<(size_t)nfft1; i++) {
        realin[i]=0.0;
    }
    
    free(buf2);
    free(buf2_int16);
    free(buf2_float);
    fftwf_execute(PLAN1);
    fftwf_free(realin);

    fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*nfft2);

    for (i=0; i<(size_t)nfft2; i++) {
      j=i0+i;                            /* i0===1500Hz; -nh2----->i0----->nh2,   nfft2=2nh2  */
      if( i>(size_t)nh2 ) j=j-nfft2;
      fftin[i][0]=fftout[j][0];
      fftin[i][1]=fftout[j][1];
    }
    fftwf_free(fftout);

    fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*nfft2);
    PLAN2 = fftwf_plan_dft_1d(nfft2, fftin, fftout, FFTW_BACKWARD, PATIENCE);
    fftwf_execute(PLAN2);
    
    for (i=0; i<(size_t)nfft2; i++) {
        idat[i]=fftout[i][0]/1000.0;
        qdat[i]=fftout[i][1]/1000.0;
    }

    fftwf_free(fftin);
    fftwf_free(fftout);
    return nfft2;
}

//***************************************************************************
void sync_and_demodulate(float *id, float *qd, long np,
                         unsigned char *symbols, float *f1, int ifmin, int ifmax, float fstep,
                         int *shift1, int lagmin, int lagmax, int lagstep,
                         float *drift1, int symfac, float *sync, int mode)
{
  /*********************************************************************************
     * mode = 0: tanpa pencarian frekuensi atau drift. Dapatkan time lag terbaik    *
     *        1: tanpa pencarian time lag dan  drift. Dapatkan frekuensi terbaik    *
     *******************************************************************************/
    
    static float fplast=-10000.0;
    static float dt=1.0/375.0, df=375.0/256.0;
    static float pi=3.14159265358979323846;
    float twopidt, df15=df*1.5, df05=df*0.5;

    int i, j, k, lag;
    float i0[162],q0[162],i1[162],q1[162],i2[162],q2[162],i3[162],q3[162];
    float p0,p1,p2,p3,cmet,totp,syncmax;
    float c0[256],s0[256],c1[256],s1[256],c2[256],s2[256],c3[256],s3[256];
    float dphi0, cdphi0, sdphi0, dphi1, cdphi1, sdphi1, dphi2, cdphi2, sdphi2, dphi3, cdphi3, sdphi3;
    float f0=0.0, fp, ss, fbest=0.0;
    int best_shift = 0, ifreq;

    syncmax=-1e30;
    if( mode == 0 ) {ifmin=0; ifmax=0; fstep=0.0; f0=*f1;}
    if( mode == 1 ) {lagmin=*shift1;lagmax=*shift1;f0=*f1;}
    
    twopidt=2*pi*dt;
    for(ifreq=ifmin; ifreq<=ifmax; ifreq++) {
        f0=*f1+ifreq*fstep;
        for(lag=lagmin; lag<=lagmax; lag=lag+lagstep) {
            ss=0.0;
            totp=0.0;
            for (i=0; i<162; i++) {
	      fp = f0 + (*drift1/2.0)*((float)i-81.0)/81.0;  //ini rumus frekuensi drift linier
                if( i==0 || (fp != fplast) ) {  // generator sin/cos cepat
                    dphi0=twopidt*(fp-df15);
                    cdphi0=cos(dphi0);
                    sdphi0=sin(dphi0);
                    
                    dphi1=twopidt*(fp-df05);
                    cdphi1=cos(dphi1);
                    sdphi1=sin(dphi1);
                    
                    dphi2=twopidt*(fp+df05);
                    cdphi2=cos(dphi2);
                    sdphi2=sin(dphi2);
                    
                    dphi3=twopidt*(fp+df15);
                    cdphi3=cos(dphi3);
                    sdphi3=sin(dphi3);
                    
                    c0[0]=1; s0[0]=0;
                    c1[0]=1; s1[0]=0;
                    c2[0]=1; s2[0]=0;
                    c3[0]=1; s3[0]=0;


                    for (j=1; j<256; j++) {
                        c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0;
                        s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0;
                        c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1;
                        s1[j]=c1[j-1]*sdphi1 + s1[j-1]*cdphi1;
                        c2[j]=c2[j-1]*cdphi2 - s2[j-1]*sdphi2;
                        s2[j]=c2[j-1]*sdphi2 + s2[j-1]*cdphi2;
                        c3[j]=c3[j-1]*cdphi3 - s3[j-1]*sdphi3;
                        s3[j]=c3[j-1]*sdphi3 + s3[j-1]*cdphi3;
                    }

		      /* generator sin/cos standar */

		    /*
                    for (j=1; j<256; j++) {
			c0[j]=cos(dphi0*(double)j);
			s0[j]=sin(dphi0*(double)j);
			c1[j]=cos(dphi1*(double)j);
			s1[j]=sin(dphi1*(double)j);
			c2[j]=cos(dphi2*(double)j);
			s2[j]=sin(dphi2*(double)j);
			c3[j]=cos(dphi3*(double)j);
			s3[j]=sin(dphi3*(double)j);
                    }
		    */
                    fplast = fp;
                }
                
                i0[i]=0.0; q0[i]=0.0;
                i1[i]=0.0; q1[i]=0.0;
                i2[i]=0.0; q2[i]=0.0;
                i3[i]=0.0; q3[i]=0.0;
                
                for (j=0; j<256; j++) {
                    k=lag+i*256+j;
                    if( (k>0) && (k<np) ) {
                        i0[i]=i0[i] + id[k]*c0[j] + qd[k]*s0[j];
                        q0[i]=q0[i] - id[k]*s0[j] + qd[k]*c0[j];
                        i1[i]=i1[i] + id[k]*c1[j] + qd[k]*s1[j];
                        q1[i]=q1[i] - id[k]*s1[j] + qd[k]*c1[j];
                        i2[i]=i2[i] + id[k]*c2[j] + qd[k]*s2[j];
                        q2[i]=q2[i] - id[k]*s2[j] + qd[k]*c2[j];
                        i3[i]=i3[i] + id[k]*c3[j] + qd[k]*s3[j];
                        q3[i]=q3[i] - id[k]*s3[j] + qd[k]*c3[j];
                    }
                }
                p0=i0[i]*i0[i] + q0[i]*q0[i];
                p1=i1[i]*i1[i] + q1[i]*q1[i];
                p2=i2[i]*i2[i] + q2[i]*q2[i];
                p3=i3[i]*i3[i] + q3[i]*q3[i];

                p0=sqrt(p0);
                p1=sqrt(p1);
                p2=sqrt(p2);
                p3=sqrt(p3);
                
                totp=totp+p0+p1+p2+p3;
                cmet=(p1+p3)-(p0+p2);
                ss = (pr3[i] == 1) ? ss+cmet : ss-cmet;
	    }
            ss=ss/totp;
            if( ss > syncmax ) {          //simpan paramater terbaik
                syncmax=ss;
                best_shift=lag;
                fbest=f0;
            }
        } // lag loop
    } //freq loop
    
    if( mode <=1 ) {                       // kirim kembali parameter terbaik kepada caller 
        *sync=syncmax;
        *shift1=best_shift;
        *f1=fbest;
        return;
    }
    
    return;
}

void noncoherent_sequence_detection(float *id, float *qd, long np,
                         unsigned char *symbols, float *f1,  int *shift1,
                         float *drift1, int symfac, int *nblocksize)
{

    static float fplast=-10000.0;
    static float dt=1.0/375.0, df=375.0/256.0;
    static float pi=3.14159265358979323846;
    float twopidt, df15=df*1.5, df05=df*0.5;
    
    int i, j, k, lag, itone, b, nblock, nseq;
    float xi[512],xq[512];
    float is[4][162],qs[4][162];
    float p[512],fac;
    float c0[257],s0[257],c1[257],s1[257],c2[257],s2[257],c3[257],s3[257];
    float dphi0, cdphi0, sdphi0, dphi1, cdphi1, sdphi1, dphi2, cdphi2, sdphi2,
    dphi3, cdphi3, sdphi3;
    float f0, fp, fsum=0.0, f2sum=0.0, fsymb[162];
    
    twopidt=2*pi*dt;
    f0=*f1;
    lag=*shift1;
    nblock=*nblocksize;
    nseq=1<<nblock;

    for (i=0; i<162; i++) {
        fp = f0 + (*drift1/2.0)*((float)i-81.0)/81.0;
        if( i==0 || (fp != fplast) ) {  // generator sin/cos cepat
            dphi0=twopidt*(fp-df15);
            cdphi0=cos(dphi0);
            sdphi0=sin(dphi0);
                    
            dphi1=twopidt*(fp-df05);
            cdphi1=cos(dphi1);
            sdphi1=sin(dphi1);
                    
            dphi2=twopidt*(fp+df05);
            cdphi2=cos(dphi2);
            sdphi2=sin(dphi2);
                    
            dphi3=twopidt*(fp+df15);
            cdphi3=cos(dphi3);
            sdphi3=sin(dphi3);
                    
            c0[0]=1; s0[0]=0;
            c1[0]=1; s1[0]=0;
            c2[0]=1; s2[0]=0;
            c3[0]=1; s3[0]=0;
                    

            for (j=1; j<257; j++) {
                c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0;
                s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0;
                c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1;
                s1[j]=c1[j-1]*sdphi1 + s1[j-1]*cdphi1;
                c2[j]=c2[j-1]*cdphi2 - s2[j-1]*sdphi2;
                s2[j]=c2[j-1]*sdphi2 + s2[j-1]*cdphi2;
                c3[j]=c3[j-1]*cdphi3 - s3[j-1]*sdphi3;
                s3[j]=c3[j-1]*sdphi3 + s3[j-1]*cdphi3;
            }

	    /*generator sin/cos standar */

	    /*
	      for (j=0; j<257; j++) {
	      c0[j]=cos(dphi0*(double)j);
	      s0[j]=sin(dphi0*(double)j);
	      c1[j]=cos(dphi1*(double)j);
	      s1[j]=sin(dphi1*(double)j);
	      c2[j]=cos(dphi2*(double)j);
	      s2[j]=sin(dphi2*(double)j);
	      c3[j]=cos(dphi3*(double)j);
	      s3[j]=sin(dphi3*(double)j);
	      }
	    */
            fplast = fp;
        }

        is[0][i]=0.0; qs[0][i]=0.0;
        is[1][i]=0.0; qs[1][i]=0.0;
        is[2][i]=0.0; qs[2][i]=0.0;
        is[3][i]=0.0; qs[3][i]=0.0;

	
        for (j=0; j<256; j++) {
            k=lag+i*256+j;
            if( (k>0) && (k<np) ) {
                is[0][i]=is[0][i] + id[k]*c0[j] + qd[k]*s0[j];
                qs[0][i]=qs[0][i] - id[k]*s0[j] + qd[k]*c0[j];
                is[1][i]=is[1][i] + id[k]*c1[j] + qd[k]*s1[j];
                qs[1][i]=qs[1][i] - id[k]*s1[j] + qd[k]*c1[j];
                is[2][i]=is[2][i] + id[k]*c2[j] + qd[k]*s2[j];
                qs[2][i]=qs[2][i] - id[k]*s2[j] + qd[k]*c2[j];
                is[3][i]=is[3][i] + id[k]*c3[j] + qd[k]*s3[j];
                qs[3][i]=qs[3][i] - id[k]*s3[j] + qd[k]*c3[j];
            }
	}
    }

    /* ========== nseq=2, nblock=1 =========  */

    for (i=0; i<162; i=i+nblock) {
      for (j=0;j<nseq;j++) {
	xi[j]=0.0; xq[j]=0.0;

	b=j;
	itone=pr3[i]+2*b;
	xi[j]=xi[j]+is[itone][i] ; //nilai korelasi infase
	xq[j]=xq[j]+qs[itone][i] ; //nilai korelasi quadratur
            
	p[j]=xi[j]*xi[j]+xq[j]*xq[j];
	p[j]=sqrt(p[j]);
      }

      /**************************************************/
      /* hipotesis --> fsymb[i]== besar == 1; kecil== 0 */
      fsymb[i]=p[1]-p[0];   
      /**************************************************/
    }
    for (i=0; i<162; i++) {              //Normalize the soft symbols
      fsum=fsum+fsymb[i]/162.0;
      f2sum=f2sum+fsymb[i]*fsymb[i]/162.0;
    }
    fac=sqrt(f2sum-fsum*fsum);
    for (i=0; i<162; i++) {
      fsymb[i]=symfac*fsymb[i]/fac;
      if( fsymb[i] > 127) fsymb[i]=127.0;
      if( fsymb[i] < -128 ) fsymb[i]=-128.0;
      symbols[i]=fsymb[i] + 128;
      //      printf("[%d]= %u; ",i,symbols[i]);
    }
    //    getchar();
    iterasi++;

    /* plot simbol */ 

    /*

      if(iterasi==1){
	FILE *fgambar, *gx;
	char gambar_fname7[200];
	strcpy(gambar_fname7,"./gambar7.tmp");
	if((fgambar = fopen(gambar_fname7,"w")) == NULL) {
	  fprintf(stderr, "Cannot open data file '%s'....\n", gambar_fname7);
	  return ;
	}
	for (i=0; i<162; i++ ) {
	  fprintf(fgambar,"%d %d\n",i,symbols[i]);
	}
	fclose(fgambar);
	if((gx = popen(GNUPLOT,"w"))==NULL){ 
	  printf("Error opening pipe to GNU plot. Check if you have it! \n");
	  exit(0);
	}
	fprintf(gx, "plot './gambar7.tmp' with impulses lt 16 title 'tipikal 162 simbol terdekode' \n");
    fflush(gx);
    getchar();

      }

*/
      /* end plot simbol */

    return;
}

//***************************************************************************
int main(int argc, char *argv[])
{
    char cr[] = "(C) 2018, Steven Franke - K9AN";
    (void)cr;
    extern char *optarg;
    extern int optind;
    int i,j,k;
    unsigned char *symbols, *decdata, *channel_symbols;
    signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0};
    char *callsign, *call_loc_pow;
    //    char *ptr_to_infile,*ptr_to_infile_suffix;

    char infile[50];
    char *data_dir=NULL;
    char spots_fname[200];
    char timer_fname[200];

    char uttime[7],date[11];
    int delta,maxpts=65536;
    int writenoise=0, nblocksize;
    int maxdrift;
    int shift1, lagmin, lagmax, lagstep, ifmin, ifmax, worth_a_try, not_decoded;
    unsigned int nbits=81;
    unsigned int npoints, metric, cycles, maxnp;
    float df=375.0/256.0/2;
    float freq0[200],snr0[200],drift0[200],sync0[200];
    int shift0[200];
    float dt=1.0/375.0, dt_print;
    double dialfreq_cmdline=0.0, dialfreq, freq_print;
    double dialfreq_error=0.0;
    float fmin=-110, fmax=110;
    float f1, fstep, sync1, drift1;
    float psavg[512];
    float *idat, *qdat;
    clock_t t0,t00;
    float tfano=0.0,treadwav=0.0,tcandidates=0.0,tsync0=0.0;
    float tsync1=0.0,tsync2=0.0,ttotal=0.0;
    
    time_t t = time(NULL);
    struct tm tm = *localtime(&t);
    struct result { char date[11]; char time[6]; float sync; float snr;
                    float dt; double freq; char message[23]; float drift;
                    unsigned int cycles; int jitter; int blocksize; unsigned int metric; };
    struct result decodes[50];
    
    symbols=calloc(nbits*2,sizeof(char));
    decdata=calloc(11,sizeof(char));
    channel_symbols=calloc(nbits*2,sizeof(char));
    callsign=calloc(13,sizeof(char));
    call_loc_pow=calloc(23,sizeof(char));
    float allfreqs[100];
    char allcalls[100][13];
    for (i=0; i<100; i++) allfreqs[i]=0.0;
    memset(allcalls,0,sizeof(char)*100*13);
    
    int uniques=0, noprint=0;
    unsigned int upaya,lupaya[100];    

    // Parameters used for performance-tuning:
    unsigned int maxcycles=10000;            //Decoder timeout limit
    float minsync1=0.10;                     //First sync limit
    float minsync2=0.12;                     //Second sync limit
    int iifac=8;                             //Step size in final DT peakup
    int symfac=50;                           //Soft-symbol normalizing factor
    //    int block_demod=1;                       //Default is to use block demod on pass 2
    //    int subtraction=1;
    //    int npasses=2;
       //       int npasses=1;

    float minrms=52.0 * (symfac/64.0);      //Final test for plausible decoding
    delta=60;                                //Fano threshold step
    float bias=0.45;                        //Fano metric bias (used for both Fano and stack algorithms)
    
    t00=clock();
    fftwf_complex *fftin, *fftout;
#include "./metric_tables.c"
    
    int mettab[2][256];
    
    if(argc < 2){
      printf("ketik argumen sebagai nama file dengan ext .raw atau .WAV .....\n\n");
      exit(1);
    }

    strcpy(infile, argv[1]);
    idat=calloc(maxpts,sizeof(float));
    qdat=calloc(maxpts,sizeof(float));
    
    //ptr_to_infile=argv[1];

    FILE *fwsprd, *ftimer;
    strcpy(spots_fname,".");
    strcpy(timer_fname,".");

    if(data_dir != NULL) {
        strcpy(spots_fname,data_dir);
        strcpy(timer_fname,data_dir);
    }

    strncat(spots_fname,"/wspr_spots.txt",20);
    strncat(timer_fname,"/wspr_timer.out",20);
    
    for(i=0; i<256; i++) {
        mettab[0][i]=round( 10*(metric_tables[2][i]-bias) );
        mettab[1][i]=round( 10*(metric_tables[2][255-i]-bias) );
    }


    fwsprd=fopen(spots_fname,"w");
    //  FILE *fdiag;
    //  fdiag=fopen("wsprd_diag","a");
    
    /*
    if((ftimer=fopen(timer_fname,"r"))) {
        //Accumulate timing data
        nr=fscanf(ftimer,"%f %f %f %f %f %f %f",
                  &treadwav,&tcandidates,&tsync0,&tsync1,&tsync2,&tfano,&ttotal);
        fclose(ftimer);
    }
    */
    ftimer=fopen(timer_fname,"w");
    
    //    ptr_to_infile_suffix=strstr(ptr_to_infile,".wav");
        
    t0 = clock();
    npoints=readwavfile(infile,idat, qdat);
    treadwav += (float)(clock()-t0)/CLOCKS_PER_SEC;
    
    if( npoints == 1 ) {
      return 1;
    }
    dialfreq=dialfreq_cmdline - (dialfreq_error*1.0e-06);


    // Parse date and time from given filename
    //    strncpy(date,ptr_to_infile_suffix-11,6);
    //    strncpy(uttime,ptr_to_infile_suffix-4,4);
    //    date[6]='\0';
    //    uttime[4]='\0';

    // Do windowed ffts over 2 symbols, stepped by half symbols
    int nffts=4*floor(npoints/512)-1;

    fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512);
    fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512);
    PLAN3 = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE);
    
    float ps[512][nffts];
    float w[512];
    for(i=0; i<512; i++) {
        w[i]=sin(0.006147931*i);
    }
    

    //*************** main loop starts here *****************

    nblocksize=1;
    maxdrift=4;
    minsync2=0.12;

    // nffts=359
    //1 blok = 256
    //2 blok = 512
    // 1/2 blok = 128
    

    for (i=0; i<nffts; i++) {

      for(j=0; j<512; j++ ) {
	k=i*128+j;
	fftin[j][0]=idat[k] * w[j];
	fftin[j][1]=qdat[k] * w[j];
      }
      fftwf_execute(PLAN3);
      
      /* spectrum swap   */
      for (j=0; j<512; j++ ) {
	k=j+256;
	if( k>511 )
	  k=k-512;
	ps[j][i]=fftout[k][0]*fftout[k][0]+fftout[k][1]*fftout[k][1];
      }

    }


    // hitung  average spektrum
    //	double logpsavg;

        for (i=0; i<512; i++) psavg[i]=0.0;
        for (i=0; i<nffts; i++) {
            for (j=0; j<512; j++) {
                psavg[j]=psavg[j]+ps[j][i];
            }
	}

	/***  Smoothing  dengan  7-point window dan dibatasi spectrum to +/-150 Hz ***/
        int window[7]={1,1,1,1,1,1,1};
        float smspec[411];

        for (i=0; i<411; i++) {
            smspec[i]=0.0;
            for(j=-3; j<=3; j++) {
                k=256-205+i+j;
		smspec[i]=smspec[i]+(window[j+3]*psavg[k]);
            }
	}

        /**** urutkan nilai  spektrum, kemudian ambil level noise dalam betuk prosentase  **/
        float tmpsort[411];
        for (j=0; j<411; j++) {
            tmpsort[j]=smspec[j];
        }
        qsort(tmpsort, 411, sizeof(float), floatcomp);

        /*** estimasi level Noise dari spektrum adalah prosentase harga ke 123/411= 30'th ***/
        float noise_level = tmpsort[122];
        /*************************************************************************/        
        /* Renormalize spectrum so that (large) peaks represent an estimate of snr.
         * We know from experience that threshold snr is near -7dB (-6dB..?) in wspr bandwidth,
         * corresponding to -7-26.3=-33.3dB in 2500 Hz (3000Hz ..?)bandwidth.
         * The corresponding threshold is -42.3 dB in 2500 Hz bandwidth for WSPR-15. */
        
        float min_snr, snr_scaling_factor;
	min_snr = pow(10.0,-7.0/10.0); //this is min snr in wspr bw
	//        min_snr = pow(10.0,-8.0/10.0); //this is min snr in wspr bw

	//	snr_scaling_factor=26.3; /* angka heuristic asli program untuk BW 2500 */
	snr_scaling_factor=27.3;     /* angka heuristic untuk BW 3000 hz */

        for (j=0; j<411; j++) {
	  smspec[j]=smspec[j]/noise_level - 1.0;
            if( smspec[j] < min_snr) smspec[j]=0.1*min_snr;
            continue;
        }
        
        /*** mencari lokal maksimum dalam smoothed spectrum.  ***/
	/*** menghitung SNR untuk setiap lokal maksimum       ***/

        for (i=0; i<200; i++) {
            freq0[i]=0.0;
            snr0[i]=0.0;
            drift0[i]=0.0;
            shift0[i]=0;
            sync0[i]=0.0;
        }
        
        /* min_snr = pow(10.0,-8.0/10.0);  */
        int npk=0;
        unsigned char candidate;
	for(j=1; j<410; j++) {
	  candidate = (smspec[j]>smspec[j-1]) &&
	    (smspec[j]>smspec[j+1]) &&
	    (npk<200);
	  if ( candidate ) {
	    freq0[npk]=(j-205)*df; /* freq0[npk] terhadap 150 Hz (i=205) kekiri(-) 
				      atau kekanan (+) */
	    snr0[npk]=10*log10(smspec[j])-snr_scaling_factor;
	    npk++;
	  }
	}
	
	/*********************************************************/

        /*** jangan buang waktu untuk sinyal diluar rentang [fmin,fmax].  ***/
	/*              fmin = -110 Hz; fmax = +110 Hz                      */
	printf("\n\nnpk= %d\n",npk);

        i=0;
	for( j=0; j<npk; j++) {
	  if( freq0[j] >= fmin && freq0[j] <= fmax ) {
	    freq0[i]=freq0[j];
	    snr0[i]=snr0[j];
	    i++;
	  }
        }
        npk=i;
        
        /** sortir berdasarkan snr, selebar wilayah frekuensi observasi  ***/
        int pass;
        float tmp;
        for (pass = 1; pass <= npk - 1; pass++) {
            for (k = 0; k < npk - pass ; k++) {
                if (snr0[k] < snr0[k+1]) {
                    tmp = snr0[k];
                    snr0[k] = snr0[k+1];
                    snr0[k+1] = tmp;
                    tmp = freq0[k];
                    freq0[k] = freq0[k+1];
                    freq0[k+1] = tmp;
                }
	    }
	}
        t0=clock();


        /* Make coarse estimates of shift (DT), freq, and drift
         
         * Look for time offsets up to +/- 8 symbols (about +/- 5.4 s) relative
         to nominal start time, which is 2 seconds into the file
         
         * Calculates shift relative to the beginning of the file
         
         * Negative shifts mean that signal started before start of file
         
         * The program prints DT = shift-2 s
         
         * Shifts that cause sync vector to fall off of either end of the data
         vector are accommodated by "partial decoding", such that missing
         symbols produce a soft-decision symbol value of 128
         
         * The frequency drift model is linear, deviation of +/- drift/2 over the
         span of 162 symbols, with deviation equal to 0 at the center of the
         signal vector.
         */

        int idrift,ifr,if0,ifd,k0;
        int kindex;
        float smax,ss,pow,p0,p1,p2,p3;

        for(j=0; j<npk; j++) {                              /* untuk semua kandidat frekuensi */
            smax=-1e30;
            if0=freq0[j]/df+256;
            for (ifr=if0-2; ifr<=if0+2; ifr++) {                      //Freq search
                for( k0=-10; k0<22; k0++) {                             //Time search
                    for (idrift=-maxdrift; idrift<=maxdrift; idrift++) {  //Drift search
                        ss=0.0;
                        pow=0.0;
                         for (k=0; k<162; k++) {                             //Sum over symbols
                            ifd=ifr+((float)k-81.0)/81.0*( (float)idrift )/(2.0*df);
                            kindex=k0+2*k;
                            if( kindex < nffts ) {
                                p0=ps[ifd-3][kindex];
                                p1=ps[ifd-1][kindex];
                                p2=ps[ifd+1][kindex];
                                p3=ps[ifd+3][kindex];
                                
                                p0=sqrt(p0);
                                p1=sqrt(p1);
                                p2=sqrt(p2);
                                p3=sqrt(p3);
                                
                                ss=ss+(2*pr3[k]-1)*((p1+p3)-(p0+p2));
                                pow=pow+p0+p1+p2+p3;
                            }
			 }

                        sync1=ss/pow;

                        if( sync1 > smax ) {                  /** simpan parameter kasar **/
                            smax=sync1;
                            shift0[j]=128*(k0+1);
                            drift0[j]=idrift;
                            freq0[j]=(ifr-256)*df;
                            sync0[j]=sync1;
			}
		    }
                }
            }
	}
        tcandidates += (float)(clock()-t0)/CLOCKS_PER_SEC;
	//	tcandidates = (float)(clock()-t0)/CLOCKS_PER_SEC;

        /*
         Rafinasi dari estimasi freq, shift dengan mamakai sync sebagai ukuran.
         Sync dihitung dengan rentang harga mengambang [0.0,1.0].
         
         Fungsi sync_and_demodulate memiliki 3 moda operasi
         moda ada pada argumen terakhir :
         
         0 = tanpa pencarian frequensi atau drift. Hanya mencari ingsut waktu (lag).
         1 = tanpa pencarian ingsut waktu (lag) atau  drift frekuensi. 
	 Hanya mencari frequensi terbaik.
         */

        for (j=0; j<npk; j++){
	  printf("frek[%d]= %f; shift[%d]= %d; sync[%d]= %f\n",j,freq0[j],j,shift0[j],j, sync0[j]);
	}


        for (j=0; j<npk; j++) {
	  
	  printf("\nKANDIDAT FREK KE: %d [%f];\n",j+1,freq0[j]);

            memset(symbols,0,sizeof(char)*nbits*2);
            memset(callsign,0,sizeof(char)*13);
            memset(call_loc_pow,0,sizeof(char)*23);

            f1=freq0[j];
            drift1=drift0[j];
            shift1=shift0[j];
            sync1=sync0[j];

            /* pencarian-kasar ingsut waktu (lag) dan frekuensi, 
	       kemudian sync>minsync1 continue **/
            fstep=0.0; ifmin=0; ifmax=0;
            lagmin=shift1-128;
            lagmax=shift1+128;
            lagstep=64;
            t0 = clock();
            sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
                                lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0);
            tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC;
	    //            tsync0 = (float)(clock()-t0)/CLOCKS_PER_SEC;

            fstep=0.25; ifmin=-2; ifmax=2;
            t0 = clock();
            sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
                                lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1);

	      /** rafinasi estimasi frekuensi  drift **/
                fstep=0.0; ifmin=0; ifmax=0;
                float driftp,driftm,syncp,syncm;
                driftp=drift1+0.5;
                sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
                                lagmin, lagmax, lagstep, &driftp, symfac, &syncp, 1);

                driftm=drift1-0.5;
                sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
                                lagmin, lagmax, lagstep, &driftm, symfac, &syncm, 1);
            

                if(syncp>sync1) {
                    drift1=driftp;
                    sync1=syncp;
                } else if (syncm>sync1) {
                    drift1=driftm;
                    sync1=syncm;
                }


            tsync1 += (float)(clock()-t0)/CLOCKS_PER_SEC;
	    //            tsync1 = (float)(clock()-t0)/CLOCKS_PER_SEC;

	    /*************************************************/
            /** pencarian-halus frekuensi dan ingsut waktu  **/
	    /*************************************************/

            if( sync1 > minsync1 ) {         //minsync1=0.1
        
                lagmin=shift1-32; lagmax=shift1+32; lagstep=16;
                t0 = clock();
                sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
                                    lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0);

		tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC;
		//		                tsync0 = (float)(clock()-t0)/CLOCKS_PER_SEC;
            
                /* pencarian halus  frequensi  **/
                fstep=0.05; ifmin=-2; ifmax=2;
                t0 = clock();
                sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
                                lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1);

                tsync1 += (float)(clock()-t0)/CLOCKS_PER_SEC;
		//            tsync1 = (float)(clock()-t0)/CLOCKS_PER_SEC;
		
                worth_a_try = 1;
            }          /* apakah masih bisa lebih besar korelasinya */
	    else{
	      worth_a_try = 0;
	    }
            
	    // ada update harga shift1 dan drift1

            int idt, ii, jittered_shift;
            float y,sq,rms;
            not_decoded=1;
            int blocksize;  
	    /* iifac = 8 --> Step size in final DT peakup **/
	    /* nblocksize = 1 */

                blocksize=nblocksize;
                idt=0; ii=0;
		int FinalStep = 128/iifac;         /* FinalStep=16 */
                while ( worth_a_try && not_decoded && idt<=FinalStep) {
                    ii=(idt+1)/2;
                    if( idt%2 == 1 ) ii=-ii;   /* jika idt ganjil */
                    ii=iifac*ii;
                    jittered_shift=shift1+ii;
		    //		    printf("WAT= %d; ND= %d; idt=%d ii=%d jitter_shift=%d\n", 
		    //			   worth_a_try,not_decoded,idt,ii,jittered_shift);
		    //		    printf("idt=%d ii=%d jitter_shift=%d\n", idt,ii,jittered_shift);

                    t0 = clock();
                    noncoherent_sequence_detection(idat, qdat, npoints, symbols, &f1,
                                    &jittered_shift, &drift1, symfac, &blocksize);

                    tsync2 += (float)(clock()-t0)/CLOCKS_PER_SEC;
		    //                    tsync2 = (float)(clock()-t0)/CLOCKS_PER_SEC;
                
		    /* menghitung nilai RMS  */
                    sq=0.0;
                    for(i=0; i<162; i++) {
                        y=(float)symbols[i] - 128.0;
                        sq += y*y;
                    }
                    rms=sqrt(sq/162.0);
		    /*************************/
		    
                    if((sync1 > minsync2) && (rms > minrms)) {    // minsync2=0.12
                        deinterleave(symbols);
                        t0 = clock();

			not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
                                           mettab,delta,maxcycles, &upaya);

			if(not_decoded==0){
			  lupaya[uniques]=upaya;
			}
			tfano += (float)(clock()-t0)/CLOCKS_PER_SEC;
                    
		    }
                    idt++;
		}
           
		if( worth_a_try && !not_decoded ) {  // worth_a_try=1 and !not_decoded=1

		  for(i=0; i<11; i++) {
		      message[i]=decdata[i];
		  }


		  // sanity checks on grid and power, and return
		  // call_loc_pow string and also callsign (for de-duping).
		  noprint=unpk_(message,call_loc_pow,callsign);

                // Remove dupes (same callsign and freq within 3 Hz)
                int dupe=0;
                for (i=0; i<uniques; i++) {
                    if(!strcmp(callsign,allcalls[i]) &&
                       (fabs(f1-allfreqs[i]) <3.0)) dupe=1;
                }
                if( !dupe && !noprint) {
		  strcpy(allcalls[uniques],callsign);
		  allfreqs[uniques]=f1;
		  uniques++;
                    
		  // Add an extra space at the end of each line so that wspr-x doesn't
		  // truncate the power (TNX to DL8FCL!)
                    

		  freq_print=dialfreq+(1500+f1)/1e6; //dialfreq == 0
		  dt_print=shift1*dt-1.0;

		  t = time(NULL);
		  tm = *localtime(&t);
		  sprintf(date,"%d:%d:%d", tm.tm_year+1900, tm.tm_mon + 1, tm.tm_mday);
		  sprintf(uttime,"%d:%d", tm.tm_hour, tm.tm_min);
		  
		  strcpy(decodes[uniques-1].date,date);
		  strcpy(decodes[uniques-1].time,uttime);
		  decodes[uniques-1].sync=sync1;
		  decodes[uniques-1].snr=snr0[j];
		  decodes[uniques-1].dt=dt_print;
		  decodes[uniques-1].freq=freq_print;
		  strcpy(decodes[uniques-1].message,call_loc_pow);
		  decodes[uniques-1].drift=drift1;
		  decodes[uniques-1].cycles=cycles;
		  decodes[uniques-1].jitter=ii;
		  decodes[uniques-1].blocksize=blocksize;
		  decodes[uniques-1].metric=metric;
		}
		}  // worth a try dan !not decoded

		//jika dekode sukses maka cari kandidat yang lain
	}  // batas j --> npk; jika segala upaya gagal maka cari kandidat lain
        

    // urutkan hasil berdasarkan urutan frekuensi 
    struct result temp;
    for (j = 1; j <= uniques - 1; j++) {
        for (k = 0; k < uniques - j ; k++) {
            if (decodes[k].freq > decodes[k+1].freq) {
                temp = decodes[k];
                decodes[k]=decodes[k+1];;
                decodes[k+1] = temp;
            }
        }
    }
    
    for (i=0; i<uniques; i++) {
      printf("\n%4s %3.0f %4.1f %10.6f %2d  %-s (jumlah upaya= %d)\n",
	     decodes[i].time, decodes[i].snr,decodes[i].dt, decodes[i].freq,
	     (int)decodes[i].drift, decodes[i].message, lupaya[i]);

        fprintf(fwsprd,
                "%6s %4s %3d %3.0f %4.1f %10.6f  %-22s %2d %5u %4d\n",
                decodes[i].date, decodes[i].time, (int)(10*decodes[i].sync),
                decodes[i].snr, decodes[i].dt, decodes[i].freq,
                decodes[i].message, (int)decodes[i].drift, decodes[i].cycles/81,
                decodes[i].jitter);
        
    }
    printf("<Dekode Paripurna...>\n");
    
    fftwf_free(fftin);
    fftwf_free(fftout);
    
    ttotal += (float)(clock()-t00)/CLOCKS_PER_SEC;
    //    ttotal = (float)(clock()-t00)/CLOCKS_PER_SEC;

    fprintf(ftimer,"%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n\n",
            treadwav,tcandidates,tsync0,tsync1,tsync2,tfano,ttotal);
    
    fprintf(ftimer,"Code segment        Seconds   Frac\n");
    fprintf(ftimer,"-----------------------------------\n");
    fprintf(ftimer,"readwavfile        %7.2f %7.2f\n",treadwav,treadwav/ttotal);
    fprintf(ftimer,"Coarse DT f0 f1    %7.2f %7.2f\n",tcandidates,
            tcandidates/ttotal);
    fprintf(ftimer,"sync_and_demod(0)  %7.2f %7.2f\n",tsync0,tsync0/ttotal);
    fprintf(ftimer,"sync_and_demod(1)  %7.2f %7.2f\n",tsync1,tsync1/ttotal);
    fprintf(ftimer,"sync_and_demod(2)  %7.2f %7.2f\n",tsync2,tsync2/ttotal);
    fprintf(ftimer,"Stack/Fano decoder %7.2f %7.2f\n",tfano,tfano/ttotal);
    fprintf(ftimer,"-----------------------------------\n");
    fprintf(ftimer,"Total              %7.2f %7.2f\n",ttotal,1.0);
    
    fclose(fwsprd);
    fclose(ftimer);
    fftwf_destroy_plan(PLAN1);
    fftwf_destroy_plan(PLAN2);
    fftwf_destroy_plan(PLAN3);
    

   
    free(symbols);
    free(decdata);
    free(channel_symbols);
    free(callsign);
    free(call_loc_pow);
    free(idat);
    free(qdat); 

    
    if(writenoise == 999) return -1;  //Silence compiler warning
    return 0;
}