From e2c76b59915def1cab1931a47258e56193d5fa7f Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Tue, 2 Jan 2018 02:16:44 +0000 Subject: [PATCH] Experimental noncoherent block sequence detector for wspr. Not yet completely functional. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8393 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/wsprd/wsprd_exp.c | 1431 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1431 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..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; +}