/*
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 "nhash.h"
#include "wsprd_utils.h"
#include "wsprsim_utils.h"
#define max(x,y) ((x) > (y) ? (x) : (y))
// Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT,
// FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
#define PATIENCE FFTW_ESTIMATE
fftw_plan PLAN1,PLAN2,PLAN3;
unsigned char pr3[162]=
{1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0,
0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1,
0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1,
1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1,
0,0,1,0,1,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0,
0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1,
0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1,
0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,0,1,1,0,
0,0};
unsigned long nr;
int printdata=0;
//***************************************************************************
unsigned long readc2file(char *ptr_to_infile, double *idat, double *qdat,
double *freq, int *wspr_type)
{
float buffer[2*65536];
double dfreq;
int i,ntrmin;
char *c2file[15];
FILE* fp;
fp = fopen(ptr_to_infile,"rb");
if (fp == NULL) {
fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile);
return 1;
}
unsigned long nread=fread(c2file,sizeof(char),14,fp);
nread=fread(&ntrmin,sizeof(int),1,fp);
nread=fread(&dfreq,sizeof(double),1,fp);
*freq=dfreq;
nread=fread(buffer,sizeof(float),2*45000,fp);
*wspr_type=ntrmin;
for(i=0; i<45000; i++) {
idat[i]=buffer[2*i];
qdat[i]=-buffer[2*i+1];
}
if( nread == 2*45000 ) {
return nread/2;
} else {
return 1;
}
}
//***************************************************************************
unsigned long readwavfile(char *ptr_to_infile, int ntrmin, double *idat, double *qdat )
{
unsigned long i, j, npoints;
int nfft1, nfft2, nh2, i0;
double df;
nfft2=46080; //this is the number of downsampled points that will be returned
nh2=nfft2/2;
if( ntrmin == 2 ) {
nfft1=nfft2*32; //need to downsample by a factor of 32
df=12000.0/nfft1;
i0=1500.0/df+0.5;
npoints=114*12000;
} else if ( ntrmin == 15 ) {
nfft1=nfft2*8*32;
df=12000.0/nfft1;
i0=(1500.0+112.5)/df+0.5;
npoints=8*114*12000;
} else {
fprintf(stderr,"This should not happen\n");
return 1;
}
double *realin;
fftw_complex *fftin, *fftout;
FILE *fp;
short int *buf2;
buf2 = malloc(npoints*sizeof(short int));
fp = fopen(ptr_to_infile,"rb");
if (fp == NULL) {
fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile);
return 1;
}
nr=fread(buf2,2,22,fp); //Read and ignore header
nr=fread(buf2,2,npoints,fp); //Read raw data
fclose(fp);
realin=(double*) fftw_malloc(sizeof(double)*nfft1);
fftout=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*nfft1);
PLAN1 = fftw_plan_dft_r2c_1d(nfft1, realin, fftout, PATIENCE);
for (i=0; inh2 ) j=j-nfft2;
fftin[i][0]=fftout[j][0];
fftin[i][1]=fftout[j][1];
}
fftw_free(fftout);
fftout=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*nfft2);
PLAN2 = fftw_plan_dft_1d(nfft2, fftin, fftout, FFTW_BACKWARD, PATIENCE);
fftw_execute(PLAN2);
for (i=0; i0) & (k syncmax ) { //Save best parameters
syncmax=ss/totp;
best_shift=lag;
fbest=f0;
}
} // lag loop
} //freq loop
if( mode <=1 ) { //Send best params back to caller
*sync=syncmax;
*shift1=best_shift;
*f1=fbest;
return;
}
if( mode == 2 ) {
*sync=syncmax;
for (i=0; i<162; i++) { //Normalize the soft symbols
fsum=fsum+fsymb[i]/162.0;
f2sum=f2sum+fsymb[i]*fsymb[i]/162.0;
}
fac=sqrt(f2sum-fsum*fsum);
for (i=0; i<162; i++) {
fsymb[i]=symfac*fsymb[i]/fac;
if( fsymb[i] > 127) fsymb[i]=127.0;
if( fsymb[i] < -128 ) fsymb[i]=-128.0;
symbols[i]=fsymb[i] + 128;
}
return;
}
return;
}
/***************************************************************************
symbol-by-symbol signal subtraction
****************************************************************************/
void subtract_signal(double *id, double *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols)
{
float dt=1.0/375.0, df=375.0/256.0;
int i, j, k;
double pi=4.*atan(1.0),twopidt, fp;
double i0,q0;
double c0[256],s0[256];
double 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(" -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(" -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");
}
//***************************************************************************
int main(int argc, char *argv[])
{
extern char *optarg;
extern int optind;
int i,j,k;
unsigned char *symbols, *decdata;
signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0};
char *callsign, *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;
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;
float df=375.0/256.0/2;
float freq0[200],snr0[200],drift0[200],sync0[200];
int shift0[200];
float dt=1.0/375.0, dt_print;
double dialfreq_cmdline=0.0, dialfreq, freq_print;
float dialfreq_error=0.0;
float fmin=-110, fmax=110;
float f1, fstep, sync1, drift1;
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;
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((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
double minsync2=0.12; //Second sync limit
int iifac=3; //Step size in final DT peakup
int symfac=50; //Soft-symbol normalizing factor
int maxdrift=4; //Maximum (+/-) drift
double minrms=52.0 * (symfac/64.0); //Final test for plausible decoding
delta=60; //Fano threshold step
t00=clock();
fftw_complex *fftin, *fftout;
#include "./metric_tables.c"
int mettab[2][256];
float bias=0.42;
idat=malloc(sizeof(double)*maxpts);
qdat=malloc(sizeof(double)*maxpts);
while ( (c = getopt(argc, argv, "a:ce:f:Hmqstwvz:")) !=-1 ) {
switch (c) {
case 'a':
data_dir = optarg;
break;
case 'c':
writec2=1;
break;
case 'e':
dialfreq_error = strtof(optarg,NULL); // units of Hz
// dialfreq_error = dial reading - actual, correct frequency
break;
case 'f':
dialfreq_cmdline = strtod(optarg,NULL); // units of MHz
break;
case 'H':
usehashtable = 0;
break;
case 'm':
wspr_type = 15;
break;
case 'q':
quickmode = 1;
break;
case 's':
subtraction = 0; //single pass mode (same as original wsprd)
npasses = 1;
break;
case 'v':
verbose = 1;
break;
case 'w':
fmin=-150.0;
fmax=150.0;
break;
case 'z':
bias=strtof(optarg,NULL); //fano metric bias (default is 0.42)
break;
case '?':
usage();
return 1;
}
}
if( optind+1 > argc) {
usage();
return 1;
} else {
ptr_to_infile=argv[optind];
}
// setup metric table
for(i=0; i<256; i++) {
mettab[0][i]=round( 10*(metric_tables[2][i]-bias) );
mettab[1][i]=round( 10*(metric_tables[2][255-i]-bias) );
}
FILE *fp_fftw_wisdom_file, *fall_wspr, *fwsprd, *fhash, *ftimer;
strcpy(wisdom_fname,".");
strcpy(all_fname,".");
strcpy(spots_fname,".");
strcpy(timer_fname,".");
strcpy(hash_fname,".");
if(data_dir != NULL) {
strcpy(wisdom_fname,data_dir);
strcpy(all_fname,data_dir);
strcpy(spots_fname,data_dir);
strcpy(timer_fname,data_dir);
strcpy(hash_fname,data_dir);
}
strncat(wisdom_fname,"/wspr_wisdom.dat",20);
strncat(all_fname,"/ALL_WSPR.TXT",20);
strncat(spots_fname,"/wspr_spots.txt",20);
strncat(timer_fname,"/wspr_timer.out",20);
strncat(hash_fname,"/hashtable.txt",20);
if ((fp_fftw_wisdom_file = fopen(wisdom_fname, "r"))) { //Open FFTW wisdom
fftw_import_wisdom_from_file(fp_fftw_wisdom_file);
fclose(fp_fftw_wisdom_file);
}
fall_wspr=fopen(all_fname,"a");
fwsprd=fopen(spots_fname,"w");
// FILE *fdiag;
// fdiag=fopen("wsprd_diag","a");
if((ftimer=fopen(timer_fname,"r"))) {
//Accumulate timing data
nr=fscanf(ftimer,"%lf %lf %lf %lf %lf %lf %lf",
&treadwav,&tcandidates,&tsync0,&tsync1,&tsync2,&tfano,&ttotal);
fclose(ftimer);
}
ftimer=fopen(timer_fname,"w");
if( strstr(ptr_to_infile,".wav") ) {
ptr_to_infile_suffix=strstr(ptr_to_infile,".wav");
t0 = clock();
npoints=readwavfile(ptr_to_infile, wspr_type, idat, qdat);
treadwav += (double)(clock()-t0)/CLOCKS_PER_SEC;
if( npoints == 1 ) {
return 1;
}
dialfreq=dialfreq_cmdline - (dialfreq_error*1.0e-06);
} else if ( strstr(ptr_to_infile,".c2") !=0 ) {
ptr_to_infile_suffix=strstr(ptr_to_infile,".c2");
npoints=readc2file(ptr_to_infile, idat, qdat, &dialfreq, &wspr_type);
if( npoints == 1 ) {
return 1;
}
dialfreq -= (dialfreq_error*1.0e-06);
} else {
printf("Error: Failed to open %s\n",ptr_to_infile);
printf("WSPR file must have suffix .wav or .c2\n");
return 1;
}
// Parse date and time from given filename
strncpy(date,ptr_to_infile_suffix-11,6);
strncpy(uttime,ptr_to_infile_suffix-4,4);
date[6]='\0';
uttime[4]='\0';
// Do windowed ffts over 2 symbols, stepped by half symbols
int nffts=4*floor(npoints/512)-1;
fftin=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*512);
fftout=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*512);
PLAN3 = fftw_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE);
float ps[512][nffts];
float w[512];
for(i=0; i<512; i++) {
w[i]=sin(0.006147931*i);
}
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; ipass511 )
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; ismspec[j-1]) && (smspec[j]>smspec[j+1]) && (npk<200)) {
freq0[npk]=(j-205)*df;
snr0[npk]=10*log10(smspec[j])-snr_scaling_factor;
npk++;
}
}
// Compute corrected fmin, fmax, accounting for dial frequency error
fmin += dialfreq_error; // dialfreq_error is in units of Hz
fmax += dialfreq_error;
// Don't waste time on signals outside of the range [fmin,fmax].
i=0;
for( j=0; j= 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 += (double)(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; j minsync1 ) {
worth_a_try = 1;
} else {
worth_a_try = 0;
}
int idt=0, ii=0, jiggered_shift;
double y,sq,rms;
not_decoded=1;
while ( worth_a_try && not_decoded && idt<=(128/iifac)) {
ii=(idt+1)/2;
if( idt%2 == 1 ) ii=-ii;
ii=iifac*ii;
jiggered_shift=shift1+ii;
// Use mode 2 to get soft-decision symbols
t0 = clock();
sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep,
&jiggered_shift, lagmin, lagmax, lagstep, &drift1, symfac,
&sync1, 2);
tsync2 += (double)(clock()-t0)/CLOCKS_PER_SEC;
sq=0.0;
for(i=0; i<162; i++) {
y=(double)symbols[i] - 128.0;
sq += y*y;
}
rms=sqrt(sq/162.0);
if((sync1 > minsync2) && (rms > minrms)) {
deinterleave(symbols);
t0 = clock();
not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
mettab,delta,maxcycles);
tfano += (double)(clock()-t0)/CLOCKS_PER_SEC;
/* ### Used for timing tests:
if(not_decoded) fprintf(fdiag,
"%6s %4s %4.1f %3.0f %4.1f %10.7f %-18s %2d %5u %4d %6.1f %2d\n",
date,uttime,sync1*10,snr0[j], shift1*dt-2.0, dialfreq+(1500+f1)/1e6,
"@ ", (int)drift1, cycles/81, ii, rms, maxnp);
*/
}
idt++;
if( quickmode ) break;
}
if( worth_a_try && !not_decoded ) {
for(i=0; i<11; i++) {
if( decdata[i]>127 ) {
message[i]=decdata[i]-256;
} else {
message[i]=decdata[i];
}
}
// Unpack the decoded message, update the hashtable, apply
// sanity checks on grid and power, and return
// call_loc_pow string and also callsign (for de-duping).
noprint=unpk_(message,hashtab,call_loc_pow,callsign);
if( subtraction && (ipass == 0) && !noprint ) {
unsigned char channel_symbols[162];
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");
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);
}
ttotal += (double)(clock()-t00)/CLOCKS_PER_SEC;
fprintf(ftimer,"%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n\n",
treadwav,tcandidates,tsync0,tsync1,tsync2,tfano,ttotal);
fprintf(ftimer,"Code segment Seconds Frac\n");
fprintf(ftimer,"-----------------------------------\n");
fprintf(ftimer,"readwavfile %7.2f %7.2f\n",treadwav,treadwav/ttotal);
fprintf(ftimer,"Coarse DT f0 f1 %7.2f %7.2f\n",tcandidates,
tcandidates/ttotal);
fprintf(ftimer,"sync_and_demod(0) %7.2f %7.2f\n",tsync0,tsync0/ttotal);
fprintf(ftimer,"sync_and_demod(1) %7.2f %7.2f\n",tsync1,tsync1/ttotal);
fprintf(ftimer,"sync_and_demod(2) %7.2f %7.2f\n",tsync2,tsync2/ttotal);
fprintf(ftimer,"Fano decoder %7.2f %7.2f\n",tfano,tfano/ttotal);
fprintf(ftimer,"-----------------------------------\n");
fprintf(ftimer,"Total %7.2f %7.2f\n",ttotal,1.0);
fclose(fall_wspr);
fclose(fwsprd);
// fclose(fdiag);
fclose(ftimer);
fftw_destroy_plan(PLAN1);
fftw_destroy_plan(PLAN2);
fftw_destroy_plan(PLAN3);
if( usehashtable ) {
fhash=fopen(hash_fname,"w");
for (i=0; i<32768; i++) {
if( strncmp(hashtab+i*13,"\0",1) != 0 ) {
fprintf(fhash,"%5d %s\n",i,hashtab+i*13);
}
}
fclose(fhash);
}
if(writenoise == 999) return -1; //Silence compiler warning
return 0;
}