diff --git a/CMakeLists.txt b/CMakeLists.txt index 59e46dfe1..3f2abbdd7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -470,6 +470,7 @@ set (wsprd_CSRCS lib/wsprd/wsprsim_utils.c lib/wsprd/wsprd_utils.c lib/wsprd/fano.c + lib/wsprd/jelinek.c lib/wsprd/tab.c lib/wsprd/nhash.c lib/init_random_seed.c @@ -775,7 +776,7 @@ endif () find_package (OpenMP) # -# fftw3 single precsion library +# fftw3 single precision library # find_package (FFTW3 COMPONENTS double single threads REQUIRED) diff --git a/lib/wsprd/wsprd.c b/lib/wsprd/wsprd.c index eb689de71..391661e08 100644 --- a/lib/wsprd/wsprd.c +++ b/lib/wsprd/wsprd.c @@ -37,16 +37,16 @@ #include #include "fano.h" +#include "jelinek.h" #include "nhash.h" #include "wsprd_utils.h" #include "wsprsim_utils.h" -//#include "lib/init_random_seed.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; +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, @@ -64,15 +64,18 @@ unsigned long nr; int printdata=0; //*************************************************************************** -unsigned long readc2file(char *ptr_to_infile, double *idat, double *qdat, +unsigned long readc2file(char *ptr_to_infile, float *idat, float *qdat, double *freq, int *wspr_type) { - float buffer[2*65536]; + 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); @@ -96,12 +99,13 @@ unsigned long readc2file(char *ptr_to_infile, double *idat, double *qdat, } else { return 1; } + free(buffer); } //*************************************************************************** -unsigned long readwavfile(char *ptr_to_infile, int ntrmin, double *idat, double *qdat ) +unsigned long readwavfile(char *ptr_to_infile, int ntrmin, float *idat, float *qdat ) { - int i, j, npoints; + unsigned long i, j, npoints; int nfft1, nfft2, nh2, i0; double df; @@ -123,8 +127,8 @@ unsigned long readwavfile(char *ptr_to_infile, int ntrmin, double *idat, double return 1; } - double *realin; - fftw_complex *fftin, *fftout; + float *realin; + fftwf_complex *fftin, *fftout; FILE *fp; short int *buf2; @@ -139,9 +143,9 @@ unsigned long readwavfile(char *ptr_to_infile, int ntrmin, double *idat, double 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); + 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; i0) & (k0) && (k syncmax ) { //Save best parameters - syncmax=ss/totp; + ss=ss/totp; + if( ss > syncmax ) { //Save best parameters + syncmax=ss; best_shift=lag; fbest=f0; } @@ -330,16 +333,16 @@ void sync_and_demodulate(double *id, double *qd, long np, /*************************************************************************** symbol-by-symbol signal subtraction ****************************************************************************/ -void subtract_signal(double *id, double *qd, long np, +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; - double pi=4.*atan(1.0),twopidt, fp; + float pi=4.*atan(1.0),twopidt, fp; - double i0,q0; - double c0[256],s0[256]; - double dphi, cdphi, sdphi; + float i0,q0; + float c0[256],s0[256]; + float dphi, cdphi, sdphi; twopidt=2*pi*dt; @@ -386,30 +389,30 @@ void subtract_signal(double *id, double *qd, long np, /****************************************************************************** Fully coherent signal subtraction *******************************************************************************/ -void subtract_signal2(double *id, double *qd, long np, +void subtract_signal2(float *id, float *qd, long np, float f0, int shift0, float drift0, unsigned char* channel_symbols) { - double dt=1.0/375.0, df=375.0/256.0; - double pi=4.*atan(1.0), twopidt, phi=0, dphi, cs; + float dt=1.0/375.0, df=375.0/256.0; + float pi=4.*atan(1.0), twopidt, phi=0, dphi, cs; int i, j, k, ii, nsym=162, nspersym=256, nfilt=256; //nfilt must be even number. int nsig=nsym*nspersym; int nc2=45000; - double *refi, *refq, *ci, *cq, *cfi, *cfq; + float *refi, *refq, *ci, *cq, *cfi, *cfq; - refi=malloc(sizeof(double)*nc2); - refq=malloc(sizeof(double)*nc2); - ci=malloc(sizeof(double)*nc2); - cq=malloc(sizeof(double)*nc2); - cfi=malloc(sizeof(double)*nc2); - cfq=malloc(sizeof(double)*nc2); + refi=malloc(sizeof(float)*nc2); + refq=malloc(sizeof(float)*nc2); + ci=malloc(sizeof(float)*nc2); + cq=malloc(sizeof(float)*nc2); + cfi=malloc(sizeof(float)*nc2); + cfq=malloc(sizeof(float)*nc2); - memset(refi,0,sizeof(double)*nc2); - memset(refq,0,sizeof(double)*nc2); - memset(ci,0,sizeof(double)*nc2); - memset(cq,0,sizeof(double)*nc2); - memset(cfi,0,sizeof(double)*nc2); - memset(cfq,0,sizeof(double)*nc2); + memset(refi,0,sizeof(float)*nc2); + memset(refq,0,sizeof(float)*nc2); + memset(ci,0,sizeof(float)*nc2); + memset(cq,0,sizeof(float)*nc2); + memset(cfi,0,sizeof(float)*nc2); + memset(cfq,0,sizeof(float)*nc2); twopidt=2.0*pi*dt; @@ -425,18 +428,18 @@ void subtract_signal2(double *id, double *qd, long np, // for (i=0; i0) & (k0) && (k0) & (k0) && (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.42)\n"); + printf(" -z x (x is fano metric table bias, default is 0.45)\n"); } //*************************************************************************** @@ -567,7 +573,7 @@ int main(int argc, char *argv[]) extern char *optarg; extern int optind; int i,j,k; - unsigned char *symbols, *decdata; + 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; @@ -575,38 +581,39 @@ int main(int argc, char *argv[]) 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; + int c,delta,maxpts=65536,verbose=0,quickmode=0,more_candidates=0, stackdecoder=0; int writenoise=0,usehashtable=1,wspr_type=2, ipass; int writec2=0, npasses=2, subtraction=1; - int shift1, lagmin, lagmax, lagstep, worth_a_try, not_decoded; - unsigned int nbits=81; - unsigned int npoints, metric, maxcycles, cycles, maxnp; + 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; - float dialfreq_error=0.0; + double dialfreq_error=0.0; float fmin=-110, fmax=110; float f1, fstep, sync1, drift1; float psavg[512]; - double *idat, *qdat; + float *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; + 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; + float dt; float 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((nbits+7)/8); + 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]; @@ -614,31 +621,29 @@ int main(int argc, char *argv[]) memset(allfreqs,0,sizeof(float)*100); memset(allcalls,0,sizeof(char)*100*13); - int uniques=0, noprint=0; - -// init_random_seed(); - + int uniques=0, noprint=0, ndecodes_pass=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 + 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 maxdrift=4; //Maximum (+/-) drift - double minrms=52.0 * (symfac/64.0); //Final test for plausible decoding + 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(); - fftw_complex *fftin, *fftout; + fftwf_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); + idat=malloc(sizeof(float)*maxpts); + qdat=malloc(sizeof(float)*maxpts); - while ( (c = getopt(argc, argv, "a:ce:f:Hmqstwvz:")) !=-1 ) { + while ( (c = getopt(argc, argv, "a:cC:de:f:HJmqstwvz:")) !=-1 ) { switch (c) { case 'a': data_dir = optarg; @@ -646,8 +651,14 @@ int main(int argc, char *argv[]) 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 = strtof(optarg,NULL); // units of Hz + dialfreq_error = strtod(optarg,NULL); // units of Hz // dialfreq_error = dial reading - actual, correct frequency break; case 'f': @@ -656,14 +667,17 @@ int main(int argc, char *argv[]) case 'H': usehashtable = 0; break; - case 'm': + 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': + case 'q': //no shift jittering quickmode = 1; break; - case 's': - subtraction = 0; //single pass mode (same as original wsprd) + case 's': //single pass mode (same as original wsprd) + subtraction = 0; npasses = 1; break; case 'v': @@ -674,7 +688,7 @@ int main(int argc, char *argv[]) fmax=150.0; break; case 'z': - bias=strtof(optarg,NULL); //fano metric bias (default is 0.42) + bias=strtod(optarg,NULL); //fano metric bias (default is 0.45) break; case '?': usage(); @@ -682,6 +696,10 @@ int main(int argc, char *argv[]) } } + if( stackdecoder ) { + stack=malloc(stacksize*sizeof(struct snode)); + } + if( optind+1 > argc) { usage(); return 1; @@ -695,7 +713,7 @@ int main(int argc, char *argv[]) mettab[1][i]=round( 10*(metric_tables[2][255-i]-bias) ); } - FILE *fp_fftw_wisdom_file, *fall_wspr, *fwsprd, *fhash, *ftimer; + FILE *fp_fftwf_wisdom_file, *fall_wspr, *fwsprd, *fhash, *ftimer; strcpy(wisdom_fname,"."); strcpy(all_fname,"."); strcpy(spots_fname,"."); @@ -713,9 +731,9 @@ int main(int argc, char *argv[]) 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); + 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"); @@ -725,7 +743,7 @@ int main(int argc, char *argv[]) if((ftimer=fopen(timer_fname,"r"))) { //Accumulate timing data - nr=fscanf(ftimer,"%lf %lf %lf %lf %lf %lf %lf", + nr=fscanf(ftimer,"%f %f %f %f %f %f %f", &treadwav,&tcandidates,&tsync0,&tsync1,&tsync2,&tfano,&ttotal); fclose(ftimer); } @@ -736,7 +754,7 @@ int main(int argc, char *argv[]) t0 = clock(); npoints=readwavfile(ptr_to_infile, wspr_type, idat, qdat); - treadwav += (double)(clock()-t0)/CLOCKS_PER_SEC; + treadwav += (float)(clock()-t0)/CLOCKS_PER_SEC; if( npoints == 1 ) { return 1; @@ -760,12 +778,12 @@ int main(int argc, char *argv[]) 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); + 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]; @@ -785,14 +803,12 @@ int main(int argc, char *argv[]) } 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 ) @@ -853,7 +869,7 @@ int main(int argc, char *argv[]) } for (j=0; j<411; j++) { smspec[j]=smspec[j]/noise_level - 1.0; - if( smspec[j] < min_snr) smspec[j]=0.1; + if( smspec[j] < min_snr) smspec[j]=0.1*min_snr; continue; } @@ -867,14 +883,29 @@ int main(int argc, char *argv[]) } 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++; + unsigned char candidate; + if( more_candidates ) { + for(j=0; j<411; j=j+2) { + candidate = (smspec[j]>min_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; @@ -907,6 +938,7 @@ int main(int argc, char *argv[]) } t0=clock(); + /* Make coarse estimates of shift (DT), freq, and drift * Look for time offsets up to +/- 8 symbols (about +/- 5.4 s) relative @@ -926,14 +958,14 @@ int main(int argc, char *argv[]) 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); @@ -968,8 +1000,8 @@ int main(int argc, char *argv[]) } } } - tcandidates += (double)(clock()-t0)/CLOCKS_PER_SEC; - + 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 @@ -986,43 +1018,76 @@ int main(int argc, char *argv[]) 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, fstep, &shift1, + sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, 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; + tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC; + + fstep=0.25; ifmin=-2; ifmax=2; t0 = clock(); - sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep, &shift1, + sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1, lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1); - tsync1 += (double)(clock()-t0)/CLOCKS_PER_SEC; + + // 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; - double y,sq,rms; + float y,sq,rms; not_decoded=1; while ( worth_a_try && not_decoded && idt<=(128/iifac)) { @@ -1033,38 +1098,39 @@ int main(int argc, char *argv[]) // Use mode 2 to get soft-decision symbols t0 = clock(); - sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep, + sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &jiggered_shift, lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 2); - tsync2 += (double)(clock()-t0)/CLOCKS_PER_SEC; - + tsync2 += (float)(clock()-t0)/CLOCKS_PER_SEC; + sq=0.0; for(i=0; i<162; i++) { - y=(double)symbols[i] - 128.0; + y=(float)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; + 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; - /* ### 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 ) { + ndecodes_pass++; for(i=0; i<11; i++) { @@ -1075,16 +1141,14 @@ int main(int argc, char *argv[]) } } - + // 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); - - if( subtraction && (ipass == 0) && !noprint ) { - - unsigned char channel_symbols[162]; - + + // 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 { @@ -1092,7 +1156,7 @@ int main(int argc, char *argv[]) } } - + // Remove dupes (same callsign and freq within 3 Hz) int dupe=0; for (i=0; i\n"); - fftw_free(fftin); - fftw_free(fftout); + fftwf_free(fftin); + fftwf_free(fftout); - if ((fp_fftw_wisdom_file = fopen(wisdom_fname, "w"))) { - fftw_export_wisdom_to_file(fp_fftw_wisdom_file); - fclose(fp_fftw_wisdom_file); + if ((fp_fftwf_wisdom_file = fopen(wisdom_fname, "w"))) { + fftwf_export_wisdom_to_file(fp_fftwf_wisdom_file); + fclose(fp_fftwf_wisdom_file); } - - ttotal += (double)(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); @@ -1201,7 +1256,7 @@ int main(int argc, char *argv[]) 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,"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); @@ -1209,9 +1264,9 @@ int main(int argc, char *argv[]) fclose(fwsprd); // fclose(fdiag); fclose(ftimer); - fftw_destroy_plan(PLAN1); - fftw_destroy_plan(PLAN2); - fftw_destroy_plan(PLAN3); + fftwf_destroy_plan(PLAN1); + fftwf_destroy_plan(PLAN2); + fftwf_destroy_plan(PLAN3); if( usehashtable ) { fhash=fopen(hash_fname,"w"); @@ -1223,6 +1278,10 @@ int main(int argc, char *argv[]) fclose(fhash); } + if( stackdecoder ) { + free(stack); + } + if(writenoise == 999) return -1; //Silence compiler warning return 0; } diff --git a/lib/wsprd/wsprd_exp.c b/lib/wsprd/wsprd_exp.c deleted file mode 100644 index 391661e08..000000000 --- a/lib/wsprd/wsprd_exp.c +++ /dev/null @@ -1,1287 +0,0 @@ -/* - 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); - - *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 ) -{ - 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; - } - - 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; inh2 ) 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; i0) && (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; - } - return; - } - 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[]) -{ - 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; - 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; float 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 - 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; - - 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; -}