From 57ef481004abb303dd533e3043e34b3b53167329 Mon Sep 17 00:00:00 2001
From: Steven Franke <s.j.franke@icloud.com>
Date: Sun, 14 Jun 2015 20:24:52 +0000
Subject: [PATCH] First steps toward 2-pass decoder that attempts to subtract
 strong signals.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@5611 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
---
 lib/wsprd/wsprd_exp.c | 1026 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 1026 insertions(+)
 create mode 100644 lib/wsprd/wsprd_exp.c

diff --git a/lib/wsprd/wsprd_exp.c b/lib/wsprd/wsprd_exp.c
new file mode 100644
index 000000000..02fb3e10b
--- /dev/null
+++ b/lib/wsprd/wsprd_exp.c
@@ -0,0 +1,1026 @@
+/*
+ This file is part of program wsprd, a detector/demodulator/decoder
+ for the Weak Signal Propagation Reporter (WSPR) mode.
+ 
+ File name: wsprd.c
+ 
+ Copyright 2001-2015, Joe Taylor, K1JT
+ 
+ Much of the present code is based on work by Steven Franke, K9AN,
+ which in turn was based on earlier work by K1JT.
+ 
+ Copyright 2014-2015, Steven Franke, K9AN
+ 
+ License: GNU GPL v3
+ 
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#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 "fano.h"
+#include "nhash.h"
+#include "wsprd_utils.h"
+#include "wsprsim_utils.h"
+
+#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
+fftw_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;
+
+//***************************************************************************
+unsigned long readc2file(char *ptr_to_infile, double *idat, double *qdat,
+                         double *freq, int *wspr_type)
+{
+    float buffer[2*65536];
+    double dfreq;
+    int i,ntrmin;
+    char *c2file[15];
+    FILE* fp;
+    
+    fp = fopen(ptr_to_infile,"rb");
+    if (fp == NULL) {
+        fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile);
+        return 1;
+    }
+    unsigned long nread=fread(c2file,sizeof(char),14,fp);
+    nread=fread(&ntrmin,sizeof(int),1,fp);
+    nread=fread(&dfreq,sizeof(double),1,fp);
+    *freq=dfreq;
+    nread=fread(buffer,sizeof(float),2*45000,fp);
+    
+    *wspr_type=ntrmin;
+    
+    for(i=0; i<45000; i++) {
+        idat[i]=buffer[2*i];
+        qdat[i]=-buffer[2*i+1];
+    }
+    
+    if( nread == 2*45000 ) {
+        return nread/2;
+    } else {
+        return 1;
+    }
+}
+
+//***************************************************************************
+unsigned long readwavfile(char *ptr_to_infile, int ntrmin, double *idat, double *qdat )
+{
+    unsigned long i, j, npoints;
+    int nfft1, nfft2, nh2, i0;
+    double df;
+    
+    nfft2=46080; //this is the number of downsampled points that will be returned
+    nh2=nfft2/2;
+    
+    if( ntrmin == 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;
+    } else if ( ntrmin == 15 ) {
+        nfft1=nfft2*8*32;
+        df=12000.0/nfft1;
+        i0=(1500.0+112.5)/df+0.5;
+        npoints=8*114*12000;
+    } else {
+        fprintf(stderr,"This should not happen\n");
+        return 1;
+    }
+    
+    double *realin;
+    fftw_complex *fftin, *fftout;
+    
+    FILE *fp;
+    short int *buf2;
+    buf2 = malloc(npoints*sizeof(short int));
+    
+    fp = fopen(ptr_to_infile,"rb");
+    if (fp == NULL) {
+        fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile);
+        return 1;
+    }
+    nr=fread(buf2,2,22,fp);            //Read and ignore header
+    nr=fread(buf2,2,npoints,fp);       //Read raw data
+    fclose(fp);
+    
+    realin=(double*) fftw_malloc(sizeof(double)*nfft1);
+    fftout=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*nfft1);
+    PLAN1 = fftw_plan_dft_r2c_1d(nfft1, realin, fftout, PATIENCE);
+    
+    for (i=0; i<npoints; i++) {
+        realin[i]=buf2[i]/32768.0;
+    }
+    
+    for (i=npoints; i<nfft1; i++) {
+        realin[i]=0.0;
+    }
+    
+    free(buf2);
+    fftw_execute(PLAN1);
+    fftw_free(realin);
+    
+    fftin=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*nfft2);
+    
+    for (i=0; i<nfft2; i++) {
+        j=i0+i;
+        if( i>nh2 ) j=j-nfft2;
+        fftin[i][0]=fftout[j][0];
+        fftin[i][1]=fftout[j][1];
+    }
+    
+    fftw_free(fftout);
+    fftout=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*nfft2);
+    PLAN2 = fftw_plan_dft_1d(nfft2, fftin, fftout, FFTW_BACKWARD, PATIENCE);
+    fftw_execute(PLAN2);
+    
+    for (i=0; i<nfft2; i++) {
+        idat[i]=fftout[i][0]/1000.0;
+        qdat[i]=fftout[i][1]/1000.0;
+    }
+    
+    fftw_free(fftin);
+    fftw_free(fftout);
+    return nfft2;
+}
+
+//***************************************************************************
+void sync_and_demodulate(double *id, double *qd, long np,
+                         unsigned char *symbols, float *f1, float fstep,
+                         int *shift1, int lagmin, int lagmax, int lagstep,
+                         float *drift1, int symfac, float *sync, int mode)
+{
+    /***********************************************************************
+     * mode = 0: no frequency or drift search. find best time lag.          *
+     *        1: no time lag or drift search. find best frequency.          *
+     *        2: no frequency or time lag search. calculate soft-decision   *
+     *           symbols using passed frequency and shift.                  *
+     ************************************************************************/
+    
+    float dt=1.0/375.0, df=375.0/256.0, fbest=0.0;
+    int i, j, k;
+    double pi=4.*atan(1.0),twopidt;
+    float f0=0.0,fp,fplast=-10000.0,ss;
+    int lag;
+    
+    double i0[162],q0[162],i1[162],q1[162],i2[162],q2[162],i3[162],q3[162];
+    double p0,p1,p2,p3,cmet,totp,syncmax,fac;
+    double c0[256],s0[256],c1[256],s1[256],c2[256],s2[256],c3[256],s3[256];
+    double dphi0, cdphi0, sdphi0, dphi1, cdphi1, sdphi1, dphi2, cdphi2, sdphi2,
+    dphi3, cdphi3, sdphi3;
+    float fsum=0.0, f2sum=0.0, fsymb[162];
+    int best_shift = 0, ifreq;
+    int ifmin=0, ifmax=0;
+    
+    syncmax=-1e30;
+    if( mode == 0 ) {ifmin=0; ifmax=0; fstep=0.0; f0=*f1;}
+    if( mode == 1 ) {lagmin=*shift1;lagmax=*shift1;ifmin=-5;ifmax=5;f0=*f1;}
+    if( mode == 2 ) {lagmin=*shift1;lagmax=*shift1;ifmin=0;ifmax=0;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 + ((float)*drift1/2.0)*((float)i-81.0)/81.0;
+                if( i==0 || (fp != fplast) ) {  // only calculate sin/cos if necessary
+                    dphi0=2*pi*(fp-1.5*df)*dt;
+                    cdphi0=cos(dphi0);
+                    sdphi0=sin(dphi0);
+                    
+                    dphi1=twopidt*(fp-0.5*df);
+                    cdphi1=cos(dphi1);
+                    sdphi1=sin(dphi1);
+                    
+                    dphi2=twopidt*(fp+0.5*df);
+                    cdphi2=cos(dphi2);
+                    sdphi2=sin(dphi2);
+                    
+                    dphi3=twopidt*(fp+1.5*df);
+                    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;
+                    }
+                    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=ss+cmet*(2*pr3[i]-1);
+                if( mode == 2) {                 //Compute soft symbols
+                    if(pr3[i]) {
+                        fsymb[i]=p3-p1;
+                    } else {
+                        fsymb[i]=p2-p0;
+                    }
+                }
+            }
+            
+            if( ss/totp > syncmax ) {          //Save best parameters
+                syncmax=ss/totp;
+                best_shift=lag;
+                fbest=f0;
+            }
+        } // lag loop
+    } //freq loop
+    
+    if( mode <=1 ) {                       //Send best params back to caller
+        *sync=syncmax;
+        *shift1=best_shift;
+        *f1=fbest;
+        return;
+    }
+    
+    if( mode == 2 ) {
+        *sync=syncmax;
+        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;
+        }
+        return;
+    }
+    return;
+}
+//***************************************************************************
+void subtract_signal(double *id, double *qd, long np,
+                     float f0, int shift0, float drift0, unsigned char* channel_symbols)
+{
+    float dt=1.0/375.0, df=375.0/256.0;
+    int i, j, k;
+    double pi=4.*atan(1.0),twopidt, fp;
+    
+    double i0[162],q0[162];
+    double c0[256],s0[256];
+    double dphi, cdphi, sdphi;
+    
+    twopidt=2*pi*dt;
+    
+    printf("subtract_signal f0: %f shift0: %d drift0: %f\n",f0,shift0,drift0);
+    for (i=0; i<162; i++) {
+        fp = f0 + ((float)drift0/2.0)*((float)i-81.0)/81.0;
+        
+            dphi=twopidt*(fp+((float)channel_symbols[i]-1.5)*df);
+            cdphi=cos(dphi);
+            sdphi=sin(dphi);
+        
+            c0[0]=1; s0[0]=0;
+        
+            for (j=1; j<256; j++) {
+                c0[j]=c0[j-1]*cdphi - s0[j-1]*sdphi;
+                s0[j]=c0[j-1]*sdphi + s0[j-1]*cdphi;
+            }
+        
+        i0[i]=0.0; q0[i]=0.0;
+        
+        for (j=0; j<256; j++) {
+            k=shift0+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];
+            }
+        }
+
+// subtract the signal here.
+        
+    }
+    
+    return;
+}
+
+unsigned long writec2file(char *c2filename, int trmin, double freq
+                          , double *idat, double *qdat)
+{
+    int i;
+    float buffer[2*45000];
+    memset(buffer,0,sizeof(float)*2*45000);
+    FILE *fp;
+    
+    fp = fopen(c2filename,"wb");
+    if( fp == NULL ) {
+        fprintf(stderr, "Could not open c2 file '%s'\n", c2filename);
+        return 0;
+    }
+    unsigned long nwrite = fwrite(c2filename,sizeof(char),14,fp);
+    nwrite = fwrite(&trmin, sizeof(int), 1, fp);
+    nwrite = fwrite(&freq, sizeof(double), 1, fp);
+    
+    for(i=0; i<45000; i++) {
+        buffer[2*i]=idat[i];
+        buffer[2*i+1]=-qdat[i];
+    }
+    
+    nwrite = fwrite(buffer, sizeof(float), 2*45000, fp);
+    if( nwrite == 2*45000 ) {
+        return nwrite;
+    } else {
+        return 0;
+    }
+}
+
+//***************************************************************************
+void usage(void)
+{
+    printf("Usage: wsprd [options...] infile\n");
+    printf("       infile must have suffix .wav or .c2\n");
+    printf("\n");
+    printf("Options:\n");
+    printf("       -a <path> path to writeable data files, default=\".\"\n");
+    printf("       -e x (x is transceiver dial frequency error in Hz)\n");
+    printf("       -f x (x is transceiver dial frequency in MHz)\n");
+    // blanking is not yet implemented. The options are accepted for compatibility
+    // with development version of wsprd.
+    //    printf("       -t n (n is blanking duration in milliseconds)\n");
+    //    printf("       -b n (n is pct of time that is blanked)\n");
+    printf("       -H do not use (or update) the hash table\n");
+    printf("       -m decode wspr-15 .wav file\n");
+    printf("       -n write noise estimates to file noise.dat\n");
+    printf("       -q quick mode - doesn't dig deep for weak signals\n");
+    printf("       -s slow mode - much slower, yields a few more decodes\n");
+    printf("       -v verbose mode\n");
+    printf("       -w wideband mode - decode signals within +/- 150 Hz of center\n");
+    printf("       -z x (x is fano metric table bias, default is 0.42)\n");
+}
+
+//***************************************************************************
+int main(int argc, char *argv[])
+{
+    extern char *optarg;
+    extern int optind;
+    int i,j,k;
+    unsigned char *symbols, *decdata;
+    signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0};
+    char *callsign,*grid,*grid6, *call_loc_pow, *cdbm;
+    char *ptr_to_infile,*ptr_to_infile_suffix;
+    char *data_dir=NULL;
+    char wisdom_fname[200],all_fname[200],spots_fname[200];
+    char timer_fname[200],hash_fname[200];
+    char uttime[5],date[7];
+    int c,delta,maxpts=65536,verbose=0,quickmode=0,writenoise=0,usehashtable=1,wspr_type=2;
+    int shift1, lagmin, lagmax, lagstep, worth_a_try, not_decoded;
+    unsigned int nbits;
+    unsigned int npoints, metric, maxcycles, 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;
+    float dialfreq_error=0.0;
+    float fmin=-110, fmax=110;
+    float f1, fstep, sync1, drift1, tblank=0, fblank=0;
+    double *idat, *qdat;
+    clock_t t0,t00;
+    double tfano=0.0,treadwav=0.0,tcandidates=0.0,tsync0=0.0;
+    double tsync1=0.0,tsync2=0.0,ttotal=0.0;
+    
+    // Parameters used for performance-tuning:
+    maxcycles=10000;                         //Fano timeout limit
+    double minsync1=0.10;                    //First sync limit
+    double minsync2=0.12;                    //Second sync limit
+    int iifac=3;                             //Step size in final DT peakup
+    int symfac=50;                           //Soft-symbol normalizing factor
+    int maxdrift=4;                          //Maximum (+/-) drift
+    double minrms=52.0 * (symfac/64.0);      //Final test for plausible decoding
+    delta=60;                                //Fano threshold step
+    
+    t00=clock();
+    fftw_complex *fftin, *fftout;
+#include "./metric_tables.c"
+    
+    int mettab[2][256];
+    float bias=0.42;
+    
+    idat=malloc(sizeof(double)*maxpts);
+    qdat=malloc(sizeof(double)*maxpts);
+    
+    while ( (c = getopt(argc, argv, "a:b:e:f:Hmnqst:wvz:")) !=-1 ) {
+        switch (c) {
+            case 'a':
+                data_dir = optarg;
+                break;
+            case 'b':
+                fblank = strtof(optarg,NULL);
+                break;
+            case 'e':
+                dialfreq_error = strtof(optarg,NULL);   // units of Hz
+                // dialfreq_error = dial reading - actual, correct frequency
+                break;
+            case 'f':
+                dialfreq_cmdline = strtod(optarg,NULL); // units of MHz
+                break;
+            case 'H':
+                usehashtable = 0;
+                break;
+            case 'm':
+                wspr_type = 15;
+                break;
+            case 'n':
+                writenoise = 1;
+                break;
+            case 'q':
+                quickmode = 1;
+                break;
+            case 's':
+                maxcycles=20000;
+                iifac=1;
+                break;
+            case 't':
+                tblank = strtof(optarg,NULL);
+                break;
+            case 'v':
+                verbose = 1;
+                break;
+            case 'w':
+                fmin=-150.0;
+                fmax=150.0;
+                break;
+            case 'z':
+                bias=strtof(optarg,NULL); //fano metric bias (default is 0.42)
+                break;
+            case '?':
+                usage();
+                return 1;
+        }
+    }
+    
+    if( optind+1 > argc) {
+        usage();
+        return 1;
+    } else {
+        ptr_to_infile=argv[optind];
+    }
+    
+    // setup metric table
+    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) );
+    }
+    
+    FILE *fp_fftw_wisdom_file, *fall_wspr, *fwsprd, *fhash, *ftimer;
+    strcpy(wisdom_fname,".");
+    strcpy(all_fname,".");
+    strcpy(spots_fname,".");
+    strcpy(timer_fname,".");
+    strcpy(hash_fname,".");
+    if(data_dir != NULL) {
+        strcpy(wisdom_fname,data_dir);
+        strcpy(all_fname,data_dir);
+        strcpy(spots_fname,data_dir);
+        strcpy(timer_fname,data_dir);
+        strcpy(hash_fname,data_dir);
+    }
+    strncat(wisdom_fname,"/wspr_wisdom.dat",20);
+    strncat(all_fname,"/ALL_WSPR.TXT",20);
+    strncat(spots_fname,"/wspr_spots.txt",20);
+    strncat(timer_fname,"/wspr_timer.out",20);
+    strncat(hash_fname,"/hashtable.txt",20);
+    if ((fp_fftw_wisdom_file = fopen(wisdom_fname, "r"))) {  //Open FFTW wisdom
+        fftw_import_wisdom_from_file(fp_fftw_wisdom_file);
+        fclose(fp_fftw_wisdom_file);
+    }
+    
+    fall_wspr=fopen(all_fname,"a");
+    fwsprd=fopen(spots_fname,"w");
+    //  FILE *fdiag;
+    //  fdiag=fopen("wsprd_diag","a");
+    
+    if((ftimer=fopen(timer_fname,"r"))) {
+        //Accumulate timing data
+        nr=fscanf(ftimer,"%lf %lf %lf %lf %lf %lf %lf",
+                  &treadwav,&tcandidates,&tsync0,&tsync1,&tsync2,&tfano,&ttotal);
+        fclose(ftimer);
+    }
+    ftimer=fopen(timer_fname,"w");
+    
+    if( strstr(ptr_to_infile,".wav") ) {
+        ptr_to_infile_suffix=strstr(ptr_to_infile,".wav");
+        
+        t0 = clock();
+        npoints=readwavfile(ptr_to_infile, wspr_type, idat, qdat);
+        treadwav += (double)(clock()-t0)/CLOCKS_PER_SEC;
+        
+        if( npoints == 1 ) {
+            return 1;
+        }
+        dialfreq=dialfreq_cmdline - (dialfreq_error*1.0e-06);
+    } else if ( strstr(ptr_to_infile,".c2") !=0 )  {
+        ptr_to_infile_suffix=strstr(ptr_to_infile,".c2");
+        npoints=readc2file(ptr_to_infile, idat, qdat, &dialfreq, &wspr_type);
+        if( npoints == 1 ) {
+            return 1;
+        }
+        dialfreq -= (dialfreq_error*1.0e-06);
+    } else {
+        printf("Error: Failed to open %s\n",ptr_to_infile);
+        printf("WSPR file must have suffix .wav or .c2\n");
+        return 1;
+    }
+    
+    // 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=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*512);
+    fftout=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*512);
+    PLAN3 = fftw_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);
+    }
+    
+    memset(ps,0.0, sizeof(float)*512*nffts);
+    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];
+        }
+        fftw_execute(PLAN3);
+        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];
+        }
+    }
+    
+    fftw_free(fftin);
+    fftw_free(fftout);
+    
+    // Compute average spectrum
+    float psavg[512];
+    memset(psavg,0.0, sizeof(float)*512);
+    for (i=0; i<nffts; i++) {
+        for (j=0; j<512; j++) {
+            psavg[j]=psavg[j]+ps[j][i];
+        }
+    }
+    
+    // Smooth with 7-point window and limit 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];
+        }
+    }
+    
+    // Sort spectrum values, then pick off noise level as a percentile
+    float tmpsort[411];
+    for (j=0; j<411; j++) {
+        tmpsort[j]=smspec[j];
+    }
+    qsort(tmpsort, 411, sizeof(float), floatcomp);
+    
+    // Noise level of spectrum is estimated as 123/411= 30'th percentile
+    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 in wspr bandwidth,
+     * corresponding to -7-26.3=-33.3dB in 2500 Hz 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
+    if( wspr_type == 2 ) {
+        snr_scaling_factor=26.3;
+    } else {
+        snr_scaling_factor=35.3;
+    }
+    for (j=0; j<411; j++) {
+        smspec[j]=smspec[j]/noise_level - 1.0;
+        if( smspec[j] < min_snr) smspec[j]=0.1;
+        continue;
+    }
+    
+    // Find all local maxima in smoothed spectrum.
+    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;
+    }
+    
+    int npk=0;
+    for(j=1; j<410; j++) {
+        if((smspec[j]>smspec[j-1]) && (smspec[j]>smspec[j+1]) && (npk<200)) {
+            freq0[npk]=(j-205)*df;
+            snr0[npk]=10*log10(smspec[j])-snr_scaling_factor;
+            npk++;
+        }
+    }
+    
+    // Compute corrected fmin, fmax, accounting for dial frequency error
+    fmin += dialfreq_error;    // dialfreq_error is in units of Hz
+    fmax += dialfreq_error;
+    
+    // Don't waste time on signals outside of the range [fmin,fmax].
+    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;
+    
+    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++) {                              //For each candidate...
+        smax=-1e30;
+        if0=freq0[j]/df+256;
+        for (ifr=if0-1; ifr<=if0+1; 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 ) {                  //Save coarse parameters
+                        smax=sync1;
+                        shift0[j]=128*(k0+1);
+                        drift0[j]=idrift;
+                        freq0[j]=(ifr-256)*df;
+                        sync0[j]=sync1;
+                    }
+                }
+            }
+        }
+    }
+    tcandidates += (double)(clock()-t0)/CLOCKS_PER_SEC;
+    
+    nbits=81;
+    symbols=malloc(sizeof(char)*nbits*2);
+    memset(symbols,0,sizeof(char)*nbits*2);
+    decdata=malloc((nbits+7)/8);
+    grid=malloc(sizeof(char)*5);
+    grid6=malloc(sizeof(char)*7);
+    callsign=malloc(sizeof(char)*13);
+    call_loc_pow=malloc(sizeof(char)*23);
+    cdbm=malloc(sizeof(char)*3);
+    float allfreqs[npk];
+    memset(allfreqs,0,sizeof(float)*npk);
+    char allcalls[npk][13];
+    memset(allcalls,0,sizeof(char)*npk*13);
+    memset(grid,0,sizeof(char)*5);
+    memset(grid6,0,sizeof(char)*7);
+    memset(callsign,0,sizeof(char)*13);
+    memset(call_loc_pow,0,sizeof(char)*23);
+    memset(cdbm,0,sizeof(char)*3);
+    char hashtab[32768][13];
+    memset(hashtab,0,sizeof(char)*32768*13);
+    int nh;
+    
+    if( usehashtable ) {
+        char line[80], hcall[12];
+        if( (fhash=fopen(hash_fname,"r+")) ) {
+            while (fgets(line, sizeof(line), fhash) != NULL) {
+                sscanf(line,"%d %s",&nh,hcall);
+                strcpy(*hashtab+nh*13,hcall);
+            }
+        } else {
+            fhash=fopen(hash_fname,"w+");
+        }
+        fclose(fhash);
+    }
+    
+    int uniques=0, noprint=0;
+    /*
+     Refine the estimates of freq, shift using sync as a metric.
+     Sync is calculated such that it is a float taking values in the range
+     [0.0,1.0].
+     
+     Function sync_and_demodulate has three modes of operation
+     mode is the last argument:
+     
+     0 = no frequency or drift search. find best time lag.
+     1 = no time lag or drift search. find best frequency.
+     2 = no frequency or time lag search. Calculate soft-decision
+     symbols using passed frequency and shift.
+     
+     NB: best possibility for OpenMP may be here: several worker threads
+     could each work on one candidate at a time.
+     */
+    
+    for (j=0; j<npk; j++) {
+        f1=freq0[j];
+        drift1=drift0[j];
+        shift1=shift0[j];
+        sync1=sync0[j];
+        
+        // Fine search for best sync lag (mode 0)
+        fstep=0.0;
+        lagmin=shift1-144;
+        lagmax=shift1+144;
+        lagstep=8;
+        if(quickmode) lagstep=16;
+        t0 = clock();
+        sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep, &shift1,
+                            lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0);
+        tsync0 += (double)(clock()-t0)/CLOCKS_PER_SEC;
+        
+        // Fine search for frequency peak (mode 1)
+        fstep=0.1;
+        t0 = clock();
+        sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep, &shift1,
+                            lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1);
+        tsync1 += (double)(clock()-t0)/CLOCKS_PER_SEC;
+        
+        if( sync1 > minsync1 ) {
+            worth_a_try = 1;
+        } else {
+            worth_a_try = 0;
+        }
+        
+        int idt=0, ii=0, jiggered_shift;
+        double y,sq,rms;
+        not_decoded=1;
+        
+        while ( worth_a_try && not_decoded && idt<=(128/iifac)) {
+            ii=(idt+1)/2;
+            if( idt%2 == 1 ) ii=-ii;
+            ii=iifac*ii;
+            jiggered_shift=shift1+ii;
+            
+            // Use mode 2 to get soft-decision symbols
+            t0 = clock();
+            sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep,
+                                &jiggered_shift, lagmin, lagmax, lagstep, &drift1, symfac,
+                                &sync1, 2);
+            tsync2 += (double)(clock()-t0)/CLOCKS_PER_SEC;
+            
+            sq=0.0;
+            for(i=0; i<162; i++) {
+                y=(double)symbols[i] - 128.0;
+                sq += y*y;
+            }
+            rms=sqrt(sq/162.0);
+            
+            if((sync1 > minsync2) && (rms > minrms)) {
+                deinterleave(symbols);
+                t0 = clock();
+                not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
+                                   mettab,delta,maxcycles);
+                tfano += (double)(clock()-t0)/CLOCKS_PER_SEC;
+                
+                /* ### Used for timing tests:
+                 if(not_decoded) fprintf(fdiag,
+                 "%6s %4s %4.1f %3.0f %4.1f %10.7f  %-18s %2d %5u %4d %6.1f %2d\n",
+                 date,uttime,sync1*10,snr0[j], shift1*dt-2.0, dialfreq+(1500+f1)/1e6,
+                 "@                 ", (int)drift1, cycles/81, ii, rms, maxnp);
+                 */
+            }
+            idt++;
+            if( quickmode ) break;
+        }
+        
+        if( worth_a_try && !not_decoded ) {
+            for(i=0; i<11; i++) {
+                if( decdata[i]>127 ) {
+                    message[i]=decdata[i]-256;
+                } else {
+                    message[i]=decdata[i];
+                }
+            }
+            
+            // Unpack the decoded message, update the hashtable, apply
+            // sanity checks on grid and power, and return
+            // call_loc_pow string and also callsign (for de-duping).
+            noprint=unpk_(message,hashtab,call_loc_pow,callsign);
+
+            unsigned char channel_symbols[162];
+            get_wspr_channel_symbols(call_loc_pow, channel_symbols);
+            
+            subtract_signal(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
+
+            // Remove dupes (same callsign and freq within 1 Hz)
+            int dupe=0;
+            for (i=0; i<npk; i++) {
+                if(!strcmp(callsign,allcalls[i]) && 
+                   (fabs(f1-allfreqs[i]) <1.0)) dupe=1;
+            }
+            if( (verbose || !dupe) && !noprint) {
+                uniques++;
+                strcpy(allcalls[uniques],callsign);
+                allfreqs[uniques]=f1;
+                // Add an extra space at the end of each line so that wspr-x doesn't 
+                // truncate the power (TNX to DL8FCL!)
+                
+                if( wspr_type == 15 ) {
+                    freq_print=dialfreq+(1500+112.5+f1/8.0)/1e6;
+                    dt_print=shift1*8*dt-2.0;
+                } else {
+                    freq_print=dialfreq+(1500+f1)/1e6;
+                    dt_print=shift1*dt-2.0;
+                }
+                printf("%4s %3.0f %4.1f %10.6f %2d  %-s \n",
+                       uttime, snr0[j],dt_print,freq_print,
+                       (int)drift1, call_loc_pow);
+                
+                fprintf(fall_wspr,
+                        "%6s %4s %3.0f %3.0f %4.1f %10.7f  %-22s %2d %5u %4d\n",
+                        date,uttime,sync1*10,snr0[j],
+                        dt_print, freq_print,
+                        call_loc_pow, (int)drift1, cycles/81, ii);
+                
+                fprintf(fwsprd,"%6s %4s %3d %3.0f %4.1f %10.6f  %-22s %2d %5u %4d\n",
+                        date,uttime,(int)(sync1*10),snr0[j],
+                        dt_print, freq_print,
+                        call_loc_pow, (int)drift1, cycles/81, ii);
+                
+                /* For timing tests
+                 
+                 fprintf(fdiag,
+                 "%6s %4s %4.1f %3.0f %4.1f %10.7f  %-18s %2d %5u %4d %6.1f\n",
+                 date,uttime,sync1*10,snr0[j],
+                 shift1*dt-2.0, dialfreq+(1500+f1)/1e6,
+                 call_loc_pow, (int)drift1, cycles/81, ii, rms);
+                 */
+            }
+        }
+    }
+    printf("<DecodeFinished>\n");
+    
+    if ((fp_fftw_wisdom_file = fopen(wisdom_fname, "w"))) {
+        fftw_export_wisdom_to_file(fp_fftw_wisdom_file);
+        fclose(fp_fftw_wisdom_file);
+    }
+    
+    ttotal += (double)(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,"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(fall_wspr);
+    fclose(fwsprd);
+    //  fclose(fdiag);
+    fclose(ftimer);
+    fftw_destroy_plan(PLAN1);
+    fftw_destroy_plan(PLAN2);
+    fftw_destroy_plan(PLAN3);
+    
+    if( usehashtable ) {
+        fhash=fopen(hash_fname,"w");
+        for (i=0; i<32768; i++) {
+            if( strncmp(hashtab[i],"\0",1) != 0 ) {
+                fprintf(fhash,"%5d %s\n",i,*hashtab+i*13);
+            }
+        }
+        fclose(fhash);
+    }
+    
+    char c2filename[15];
+    double carrierfreq=10.1387;
+    int wsprtype=2;
+    strcpy(c2filename,"000000_0001.c2");
+    printf("Writing %s\n",c2filename);
+    writec2file(c2filename, wsprtype, carrierfreq, idat, qdat);
+
+    if(fblank+tblank+writenoise == 999) return -1;  //Silence compiler warning
+    return 0;
+}