diff --git a/lib/wsprd/wsprd_exp.c b/lib/wsprd/wsprd_exp.c
new file mode 100644
index 000000000..f22728b44
--- /dev/null
+++ b/lib/wsprd/wsprd_exp.c
@@ -0,0 +1,1431 @@
+/*
+ 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 .
+ */
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "fano.h"
+#include "jelinek.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
+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;
+
+//***************************************************************************
+unsigned long readc2file(char *ptr_to_infile, float *idat, float *qdat,
+ double *freq, int *wspr_type)
+{
+ float *buffer;
+ double dfreq;
+ int i,ntrmin;
+ char *c2file[15];
+ FILE* fp;
+
+ buffer=malloc(sizeof(float)*2*65536);
+ memset(buffer,0,sizeof(float)*2*65536);
+
+ 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);
+ fclose(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;
+ }
+ free(buffer);
+}
+
+//***************************************************************************
+unsigned long readwavfile(char *ptr_to_infile, int ntrmin, float *idat, float *qdat )
+{
+ size_t 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;
+ }
+
+ float *realin;
+ fftwf_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=(float*) fftwf_malloc(sizeof(float)*nfft1);
+ fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*nfft1);
+ PLAN1 = fftwf_plan_dft_r2c_1d(nfft1, realin, fftout, PATIENCE);
+
+ for (i=0; 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: 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. *
+ ************************************************************************/
+
+ 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,fac;
+ 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, fsum=0.0, f2sum=0.0, fsymb[162];
+ 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;}
+ 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 + (*drift1/2.0)*((float)i-81.0)/81.0;
+ if( i==0 || (fp != fplast) ) { // only calculate sin/cos if necessary
+ 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;
+ }
+ 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 syncmax ) { //Save best parameters
+ syncmax=ss;
+ 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;
+ }
+ printf("a: %f %f %f\n",fsymb[0],fsymb[1],fsymb[2]);
+ printf("a: %f %f %f\n",fsymb[3],fsymb[4],fsymb[5]);
+ printf("a: %f %f %f\n",fsymb[6],fsymb[7],fsymb[8]);
+ return;
+ }
+ return;
+}
+
+void noncoherent_sequence_detection(float *id, float *qd, long np,
+ unsigned char *symbols, float *f1, int *shift1,
+ float *drift1, int symfac, float *sync, int *nblocksize)
+{
+ /************************************************************************
+ * Noncoherent sequence detection for wspr. *
+ * Allowed block lengths are nblock=1,2,3,6, or 9 symbols. *
+ * Longer block lengths require longer channel coherence time. *
+ * The whole block is estimated at once. *
+ * nblock=1 corresponds to noncoherent detection of individual symbols *
+ * like the original wsprd symbol demodulator. *
+ ************************************************************************/
+ 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, k0, isign, itone, ib, b, nblock, nseq, imask;
+ float xi[512],xq[512];
+ float is[4][162],qs[4][162];
+ float p[512],totp,fac,xm1,xm0;
+ 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, fp, ss, fsum=0.0, f2sum=0.0, fsymb[162];
+
+ twopidt=2*pi*dt;
+ f0=*f1;
+ lag=*shift1;
+ nblock=*nblocksize;
+ nseq=1<0) && (k>(nblock-1-ib);
+ itone=pr3[i+ib]+2*b;
+ xi[j]=xi[j]+isign*is[itone][i+ib];
+ xq[j]=xq[j]+isign*qs[itone][i+ib];
+ isign=-isign;
+ }
+ p[j]=xi[j]*xi[j]+xq[j]*xq[j];
+ p[j]=sqrt(p[j]);
+ }
+ for (ib=0; ib xm1) xm1=p[j];
+ }
+ if((j & imask)==0) {
+ if(p[j]>xm0) xm0=p[j];
+ }
+ }
+ fsymb[i+ib]=xm1-xm0;
+ }
+ }
+ 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;
+}
+
+/***************************************************************************
+ symbol-by-symbol signal subtraction
+ ****************************************************************************/
+void subtract_signal(float *id, float *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;
+ float pi=4.*atan(1.0),twopidt, fp;
+
+ float i0,q0;
+ float c0[256],s0[256];
+ float dphi, cdphi, sdphi;
+
+ twopidt=2*pi*dt;
+
+ 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=0.0; q0=0.0;
+
+ for (j=0; j<256; j++) {
+ k=shift0+i*256+j;
+ if( (k>0) & (k0) & (k0) && (k(nsig-1-nfilt/2) ) {
+ norm=partialsum[nfilt/2+nsig-1-i];
+ } else {
+ norm=1.0;
+ }
+ k=shift0+i;
+ j=i+nfilt;
+ if( (k>0) && (k path to writeable data files, default=\".\"\n");
+ printf(" -c write .c2 file at the end of the first pass\n");
+ printf(" -C maximum number of decoder cycles per bit, default 10000\n");
+ printf(" -d deeper search. Slower, a few more decodes\n");
+ printf(" -e x (x is transceiver dial frequency error in Hz)\n");
+ printf(" -f x (x is transceiver dial frequency in MHz)\n");
+ printf(" -H do not use (or update) the hash table\n");
+ printf(" -J use the stack decoder instead of Fano decoder\n");
+ printf(" -m decode wspr-15 .wav file\n");
+ printf(" -q quick mode - doesn't dig deep for weak signals\n");
+ printf(" -s single pass mode, no subtraction (same as original wsprd)\n");
+ printf(" -v verbose mode (shows dupes)\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.45)\n");
+}
+
+//***************************************************************************
+int main(int argc, char *argv[])
+{
+ char cr[] = "(C) 2016, 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 *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,more_candidates=0, stackdecoder=0;
+ int writenoise=0,usehashtable=1,wspr_type=2, ipass, nblocksize=1;
+ int writec2=0, npasses=2, subtraction=1;
+ int shift1, lagmin, lagmax, lagstep, ifmin, ifmax, worth_a_try, not_decoded;
+ unsigned int nbits=81, stacksize=200000;
+ 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;
+
+ struct result { char date[7]; char time[5]; float sync; float snr;
+ float dt; double freq; char message[23]; float drift;
+ unsigned int cycles; int jitter; };
+ struct result decodes[50];
+
+ char *hashtab;
+ hashtab=malloc(sizeof(char)*32768*13);
+ memset(hashtab,0,sizeof(char)*32768*13);
+ int nh;
+ symbols=malloc(sizeof(char)*nbits*2);
+ decdata=malloc(sizeof(char)*11);
+ channel_symbols=malloc(sizeof(char)*nbits*2);
+
+ callsign=malloc(sizeof(char)*13);
+ call_loc_pow=malloc(sizeof(char)*23);
+ float allfreqs[100];
+ char allcalls[100][13];
+ memset(allfreqs,0,sizeof(float)*100);
+ memset(allcalls,0,sizeof(char)*100*13);
+
+ int uniques=0, noprint=0, ndecodes_pass=0;
+
+ // 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
+ float minsync2=0.10; //Second sync limit
+ int iifac=8; //Step size in final DT peakup
+ int symfac=50; //Soft-symbol normalizing factor
+ int maxdrift=4; //Maximum (+/-) drift
+ 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];
+
+ idat=malloc(sizeof(float)*maxpts);
+ qdat=malloc(sizeof(float)*maxpts);
+
+ while ( (c = getopt(argc, argv, "a:cC:de:f:HJmqstwvz:")) !=-1 ) {
+ switch (c) {
+ case 'a':
+ data_dir = optarg;
+ break;
+ case 'c':
+ writec2=1;
+ break;
+ case 'C':
+ maxcycles=(unsigned int) strtoul(optarg,NULL,10);
+ break;
+ case 'd':
+ more_candidates=1;
+ break;
+ case 'e':
+ dialfreq_error = strtod(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 'J': //Stack (Jelinek) decoder, Fano decoder is the default
+ stackdecoder = 1;
+ break;
+ case 'm': //15-minute wspr mode
+ wspr_type = 15;
+ break;
+ case 'q': //no shift jittering
+ quickmode = 1;
+ break;
+ case 's': //single pass mode (same as original wsprd)
+ subtraction = 0;
+ npasses = 1;
+ break;
+ case 'v':
+ verbose = 1;
+ break;
+ case 'w':
+ fmin=-150.0;
+ fmax=150.0;
+ break;
+ case 'z':
+ bias=strtod(optarg,NULL); //fano metric bias (default is 0.45)
+ break;
+ case '?':
+ usage();
+ return 1;
+ }
+ }
+
+ if( stackdecoder ) {
+ stack=malloc(stacksize*sizeof(struct snode));
+ }
+
+ 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_fftwf_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_fftwf_wisdom_file = fopen(wisdom_fname, "r"))) { //Open FFTW wisdom
+ fftwf_import_wisdom_from_file(fp_fftwf_wisdom_file);
+ fclose(fp_fftwf_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,"%f %f %f %f %f %f %f",
+ &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 += (float)(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=(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);
+ }
+
+ 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);
+ }
+
+ //*************** main loop starts here *****************
+ for (ipass=0; ipass 0 && ndecodes_pass == 0 ) break;
+ ndecodes_pass=0;
+
+ memset(ps,0.0, sizeof(float)*512*nffts);
+ for (i=0; i511 )
+ k=k-512;
+ ps[j][i]=fftout[k][0]*fftout[k][0]+fftout[k][1]*fftout[k][1];
+ }
+ }
+
+ // Compute average spectrum
+ memset(psavg,0.0, sizeof(float)*512);
+ for (i=0; imin_snr) && (npk<200);
+ if ( candidate ) {
+ freq0[npk]=(j-205)*df;
+ snr0[npk]=10*log10(smspec[j])-snr_scaling_factor;
+ npk++;
+ }
+ }
+ } else {
+ 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;
+ 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= fmin && freq0[j] <= fmax ) {
+ freq0[i]=freq0[j];
+ snr0[i]=snr0[j];
+ i++;
+ }
+ }
+ npk=i;
+
+ // bubble sort on snr, bringing freq along for the ride
+ 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 smax ) { //Save coarse parameters
+ 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;
+
+ /*
+ 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; jminsync1 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;
+
+ 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);
+
+ // refine drift estimate
+ 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;
+
+ // fine-grid lag and freq search
+ if( sync1 > minsync1 ) {
+
+ 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;
+
+ // fine search over frequency
+ 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;
+
+ worth_a_try = 1;
+ } else {
+ worth_a_try = 0;
+ }
+
+ int idt=0, ii=0, jiggered_shift;
+ float 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, ifmin, ifmax, fstep,
+// &jiggered_shift, lagmin, lagmax, lagstep, &drift1, symfac,
+// &sync1, 2);
+ tsync2 += (float)(clock()-t0)/CLOCKS_PER_SEC;
+
+ noncoherent_sequence_detection(idat, qdat, npoints, symbols, &f1,
+ &jiggered_shift, &drift1, symfac,
+ &sync1, &nblocksize);
+
+ 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)) {
+ deinterleave(symbols);
+ t0 = clock();
+
+ if ( stackdecoder ) {
+ not_decoded = jelinek(&metric, &cycles, decdata, symbols, nbits,
+ stacksize, stack, mettab,maxcycles);
+ } else {
+ not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
+ mettab,delta,maxcycles);
+ }
+
+ tfano += (float)(clock()-t0)/CLOCKS_PER_SEC;
+
+ }
+ idt++;
+ if( quickmode ) break;
+ }
+
+ if( worth_a_try && !not_decoded ) {
+ ndecodes_pass++;
+
+ 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);
+
+ // subtract even on last pass
+ if( subtraction && (ipass < npasses ) && !noprint ) {
+ if( get_wspr_channel_symbols(call_loc_pow, hashtab, channel_symbols) ) {
+ subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
+ } else {
+ break;
+ }
+
+ }
+
+ // Remove dupes (same callsign and freq within 3 Hz)
+ int dupe=0;
+ for (i=0; i decodes[k+1].freq) {
+ temp = decodes[k];
+ decodes[k]=decodes[k+1];;
+ decodes[k+1] = temp;
+ }
+ }
+ }
+
+ for (i=0; i\n");
+
+ fftwf_free(fftin);
+ fftwf_free(fftout);
+
+ if ((fp_fftwf_wisdom_file = fopen(wisdom_fname, "w"))) {
+ fftwf_export_wisdom_to_file(fp_fftwf_wisdom_file);
+ fclose(fp_fftwf_wisdom_file);
+ }
+
+ 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(fall_wspr);
+ fclose(fwsprd);
+ // fclose(fdiag);
+ fclose(ftimer);
+ fftwf_destroy_plan(PLAN1);
+ fftwf_destroy_plan(PLAN2);
+ fftwf_destroy_plan(PLAN3);
+
+ if( usehashtable ) {
+ fhash=fopen(hash_fname,"w");
+ for (i=0; i<32768; i++) {
+ if( strncmp(hashtab+i*13,"\0",1) != 0 ) {
+ fprintf(fhash,"%5d %s\n",i,hashtab+i*13);
+ }
+ }
+ fclose(fhash);
+ }
+
+ if( stackdecoder ) {
+ free(stack);
+ }
+
+ if(writenoise == 999) return -1; //Silence compiler warning
+ return 0;
+}