Further work on signal subtraction. wsprd_exp now has subtraction and two-pass options.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@5617 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/wsprd/wsprd_exp.c

@@ -194,9 +194,9 @@ void sync_and_demodulate(double *id, double *qd, long np,
     float dt=1.0/375.0, df=375.0/256.0, fbest=0.0;
     int i, j, k;
     double pi=4.*atan(1.0),twopidt;
-    float f0=0.0,fp,fplast=-10000.0,ss;
+    float f0=0.0,fp,ss;
     int lag;
-    
+    static float fplast=-10000.0;
     double i0[162],q0[162],i1[162],q1[162],i2[162],q2[162],i3[162],q3[162];
     double p0,p1,p2,p3,cmet,totp,syncmax,fac;
     double c0[256],s0[256],c1[256],s1[256],c2[256],s2[256],c3[256],s3[256];
@@ -220,7 +220,7 @@ void sync_and_demodulate(double *id, double *qd, long np,
             for (i=0; i<162; i++) {
                 fp = f0 + ((float)*drift1/2.0)*((float)i-81.0)/81.0;
                 if( i==0 || (fp != fplast) ) {  // only calculate sin/cos if necessary
-                    dphi0=2*pi*(fp-1.5*df)*dt;
+                    dphi0=twopidt*(fp-1.5*df);
                     cdphi0=cos(dphi0);
                     sdphi0=sin(dphi0);
                     
@@ -366,27 +366,17 @@ void subtract_signal(double *id, double *qd, long np,
         
         // subtract the signal here.
         
-        i0=i0/256.0;
+        i0=i0/256.0; //will be wrong for partial symbols at the edges...
         q0=q0/256.0;
-        double p0=i0*i0+q0*q0;
-        double is=0, qs=0;
         
         for (j=0; j<256; j++) {
             k=shift0+i*256+j;
             if( (k>0) & (k<np) ) {
                 id[k]=id[k]- (i0*c0[j] - q0*s0[j]);
                 qd[k]=qd[k]- (q0*c0[j] + i0*s0[j]);
-                is=is+id[k];
-                qs=qs+qd[k];
             }
         }
-        is=is/256.0;
-        qs=qs/256.0;
-        double p1=is*is+qs*qs;
-//        printf("symbol %d i0: %f q0: %f is: %f qs: %f\n",i,i0,q0,is,qs);
-        
     }
-    
     return;
 }
 
@@ -438,7 +428,8 @@ void usage(void)
     printf("       -m decode wspr-15 .wav file\n");
     printf("       -n write noise estimates to file noise.dat\n");
     printf("       -q quick mode - doesn't dig deep for weak signals\n");
-    printf("       -s slow mode - much slower, yields a few more decodes\n");
+    printf("       -s signal subtraction mode\n");
+    printf("       -t signal subtraction followed by a second pass\n");
     printf("       -v verbose mode\n");
     printf("       -w wideband mode - decode signals within +/- 150 Hz of center\n");
     printf("       -z x (x is fano metric table bias, default is 0.42)\n");
@@ -452,15 +443,16 @@ int main(int argc, char *argv[])
     int i,j,k;
     unsigned char *symbols, *decdata;
     signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0};
-    char *callsign,*grid,*grid6, *call_loc_pow, *cdbm;
+    char *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,writenoise=0,usehashtable=1,wspr_type=2;
+    int c,delta,maxpts=65536,verbose=0,quickmode=0;
+    int writenoise=0,usehashtable=1,wspr_type=2, subtraction=0, ipass, npasses=1;
     int shift1, lagmin, lagmax, lagstep, worth_a_try, not_decoded;
-    unsigned int nbits;
+    unsigned int nbits=81;
     unsigned int npoints, metric, maxcycles, cycles, maxnp;
     float df=375.0/256.0/2;
     float freq0[200],snr0[200],drift0[200],sync0[200];
@@ -469,12 +461,27 @@ int main(int argc, char *argv[])
     double dialfreq_cmdline=0.0, dialfreq, freq_print;
     float dialfreq_error=0.0;
     float fmin=-110, fmax=110;
-    float f1, fstep, sync1, drift1, tblank=0, fblank=0;
+    float f1, fstep, sync1, drift1;
+    float psavg[512];
     double *idat, *qdat;
     clock_t t0,t00;
     double tfano=0.0,treadwav=0.0,tcandidates=0.0,tsync0=0.0;
     double tsync1=0.0,tsync2=0.0,ttotal=0.0;
     
+    char hashtab[32768][13];
+    memset(hashtab,0,sizeof(char)*32768*13);
+    int nh;
+    symbols=malloc(sizeof(char)*nbits*2);
+    decdata=malloc((nbits+7)/8);
+    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;
+    
     // Parameters used for performance-tuning:
     maxcycles=10000;                         //Fano timeout limit
     double minsync1=0.10;                    //First sync limit
@@ -495,14 +502,11 @@ int main(int argc, char *argv[])
     idat=malloc(sizeof(double)*maxpts);
     qdat=malloc(sizeof(double)*maxpts);
     
-    while ( (c = getopt(argc, argv, "a:b:e:f:Hmnqst:wvz:")) !=-1 ) {
+    while ( (c = getopt(argc, argv, "a:b:e:f:Hmnqstwvz:")) !=-1 ) {
         switch (c) {
             case 'a':
                 data_dir = optarg;
                 break;
-            case 'b':
-                fblank = strtof(optarg,NULL);
-                break;
             case 'e':
                 dialfreq_error = strtof(optarg,NULL);   // units of Hz
                 // dialfreq_error = dial reading - actual, correct frequency
@@ -523,11 +527,11 @@ int main(int argc, char *argv[])
                 quickmode = 1;
                 break;
             case 's':
-                maxcycles=20000;
+                subtraction = 1;
-                iifac=1;
                 break;
             case 't':
-                tblank = strtof(optarg,NULL);
+                subtraction = 1;
+                npasses = 2; //npasses defaults to 1
                 break;
             case 'v':
                 verbose = 1;
@@ -636,6 +640,22 @@ int main(int argc, char *argv[])
         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<npasses; ipass++) {
+        
         memset(ps,0.0, sizeof(float)*512*nffts);
         for (i=0; i<nffts; i++) {
             for(j=0; j<512; j++ ) {
@@ -652,11 +672,7 @@ int main(int argc, char *argv[])
             }
         }
         
-    fftw_free(fftin);
-    fftw_free(fftout);
-    
         // Compute average spectrum
-    float psavg[512];
         memset(psavg,0.0, sizeof(float)*512);
         for (i=0; i<nffts; i++) {
             for (j=0; j<512; j++) {
@@ -736,6 +752,22 @@ int main(int argc, char *argv[])
         }
         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
          
@@ -800,42 +832,6 @@ int main(int argc, char *argv[])
         }
         tcandidates += (double)(clock()-t0)/CLOCKS_PER_SEC;
         
-    nbits=81;
-    symbols=malloc(sizeof(char)*nbits*2);
-    memset(symbols,0,sizeof(char)*nbits*2);
-    decdata=malloc((nbits+7)/8);
-    grid=malloc(sizeof(char)*5);
-    grid6=malloc(sizeof(char)*7);
-    callsign=malloc(sizeof(char)*13);
-    call_loc_pow=malloc(sizeof(char)*23);
-    cdbm=malloc(sizeof(char)*3);
-    float allfreqs[npk];
-    memset(allfreqs,0,sizeof(float)*npk);
-    char allcalls[npk][13];
-    memset(allcalls,0,sizeof(char)*npk*13);
-    memset(grid,0,sizeof(char)*5);
-    memset(grid6,0,sizeof(char)*7);
-    memset(callsign,0,sizeof(char)*13);
-    memset(call_loc_pow,0,sizeof(char)*23);
-    memset(cdbm,0,sizeof(char)*3);
-    char hashtab[32768][13];
-    memset(hashtab,0,sizeof(char)*32768*13);
-    int nh;
-    
-    if( usehashtable ) {
-        char line[80], hcall[12];
-        if( (fhash=fopen(hash_fname,"r+")) ) {
-            while (fgets(line, sizeof(line), fhash) != NULL) {
-                sscanf(line,"%d %s",&nh,hcall);
-                strcpy(*hashtab+nh*13,hcall);
-            }
-        } else {
-            fhash=fopen(hash_fname,"w+");
-        }
-        fclose(fhash);
-    }
-    
-    int uniques=0, noprint=0;
         /*
          Refine the estimates of freq, shift using sync as a metric.
          Sync is calculated such that it is a float taking values in the range
@@ -854,6 +850,11 @@ int main(int argc, char *argv[])
          */
         
         for (j=0; j<npk; j++) {
+            
+            memset(symbols,0,sizeof(char)*nbits*2);
+            memset(callsign,0,sizeof(char)*13);
+            memset(call_loc_pow,0,sizeof(char)*23);
+
             f1=freq0[j];
             drift1=drift0[j];
             shift1=shift0[j];
@@ -926,12 +927,15 @@ int main(int argc, char *argv[])
             }
             
             if( worth_a_try && !not_decoded ) {
+ 
                 for(i=0; i<11; i++) {
+               
                     if( decdata[i]>127 ) {
                         message[i]=decdata[i]-256;
                     } else {
                         message[i]=decdata[i];
                     }
+                    
                 }
                 
                 // Unpack the decoded message, update the hashtable, apply
@@ -939,16 +943,23 @@ int main(int argc, char *argv[])
                 // call_loc_pow string and also callsign (for de-duping).
                 noprint=unpk_(message,hashtab,call_loc_pow,callsign);
                 
+                if( subtraction ) {
+
                     unsigned char channel_symbols[162];
-            get_wspr_channel_symbols(call_loc_pow, channel_symbols);
     
+                    if( get_wspr_channel_symbols(call_loc_pow, channel_symbols) ) {
                         subtract_signal(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
+                    } else {
+                        break;
+                    }
                     
-            // Remove dupes (same callsign and freq within 1 Hz)
+                }
+                
+                // Remove dupes (same callsign and freq within 3 Hz)
                 int dupe=0;
-            for (i=0; i<npk; i++) {
+                for (i=0; i<uniques; i++) {
                     if(!strcmp(callsign,allcalls[i]) &&
-                   (fabs(f1-allfreqs[i]) <1.0)) dupe=1;
+                       (fabs(f1-allfreqs[i]) <3.0)) dupe=1;
                 }
                 if( (verbose || !dupe) && !noprint) {
                     uniques++;
@@ -990,8 +1001,12 @@ int main(int argc, char *argv[])
                 }
             }
         }
+    }
     printf("<DecodeFinished>\n");
     
+    fftw_free(fftin);
+    fftw_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);
@@ -1032,13 +1047,13 @@ int main(int argc, char *argv[])
         fclose(fhash);
     }
     
-    char c2filename[15];
+//    char c2filename[15];
-    double carrierfreq=dialfreq;
+//    double carrierfreq=dialfreq;
-    int wsprtype=2;
+//    int wsprtype=2;
-    strcpy(c2filename,"000000_0001.c2");
+//    strcpy(c2filename,"000000_0001.c2");
-    printf("Writing %s\n",c2filename);
+//    printf("Writing %s\n",c2filename);
-    writec2file(c2filename, wsprtype, carrierfreq, idat, qdat);
+//    writec2file(c2filename, wsprtype, carrierfreq, idat, qdat);
     
-    if(fblank+tblank+writenoise == 999) return -1;  //Silence compiler warning
+    if(writenoise == 999) return -1;  //Silence compiler warning
     return 0;
 }

lib/wsprd/wsprsim.c

@@ -124,22 +124,27 @@ int main(int argc, char *argv[])
 {
     extern char *optarg;
     extern int optind;
-    int i, c, printchannel=0;
+    int i, c, printchannel=0, writec2=0;
     float snr=50.0;
     char *message, *c2filename;
+    c2filename=malloc(sizeof(char)*15);
+
     // message length is 22 characters
     message=malloc(sizeof(char)*23);
-    c2filename=malloc(sizeof(char)*15);
-    memset(c2filename,0,sizeof(char)*15);
     
+    strcpy(c2filename,"000000_0001.c2");
+    printf("%s\n",c2filename);
     while ( (c = getopt(argc, argv, "cdo:s:")) !=-1 ) {
         switch (c) {
             case 'c':
                 printchannel=1;
+                break;
             case 'd':
                 printdata=1;
+                break;
             case 'o':
                 c2filename = optarg;
+                writec2=1;
                 break;
             case 's':
                 snr = (float)atoi(optarg);
@@ -187,14 +192,10 @@ int main(int argc, char *argv[])
     f0=0.0;
     t0=1.0;
     add_signal_vector(f0, t0, snr, channel_symbols, isig, qsig);
-    
+    if( writec2) {
-    if( strlen(c2filename) >0 ) {
         // write a .c2 file
         double carrierfreq=10.1387;
         int wsprtype=2;
-        if( strlen(c2filename) != 14 ) {
-            strcpy(c2filename,"000000_0001.c2");
-        }
         printf("Writing %s\n",c2filename);
         writec2file(c2filename, wsprtype, carrierfreq, isig, qsig);
     }

lib/wsprd/wsprsim_utils.c

@@ -236,7 +236,6 @@ int get_wspr_channel_symbols(char* rawmessage, unsigned char* symbols) {
         }
         grid6[5]=grid[0];
         n=pack_call(grid6);
-        printf("Callsign %s hash %d\n",callsign,ihash);
     } else if ( i2 < mlen ) {  // just looks for a right slash
         // Type 2: PJ4/K1ABC 37
         callsign=strtok(message," ");
@@ -251,7 +250,7 @@ int get_wspr_channel_symbols(char* rawmessage, unsigned char* symbols) {
         m=128*ng+ntype+64;
         n=n1;
     } else {
-        printf("Error: bad message format\n");
+//        printf("Error: bad message format\n");
         return 0;
     }