1 /*
2 This file is part of program wsprd, a detector/demodulator/decoder
3 for the Weak Signal Propagation Reporter (WSPR) mode.
4
5 File name: wsprd.c
6
7 Copyright 2001-2018, Joe Taylor, K1JT
8
9 Much of the present code is based on work by Steven Franke, K9AN,
10 which in turn was based on earlier work by K1JT.
11
12 Copyright 2014-2018, Steven Franke, K9AN
13
14 License: GNU GPL v3
15
16 This program is free software: you can redistribute it and/or modify
17 it under the terms of the GNU General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 This program is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU General Public License for more details.
25
26 You should have received a copy of the GNU General Public License
27 along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29
30 #include <stdio.h>
31 #include <unistd.h>
32 #include <stdlib.h>
33 #include <math.h>
34 #include <string.h>
35 #include <stdint.h>
36 #include <time.h>
37 #include <fftw3.h>
38
39 #include "fano.h"
40 #include "jelinek.h"
41 #include "nhash.h"
42 #include "wsprd_utils.h"
43 #include "wsprsim_utils.h"
44
45 #define max(x,y) ((x) > (y) ? (x) : (y))
46
47 extern void osdwspr_ (float [], unsigned char [], int *, unsigned char [], int *, float *);
48
49 // Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT,
50 // FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
51 #define PATIENCE FFTW_ESTIMATE
52 fftwf_plan PLAN1,PLAN2,PLAN3;
53
54 unsigned char pr3[162]=
55 {1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0,
56 0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1,
57 0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1,
58 1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1,
59 0,0,1,0,1,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0,
60 0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1,
61 0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1,
62 0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,0,1,1,0,
63 0,0};
64
65 int printdata=0;
66
67 //***************************************************************************
readc2file(char * ptr_to_infile,float * idat,float * qdat,double * freq,int * wspr_type)68 unsigned long readc2file(char *ptr_to_infile, float *idat, float *qdat,
69 double *freq, int *wspr_type)
70 {
71 float *buffer;
72 double dfreq;
73 int i,ntrmin;
74 char c2file[15];
75 size_t nr;
76 FILE* fp;
77
78 fp = fopen(ptr_to_infile,"rb");
79 if (fp == NULL) {
80 fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile);
81 return 1;
82 }
83 nr=fread(c2file,sizeof(char),14,fp);
84 nr=fread(&ntrmin,sizeof(int),1,fp);
85 nr=fread(&dfreq,sizeof(double),1,fp);
86 *freq=dfreq;
87
88 buffer=calloc(2*65536,sizeof(float));
89 nr=fread(buffer,sizeof(float),2*45000,fp);
90 fclose(fp);
91
92 *wspr_type=ntrmin;
93
94 for(i=0; i<45000; i++) {
95 idat[i]=buffer[2*i];
96 qdat[i]=-buffer[2*i+1];
97 }
98 free(buffer);
99
100 if( nr == 2*45000 ) {
101 return (unsigned long) nr/2;
102 } else {
103 return 1;
104 }
105 }
106
107 //***************************************************************************
readwavfile(char * ptr_to_infile,int ntrmin,float * idat,float * qdat)108 unsigned long readwavfile(char *ptr_to_infile, int ntrmin, float *idat, float *qdat )
109 {
110 size_t i, j, npoints, nr;
111 int nfft1, nfft2, nh2, i0;
112 double df;
113
114 nfft2=46080; //this is the number of downsampled points that will be returned
115 nh2=nfft2/2;
116
117 if( ntrmin == 2 ) {
118 nfft1=nfft2*32; //need to downsample by a factor of 32
119 df=12000.0/nfft1;
120 i0=1500.0/df+0.5;
121 npoints=114*12000;
122 } else if ( ntrmin == 15 ) {
123 nfft1=nfft2*8*32;
124 df=12000.0/nfft1;
125 i0=(1500.0+112.5)/df+0.5;
126 npoints=8*114*12000;
127 } else {
128 fprintf(stderr,"This should not happen\n");
129 return 1;
130 }
131
132 float *realin;
133 fftwf_complex *fftin, *fftout;
134
135 FILE *fp;
136 short int *buf2;
137
138 fp = fopen(ptr_to_infile,"rb");
139 if (fp == NULL) {
140 fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile);
141 return 1;
142 }
143
144 buf2 = calloc(npoints,sizeof(short int));
145 nr=fread(buf2,2,22,fp); //Read and ignore header
146 nr=fread(buf2,2,npoints,fp); //Read raw data
147 fclose(fp);
148 if( nr == 0 ) {
149 fprintf(stderr, "No data in file '%s'\n", ptr_to_infile);
150 return 1;
151 }
152
153 realin=(float*) fftwf_malloc(sizeof(float)*nfft1);
154 fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*(nfft1/2+1));
155 PLAN1 = fftwf_plan_dft_r2c_1d(nfft1, realin, fftout, PATIENCE);
156
157 for (i=0; i<npoints; i++) {
158 realin[i]=buf2[i]/32768.0;
159 }
160
161 for (i=npoints; i<(size_t)nfft1; i++) {
162 realin[i]=0.0;
163 }
164 free(buf2);
165
166 fftwf_execute(PLAN1);
167 fftwf_free(realin);
168
169 fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*nfft2);
170
171 for (i=0; i<(size_t)nfft2; i++) {
172 j=i0+i;
173 if( i>(size_t)nh2 ) j=j-nfft2;
174 fftin[i][0]=fftout[j][0];
175 fftin[i][1]=fftout[j][1];
176 }
177
178 fftwf_free(fftout);
179 fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*nfft2);
180 PLAN2 = fftwf_plan_dft_1d(nfft2, fftin, fftout, FFTW_BACKWARD, PATIENCE);
181 fftwf_execute(PLAN2);
182
183 for (i=0; i<(size_t)nfft2; i++) {
184 idat[i]=fftout[i][0]/1000.0;
185 qdat[i]=fftout[i][1]/1000.0;
186 }
187
188 fftwf_free(fftin);
189 fftwf_free(fftout);
190 return nfft2;
191 }
192
193 //***************************************************************************
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)194 void sync_and_demodulate(float *id, float *qd, long np,
195 unsigned char *symbols, float *f1, int ifmin, int ifmax, float fstep,
196 int *shift1, int lagmin, int lagmax, int lagstep,
197 float *drift1, int symfac, float *sync, int mode)
198 {
199 /***********************************************************************
200 * mode = 0: no frequency or drift search. find best time lag. *
201 * 1: no time lag or drift search. find best frequency. *
202 * 2: no frequency or time lag search. calculate soft-decision *
203 * symbols using passed frequency and shift. *
204 ************************************************************************/
205
206 static float fplast=-10000.0;
207 static float dt=1.0/375.0, df=375.0/256.0;
208 static float pi=3.14159265358979323846;
209 float twopidt, df15=df*1.5, df05=df*0.5;
210
211 int i, j, k, lag;
212 float i0[162],q0[162],i1[162],q1[162],i2[162],q2[162],i3[162],q3[162];
213 float p0,p1,p2,p3,cmet,totp,syncmax,fac;
214 float c0[256],s0[256],c1[256],s1[256],c2[256],s2[256],c3[256],s3[256];
215 float dphi0, cdphi0, sdphi0, dphi1, cdphi1, sdphi1, dphi2, cdphi2, sdphi2,
216 dphi3, cdphi3, sdphi3;
217 float f0=0.0, fp, ss, fbest=0.0, fsum=0.0, f2sum=0.0, fsymb[162];
218 int best_shift = 0, ifreq;
219
220 syncmax=-1e30;
221 if( mode == 0 ) {ifmin=0; ifmax=0; fstep=0.0; f0=*f1;}
222 if( mode == 1 ) {lagmin=*shift1;lagmax=*shift1;f0=*f1;}
223 if( mode == 2 ) {lagmin=*shift1;lagmax=*shift1;ifmin=0;ifmax=0;f0=*f1;}
224
225 twopidt=2*pi*dt;
226 for(ifreq=ifmin; ifreq<=ifmax; ifreq++) {
227 f0=*f1+ifreq*fstep;
228 for(lag=lagmin; lag<=lagmax; lag=lag+lagstep) {
229 ss=0.0;
230 totp=0.0;
231 for (i=0; i<162; i++) {
232 fp = f0 + (*drift1/2.0)*((float)i-81.0)/81.0;
233 if( i==0 || (fp != fplast) ) { // only calculate sin/cos if necessary
234 dphi0=twopidt*(fp-df15);
235 cdphi0=cos(dphi0);
236 sdphi0=sin(dphi0);
237
238 dphi1=twopidt*(fp-df05);
239 cdphi1=cos(dphi1);
240 sdphi1=sin(dphi1);
241
242 dphi2=twopidt*(fp+df05);
243 cdphi2=cos(dphi2);
244 sdphi2=sin(dphi2);
245
246 dphi3=twopidt*(fp+df15);
247 cdphi3=cos(dphi3);
248 sdphi3=sin(dphi3);
249
250 c0[0]=1; s0[0]=0;
251 c1[0]=1; s1[0]=0;
252 c2[0]=1; s2[0]=0;
253 c3[0]=1; s3[0]=0;
254
255 for (j=1; j<256; j++) {
256 c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0;
257 s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0;
258 c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1;
259 s1[j]=c1[j-1]*sdphi1 + s1[j-1]*cdphi1;
260 c2[j]=c2[j-1]*cdphi2 - s2[j-1]*sdphi2;
261 s2[j]=c2[j-1]*sdphi2 + s2[j-1]*cdphi2;
262 c3[j]=c3[j-1]*cdphi3 - s3[j-1]*sdphi3;
263 s3[j]=c3[j-1]*sdphi3 + s3[j-1]*cdphi3;
264 }
265 fplast = fp;
266 }
267
268 i0[i]=0.0; q0[i]=0.0;
269 i1[i]=0.0; q1[i]=0.0;
270 i2[i]=0.0; q2[i]=0.0;
271 i3[i]=0.0; q3[i]=0.0;
272
273 for (j=0; j<256; j++) {
274 k=lag+i*256+j;
275 if( (k>0) && (k<np) ) {
276 i0[i]=i0[i] + id[k]*c0[j] + qd[k]*s0[j];
277 q0[i]=q0[i] - id[k]*s0[j] + qd[k]*c0[j];
278 i1[i]=i1[i] + id[k]*c1[j] + qd[k]*s1[j];
279 q1[i]=q1[i] - id[k]*s1[j] + qd[k]*c1[j];
280 i2[i]=i2[i] + id[k]*c2[j] + qd[k]*s2[j];
281 q2[i]=q2[i] - id[k]*s2[j] + qd[k]*c2[j];
282 i3[i]=i3[i] + id[k]*c3[j] + qd[k]*s3[j];
283 q3[i]=q3[i] - id[k]*s3[j] + qd[k]*c3[j];
284 }
285 }
286 p0=i0[i]*i0[i] + q0[i]*q0[i];
287 p1=i1[i]*i1[i] + q1[i]*q1[i];
288 p2=i2[i]*i2[i] + q2[i]*q2[i];
289 p3=i3[i]*i3[i] + q3[i]*q3[i];
290
291 p0=sqrt(p0);
292 p1=sqrt(p1);
293 p2=sqrt(p2);
294 p3=sqrt(p3);
295
296 totp=totp+p0+p1+p2+p3;
297 cmet=(p1+p3)-(p0+p2);
298 ss = (pr3[i] == 1) ? ss+cmet : ss-cmet;
299 if( mode == 2) { //Compute soft symbols
300 if(pr3[i]==1) {
301 fsymb[i]=p3-p1;
302 } else {
303 fsymb[i]=p2-p0;
304 }
305 }
306 }
307 ss=ss/totp;
308 if( ss > syncmax ) { //Save best parameters
309 syncmax=ss;
310 best_shift=lag;
311 fbest=f0;
312 }
313 } // lag loop
314 } //freq loop
315
316 if( mode <=1 ) { //Send best params back to caller
317 *sync=syncmax;
318 *shift1=best_shift;
319 *f1=fbest;
320 return;
321 }
322
323 if( mode == 2 ) {
324 *sync=syncmax;
325 for (i=0; i<162; i++) { //Normalize the soft symbols
326 fsum=fsum+fsymb[i]/162.0;
327 f2sum=f2sum+fsymb[i]*fsymb[i]/162.0;
328 }
329 fac=sqrt(f2sum-fsum*fsum);
330 for (i=0; i<162; i++) {
331 fsymb[i]=symfac*fsymb[i]/fac;
332 if( fsymb[i] > 127) fsymb[i]=127.0;
333 if( fsymb[i] < -128 ) fsymb[i]=-128.0;
334 symbols[i]=fsymb[i] + 128;
335 }
336 return;
337 }
338 return;
339 }
340
noncoherent_sequence_detection(float * id,float * qd,long np,unsigned char * symbols,float * f1,int * shift1,float * drift1,int symfac,int * nblocksize,int * bitmetric)341 void noncoherent_sequence_detection(float *id, float *qd, long np,
342 unsigned char *symbols, float *f1, int *shift1,
343 float *drift1, int symfac, int *nblocksize, int *bitmetric)
344 {
345 /************************************************************************
346 * Noncoherent sequence detection for wspr. *
347 * Allowed block lengths are nblock=1,2,3,6, or 9 symbols. *
348 * Longer block lengths require longer channel coherence time. *
349 * The whole block is estimated at once. *
350 * nblock=1 corresponds to noncoherent detection of individual symbols *
351 * like the original wsprd symbol demodulator. *
352 ************************************************************************/
353 static float fplast=-10000.0;
354 static float dt=1.0/375.0, df=375.0/256.0;
355 static float pi=3.14159265358979323846;
356 float twopidt, df15=df*1.5, df05=df*0.5;
357
358 int i, j, k, lag, itone, ib, b, nblock, nseq, imask;
359 float xi[512],xq[512];
360 float is[4][162],qs[4][162],cf[4][162],sf[4][162],cm,sm,cmp,smp;
361 float p[512],fac,xm1,xm0;
362 float c0[257],s0[257],c1[257],s1[257],c2[257],s2[257],c3[257],s3[257];
363 float dphi0, cdphi0, sdphi0, dphi1, cdphi1, sdphi1, dphi2, cdphi2, sdphi2,
364 dphi3, cdphi3, sdphi3;
365 float f0, fp, fsum=0.0, f2sum=0.0, fsymb[162];
366
367 twopidt=2*pi*dt;
368 f0=*f1;
369 lag=*shift1;
370 nblock=*nblocksize;
371 nseq=1<<nblock;
372 int bitbybit=*bitmetric;
373
374 for (i=0; i<162; i++) {
375 fp = f0 + (*drift1/2.0)*((float)i-81.0)/81.0;
376 if( i==0 || (fp != fplast) ) { // only calculate sin/cos if necessary
377 dphi0=twopidt*(fp-df15);
378 cdphi0=cos(dphi0);
379 sdphi0=sin(dphi0);
380
381 dphi1=twopidt*(fp-df05);
382 cdphi1=cos(dphi1);
383 sdphi1=sin(dphi1);
384
385 dphi2=twopidt*(fp+df05);
386 cdphi2=cos(dphi2);
387 sdphi2=sin(dphi2);
388
389 dphi3=twopidt*(fp+df15);
390 cdphi3=cos(dphi3);
391 sdphi3=sin(dphi3);
392
393 c0[0]=1; s0[0]=0;
394 c1[0]=1; s1[0]=0;
395 c2[0]=1; s2[0]=0;
396 c3[0]=1; s3[0]=0;
397
398 for (j=1; j<257; j++) {
399 c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0;
400 s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0;
401 c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1;
402 s1[j]=c1[j-1]*sdphi1 + s1[j-1]*cdphi1;
403 c2[j]=c2[j-1]*cdphi2 - s2[j-1]*sdphi2;
404 s2[j]=c2[j-1]*sdphi2 + s2[j-1]*cdphi2;
405 c3[j]=c3[j-1]*cdphi3 - s3[j-1]*sdphi3;
406 s3[j]=c3[j-1]*sdphi3 + s3[j-1]*cdphi3;
407 }
408
409 fplast = fp;
410 }
411
412 cf[0][i]=c0[256]; sf[0][i]=s0[256];
413 cf[1][i]=c1[256]; sf[1][i]=s1[256];
414 cf[2][i]=c2[256]; sf[2][i]=s2[256];
415 cf[3][i]=c3[256]; sf[3][i]=s3[256];
416
417 is[0][i]=0.0; qs[0][i]=0.0;
418 is[1][i]=0.0; qs[1][i]=0.0;
419 is[2][i]=0.0; qs[2][i]=0.0;
420 is[3][i]=0.0; qs[3][i]=0.0;
421
422 for (j=0; j<256; j++) {
423 k=lag+i*256+j;
424 if( (k>0) && (k<np) ) {
425 is[0][i]=is[0][i] + id[k]*c0[j] + qd[k]*s0[j];
426 qs[0][i]=qs[0][i] - id[k]*s0[j] + qd[k]*c0[j];
427 is[1][i]=is[1][i] + id[k]*c1[j] + qd[k]*s1[j];
428 qs[1][i]=qs[1][i] - id[k]*s1[j] + qd[k]*c1[j];
429 is[2][i]=is[2][i] + id[k]*c2[j] + qd[k]*s2[j];
430 qs[2][i]=qs[2][i] - id[k]*s2[j] + qd[k]*c2[j];
431 is[3][i]=is[3][i] + id[k]*c3[j] + qd[k]*s3[j];
432 qs[3][i]=qs[3][i] - id[k]*s3[j] + qd[k]*c3[j];
433 }
434 }
435 }
436
437 for (i=0; i<162; i=i+nblock) {
438 for (j=0;j<nseq;j++) {
439 xi[j]=0.0; xq[j]=0.0;
440 cm=1; sm=0;
441 for (ib=0; ib<nblock; ib++) {
442 b=(j&(1<<(nblock-1-ib)))>>(nblock-1-ib);
443 itone=pr3[i+ib]+2*b;
444 xi[j]=xi[j]+is[itone][i+ib]*cm + qs[itone][i+ib]*sm;
445 xq[j]=xq[j]+qs[itone][i+ib]*cm - is[itone][i+ib]*sm;
446 cmp=cf[itone][i+ib]*cm - sf[itone][i+ib]*sm;
447 smp=sf[itone][i+ib]*cm + cf[itone][i+ib]*sm;
448 cm=cmp; sm=smp;
449 }
450 p[j]=xi[j]*xi[j]+xq[j]*xq[j];
451 p[j]=sqrt(p[j]);
452 }
453 for (ib=0; ib<nblock; ib++) {
454 imask=1<<(nblock-1-ib);
455 xm1=0.0; xm0=0.0;
456 for (j=0; j<nseq; j++) {
457 if((j & imask)!=0) {
458 if(p[j] > xm1) xm1=p[j];
459 }
460 if((j & imask)==0) {
461 if(p[j]>xm0) xm0=p[j];
462 }
463 }
464 fsymb[i+ib]=xm1-xm0;
465 if( bitbybit == 1 ) {
466 fsymb[i+ib]=fsymb[i+ib]/(xm1 > xm0 ? xm1 : xm0);
467 }
468 }
469 }
470 for (i=0; i<162; i++) { //Normalize the soft symbols
471 fsum=fsum+fsymb[i]/162.0;
472 f2sum=f2sum+fsymb[i]*fsymb[i]/162.0;
473 }
474 fac=sqrt(f2sum-fsum*fsum);
475 for (i=0; i<162; i++) {
476 fsymb[i]=symfac*fsymb[i]/fac;
477 if( fsymb[i] > 127) fsymb[i]=127.0;
478 if( fsymb[i] < -128 ) fsymb[i]=-128.0;
479 symbols[i]=fsymb[i] + 128;
480 }
481 return;
482 }
483
484 /***************************************************************************
485 symbol-by-symbol signal subtraction
486 ****************************************************************************/
subtract_signal(float * id,float * qd,long np,float f0,int shift0,float drift0,unsigned char * channel_symbols)487 void subtract_signal(float *id, float *qd, long np,
488 float f0, int shift0, float drift0, unsigned char* channel_symbols)
489 {
490 float dt=1.0/375.0, df=375.0/256.0;
491 int i, j, k;
492 float pi=4.*atan(1.0),twopidt, fp;
493
494 float i0,q0;
495 float c0[256],s0[256];
496 float dphi, cdphi, sdphi;
497
498 twopidt=2*pi*dt;
499
500 for (i=0; i<162; i++) {
501 fp = f0 + ((float)drift0/2.0)*((float)i-81.0)/81.0;
502
503 dphi=twopidt*(fp+((float)channel_symbols[i]-1.5)*df);
504 cdphi=cos(dphi);
505 sdphi=sin(dphi);
506
507 c0[0]=1; s0[0]=0;
508
509 for (j=1; j<256; j++) {
510 c0[j]=c0[j-1]*cdphi - s0[j-1]*sdphi;
511 s0[j]=c0[j-1]*sdphi + s0[j-1]*cdphi;
512 }
513
514 i0=0.0; q0=0.0;
515
516 for (j=0; j<256; j++) {
517 k=shift0+i*256+j;
518 if( (k>0) & (k<np) ) {
519 i0=i0 + id[k]*c0[j] + qd[k]*s0[j];
520 q0=q0 - id[k]*s0[j] + qd[k]*c0[j];
521 }
522 }
523
524
525 // subtract the signal here.
526
527 i0=i0/256.0; //will be wrong for partial symbols at the edges...
528 q0=q0/256.0;
529
530 for (j=0; j<256; j++) {
531 k=shift0+i*256+j;
532 if( (k>0) & (k<np) ) {
533 id[k]=id[k]- (i0*c0[j] - q0*s0[j]);
534 qd[k]=qd[k]- (q0*c0[j] + i0*s0[j]);
535 }
536 }
537 }
538 return;
539 }
540 /******************************************************************************
541 Subtract the coherent component of a signal
542 *******************************************************************************/
subtract_signal2(float * id,float * qd,long np,float f0,int shift0,float drift0,unsigned char * channel_symbols)543 void subtract_signal2(float *id, float *qd, long np,
544 float f0, int shift0, float drift0, unsigned char* channel_symbols)
545 {
546 float dt=1.0/375.0, df=375.0/256.0;
547 float pi=4.*atan(1.0), twopidt, phi=0, dphi, cs;
548 int i, j, k, ii, nsym=162, nspersym=256, nfilt=360; //nfilt must be even number.
549 int nsig=nsym*nspersym;
550 int nc2=45000;
551
552 float *refi, *refq, *ci, *cq, *cfi, *cfq;
553
554 refi=calloc(nc2,sizeof(float));
555 refq=calloc(nc2,sizeof(float));
556 ci=calloc(nc2,sizeof(float));
557 cq=calloc(nc2,sizeof(float));
558 cfi=calloc(nc2,sizeof(float));
559 cfq=calloc(nc2,sizeof(float));
560
561 twopidt=2.0*pi*dt;
562
563 /******************************************************************************
564 Measured signal: s(t)=a(t)*exp( j*theta(t) )
565 Reference is: r(t) = exp( j*phi(t) )
566 Complex amplitude is estimated as: c(t)=LPF[s(t)*conjugate(r(t))]
567 so c(t) has phase angle theta-phi
568 Multiply r(t) by c(t) and subtract from s(t), i.e. s'(t)=s(t)-c(t)r(t)
569 *******************************************************************************/
570
571 // create reference wspr signal vector, centered on f0.
572 //
573 for (i=0; i<nsym; i++) {
574
575 cs=(float)channel_symbols[i];
576
577 dphi=twopidt*
578 (
579 f0 + (drift0/2.0)*((float)i-(float)nsym/2.0)/((float)nsym/2.0)
580 + (cs-1.5)*df
581 );
582
583 for ( j=0; j<nspersym; j++ ) {
584 ii=nspersym*i+j;
585 refi[ii]=cos(phi); //cannot precompute sin/cos because dphi is changing
586 refq[ii]=sin(phi);
587 phi=phi+dphi;
588 }
589 }
590
591 float w[nfilt], norm=0, partialsum[nfilt];
592 //lowpass filter and remove startup transient
593 for (i=0; i<nfilt; i++) partialsum[i]=0.0;
594 for (i=0; i<nfilt; i++) {
595 w[i]=sin(pi*(float)i/(float)(nfilt-1));
596 norm=norm+w[i];
597 }
598 for (i=0; i<nfilt; i++) {
599 w[i]=w[i]/norm;
600 }
601 for (i=1; i<nfilt; i++) {
602 partialsum[i]=partialsum[i-1]+w[i];
603 }
604
605 // s(t) * conjugate(r(t))
606 // beginning of first symbol in reference signal is at i=0
607 // beginning of first symbol in received data is at shift0.
608 // filter transient lasts nfilt samples
609 // leave nfilt zeros as a pad at the beginning of the unfiltered reference signal
610 for (i=0; i<nsym*nspersym; i++) {
611 k=shift0+i;
612 if( (k>0) && (k<np) ) {
613 ci[i+nfilt] = id[k]*refi[i] + qd[k]*refq[i];
614 cq[i+nfilt] = qd[k]*refi[i] - id[k]*refq[i];
615 }
616 }
617
618 // LPF
619 for (i=nfilt/2; i<45000-nfilt/2; i++) {
620 cfi[i]=0.0; cfq[i]=0.0;
621 for (j=0; j<nfilt; j++) {
622 cfi[i]=cfi[i]+w[j]*ci[i-nfilt/2+j];
623 cfq[i]=cfq[i]+w[j]*cq[i-nfilt/2+j];
624 }
625 }
626
627 // subtract c(t)*r(t) here
628 // (ci+j*cq)(refi+j*refq)=(ci*refi-cq*refq)+j(ci*refq+cq*refi)
629 // beginning of first symbol in reference signal is at i=nfilt
630 // beginning of first symbol in received data is at shift0.
631 for (i=0; i<nsig; i++) {
632 if( i<nfilt/2 ) { // take care of the end effect (LPF step response) here
633 norm=partialsum[nfilt/2+i];
634 } else if( i>(nsig-1-nfilt/2) ) {
635 norm=partialsum[nfilt/2+nsig-1-i];
636 } else {
637 norm=1.0;
638 }
639 k=shift0+i;
640 j=i+nfilt;
641 if( (k>0) && (k<np) ) {
642 id[k]=id[k] - (cfi[j]*refi[i]-cfq[j]*refq[i])/norm;
643 qd[k]=qd[k] - (cfi[j]*refq[i]+cfq[j]*refi[i])/norm;
644 }
645 }
646
647 free(refi);
648 free(refq);
649 free(ci);
650 free(cq);
651 free(cfi);
652 free(cfq);
653
654 return;
655 }
656
writec2file(char * c2filename,int trmin,double freq,float * idat,float * qdat)657 unsigned long writec2file(char *c2filename, int trmin, double freq
658 , float *idat, float *qdat)
659 {
660 int i;
661 float *buffer;
662 FILE *fp;
663
664 fp = fopen(c2filename,"wb");
665 if( fp == NULL ) {
666 fprintf(stderr, "Could not open c2 file '%s'\n", c2filename);
667 return 0;
668 }
669 unsigned long nwrite = fwrite(c2filename,sizeof(char),14,fp);
670 nwrite = fwrite(&trmin, sizeof(int), 1, fp);
671 nwrite = fwrite(&freq, sizeof(double), 1, fp);
672
673 buffer=calloc(2*45000,sizeof(float));
674 for(i=0; i<45000; i++) {
675 buffer[2*i]=idat[i];
676 buffer[2*i+1]=-qdat[i];
677 }
678 nwrite = fwrite(buffer, sizeof(float), 2*45000, fp);
679 free(buffer);
680 fclose(fp);
681 if( nwrite == 2*45000 ) {
682 return nwrite;
683 } else {
684 return 0;
685 }
686 }
687
count_hard_errors(unsigned char * symbols,unsigned char * channel_symbols)688 unsigned int count_hard_errors( unsigned char *symbols, unsigned char *channel_symbols)
689 {
690 int i,is;
691 unsigned char cw[162];
692 unsigned int nerrors;
693 for (i=0; i<162; i++) {
694 cw[i] = channel_symbols[i] >=2 ? 1:0;
695 }
696 deinterleave(cw);
697 nerrors=0;
698 for (i=0; i<162; i++) {
699 is = symbols[i] > 127 ? 1:0;
700 nerrors = nerrors + (is == cw[i] ? 0:1);
701 }
702 return nerrors;
703 }
704
705 //***************************************************************************
usage(void)706 void usage(void)
707 {
708 printf("Usage: wsprd [options...] infile\n");
709 printf(" infile must have suffix .wav or .c2\n");
710 printf("\n");
711 printf("Options:\n");
712 printf(" -a <path> path to writeable data files, default=\".\"\n");
713 printf(" -B disable block demodulation - use single-symbol noncoherent demod\n");
714 printf(" -c write .c2 file at the end of the first pass\n");
715 printf(" -C maximum number of decoder cycles per bit, default 10000\n");
716 printf(" -d deeper search. Slower, a few more decodes\n");
717 printf(" -e x (x is transceiver dial frequency error in Hz)\n");
718 printf(" -f x (x is transceiver dial frequency in MHz)\n");
719 printf(" -H do not use (or update) the hash table\n");
720 printf(" -J use the stack decoder instead of Fano decoder\n");
721 printf(" -m decode wspr-15 .wav file\n");
722 printf(" -o n (0<=n<=5), decoding depth for OSD, default is disabled\n");
723 printf(" -q quick mode - doesn't dig deep for weak signals\n");
724 printf(" -s single pass mode, no subtraction (same as original wsprd)\n");
725 printf(" -v verbose mode (shows dupes)\n");
726 printf(" -w wideband mode - decode signals within +/- 150 Hz of center\n");
727 printf(" -z x (x is fano metric table bias, default is 0.45)\n");
728 }
729
730 //***************************************************************************
main(int argc,char * argv[])731 int main(int argc, char *argv[])
732 {
733 char cr[] = "(C) 2018, Steven Franke - K9AN";
734 (void)cr;
735 extern char *optarg;
736 extern int optind;
737 int i,j,k;
738 unsigned char *symbols, *decdata, *channel_symbols, *apmask, *cw;
739 signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0};
740 char *callsign, *grid, *call_loc_pow;
741 char *ptr_to_infile,*ptr_to_infile_suffix;
742 char *data_dir=".";
743 char wisdom_fname[200],all_fname[200],spots_fname[200];
744 char timer_fname[200],hash_fname[200];
745 char uttime[5],date[7];
746 int c,delta,maxpts=65536,verbose=0,quickmode=0,more_candidates=0, stackdecoder=0;
747 int usehashtable=1,wspr_type=2, ipass, nblocksize;
748 int nhardmin,ihash;
749 int writec2=0,maxdrift;
750 int shift1, lagmin, lagmax, lagstep, ifmin, ifmax, not_decoded;
751 unsigned int nbits=81, stacksize=200000;
752 struct snode * stack=NULL;
753 unsigned int npoints, cycles, maxnp, metric;
754 float df=375.0/256.0/2;
755 float fsymbs[162];
756 float dt=1.0/375.0, dt_print;
757 double dialfreq_cmdline=0.0, dialfreq, freq_print;
758 double dialfreq_error=0.0;
759 float fmin=-110, fmax=110;
760 float f1, fstep, sync1, drift1;
761 float dmin;
762 float psavg[512];
763 float *idat, *qdat;
764 clock_t t0,t00;
765 float tfano=0.0,treadwav=0.0,tcandidates=0.0,tsync0=0.0;
766 float tsync1=0.0,tsync2=0.0,tosd=0.0,ttotal=0.0;
767
768 struct cand { float freq; float snr; int shift; float drift; float sync; };
769 struct cand candidates[200];
770
771 struct result { char date[7]; char time[5]; float sync; float snr;
772 float dt; double freq; char message[23]; float drift;
773 unsigned int cycles; int jitter; int blocksize; unsigned int metric;
774 int nhardmin; int ipass; int decodetype;};
775 struct result decodes[50];
776
777 char *hashtab;
778 hashtab=calloc(32768*13,sizeof(char));
779 char *loctab;
780 loctab=calloc(32768*5,sizeof(char));
781 int nh;
782 symbols=calloc(nbits*2,sizeof(unsigned char));
783 apmask=calloc(162,sizeof(unsigned char));
784 cw=calloc(162,sizeof(unsigned char));
785 decdata=calloc(11,sizeof(unsigned char));
786 channel_symbols=calloc(nbits*2,sizeof(unsigned char));
787 callsign=calloc(13,sizeof(char));
788 grid=calloc(5,sizeof(char));
789 call_loc_pow=calloc(23,sizeof(char));
790 float allfreqs[100];
791 char allcalls[100][13];
792 for (i=0; i<100; i++) allfreqs[i]=0.0;
793 memset(allcalls,0,sizeof(char)*100*13);
794
795 int uniques=0, noprint=0, ndecodes_pass=0;
796
797 // Parameters used for performance-tuning:
798 unsigned int maxcycles=10000; //Decoder timeout limit
799 float minsync1=0.10; //First sync limit
800 float minsync2=0.12; //Second sync limit
801 int iifac=8; //Step size in final DT peakup
802 int symfac=50; //Soft-symbol normalizing factor
803 int subtraction=1;
804 int npasses=3;
805 int ndepth=-1; //Depth for OSD
806
807 float minrms=52.0 * (symfac/64.0); //Final test for plausible decoding
808 delta=60; //Fano threshold step
809 float bias=0.45; //Fano metric bias (used for both Fano and stack algorithms)
810
811 t00=clock();
812 fftwf_complex *fftin, *fftout;
813 #include "./metric_tables.c"
814
815 int mettab[2][256];
816
817 idat=calloc(maxpts,sizeof(float));
818 qdat=calloc(maxpts,sizeof(float));
819
820 while ( (c = getopt(argc, argv, "a:BcC:de:f:HJmo:qstwvz:")) !=-1 ) {
821 switch (c) {
822 case 'a':
823 data_dir = optarg;
824 break;
825 case 'B':
826 npasses=2;
827 break;
828 case 'c':
829 writec2=1;
830 break;
831 case 'C':
832 maxcycles=(unsigned int) strtoul(optarg,NULL,10);
833 break;
834 case 'd':
835 more_candidates=1;
836 break;
837 case 'e':
838 dialfreq_error = strtod(optarg,NULL); // units of Hz
839 // dialfreq_error = dial reading - actual, correct frequency
840 break;
841 case 'f':
842 dialfreq_cmdline = strtod(optarg,NULL); // units of MHz
843 break;
844 case 'H':
845 usehashtable = 0;
846 break;
847 case 'J': //Stack (Jelinek) decoder, Fano decoder is the default
848 stackdecoder = 1;
849 break;
850 case 'm': //15-minute wspr mode
851 wspr_type = 15;
852 break;
853 case 'o': //use ordered-statistics-decoder
854 ndepth=(int) strtol(optarg,NULL,10);
855 break;
856 case 'q': //no shift jittering
857 quickmode = 1;
858 break;
859 case 's': //single pass mode
860 subtraction = 0;
861 npasses = 1;
862 break;
863 case 'v':
864 verbose = 1;
865 break;
866 case 'w':
867 fmin=-150.0;
868 fmax=150.0;
869 break;
870 case 'z':
871 bias=strtod(optarg,NULL); //fano metric bias (default is 0.45)
872 break;
873 case '?':
874 usage();
875 return 1;
876 }
877 }
878
879 if( access(data_dir, R_OK | W_OK)) {
880 fprintf(stderr, "Error: inaccessible data directory: '%s'\n", data_dir);
881 usage();
882 return EXIT_FAILURE;
883 }
884
885 if( optind+1 > argc) {
886 usage();
887 return 1;
888 } else {
889 ptr_to_infile=argv[optind];
890 }
891
892 if( stackdecoder ) {
893 stack=calloc(stacksize,sizeof(struct snode));
894 }
895
896 // setup metric table
897 for(i=0; i<256; i++) {
898 mettab[0][i]=round( 10*(metric_tables[2][i]-bias) );
899 mettab[1][i]=round( 10*(metric_tables[2][255-i]-bias) );
900 }
901
902 FILE *fp_fftwf_wisdom_file, *fall_wspr, *fwsprd, *fhash, *ftimer;
903 strcpy(wisdom_fname,".");
904 strcpy(all_fname,".");
905 strcpy(spots_fname,".");
906 strcpy(timer_fname,".");
907 strcpy(hash_fname,".");
908 if(data_dir != NULL) {
909 strncpy(wisdom_fname,data_dir, sizeof wisdom_fname);
910 strncpy(all_fname,data_dir, sizeof all_fname);
911 strncpy(spots_fname,data_dir, sizeof spots_fname);
912 strncpy(timer_fname,data_dir, sizeof timer_fname);
913 strncpy(hash_fname,data_dir, sizeof hash_fname);
914 }
915 strncat(wisdom_fname,"/wspr_wisdom.dat",20);
916 strncat(all_fname,"/ALL_WSPR.TXT",20);
917 strncat(spots_fname,"/wspr_spots.txt",20);
918 strncat(timer_fname,"/wspr_timer.out",20);
919 strncat(hash_fname,"/hashtable.txt",20);
920 if ((fp_fftwf_wisdom_file = fopen(wisdom_fname, "r"))) { //Open FFTW wisdom
921 fftwf_import_wisdom_from_file(fp_fftwf_wisdom_file);
922 fclose(fp_fftwf_wisdom_file);
923 }
924
925 fall_wspr=fopen(all_fname,"a");
926 fwsprd=fopen(spots_fname,"w");
927 // FILE *fdiag;
928 // fdiag=fopen("wsprd_diag","a");
929
930 if((ftimer=fopen(timer_fname,"r"))) {
931 //Accumulate timing data
932 int nr=fscanf(ftimer,"%f %f %f %f %f %f %f %f",
933 &treadwav,&tcandidates,&tsync0,&tsync1,&tsync2,&tfano,&tosd,&ttotal);
934 fclose(ftimer);
935 if(nr == 0) fprintf(stderr, "Empty timer file: '%s'\n", timer_fname);
936 }
937 ftimer=fopen(timer_fname,"w");
938
939 if( strstr(ptr_to_infile,".wav") ) {
940 ptr_to_infile_suffix=strstr(ptr_to_infile,".wav");
941
942 t0 = clock();
943 npoints=readwavfile(ptr_to_infile, wspr_type, idat, qdat);
944 treadwav += (float)(clock()-t0)/CLOCKS_PER_SEC;
945
946 if( npoints == 1 ) {
947 return 1;
948 }
949 dialfreq=dialfreq_cmdline - (dialfreq_error*1.0e-06);
950 } else if ( strstr(ptr_to_infile,".c2") !=0 ) {
951 ptr_to_infile_suffix=strstr(ptr_to_infile,".c2");
952 npoints=readc2file(ptr_to_infile, idat, qdat, &dialfreq, &wspr_type);
953 if( npoints == 1 ) {
954 return 1;
955 }
956 dialfreq -= (dialfreq_error*1.0e-06);
957 } else {
958 printf("Error: Failed to open %s\n",ptr_to_infile);
959 printf("WSPR file must have suffix .wav or .c2\n");
960 return 1;
961 }
962
963 // Parse date and time from given filename
964 strncpy(date,ptr_to_infile_suffix-11,6);
965 strncpy(uttime,ptr_to_infile_suffix-4,4);
966 date[6]='\0';
967 uttime[4]='\0';
968
969 // Do windowed ffts over 2 symbols, stepped by half symbols
970 int nffts=4*floor(npoints/512)-1;
971 fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512);
972 fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512);
973 PLAN3 = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE);
974
975 float ps[512][nffts];
976 float w[512];
977 for(i=0; i<512; i++) {
978 w[i]=sin(0.006147931*i);
979 }
980
981 if( usehashtable ) {
982 char line[80], hcall[13], hgrid[5];
983 if( (fhash=fopen(hash_fname,"r+")) ) {
984 while (fgets(line, sizeof(line), fhash) != NULL) {
985 hgrid[0]='\0';
986 sscanf(line,"%d %s %s",&nh,hcall,hgrid);
987 strcpy(hashtab+nh*13,hcall);
988 if(strlen(hgrid)>0) strcpy(loctab+nh*5,hgrid);
989 }
990 } else {
991 fhash=fopen(hash_fname,"w+");
992 }
993 fclose(fhash);
994 }
995
996 //*************** main loop starts here *****************
997 for (ipass=0; ipass<npasses; ipass++) {
998 if(ipass==1 && ndecodes_pass == 0 && npasses>2) ipass=2;
999 if(ipass < 2) {
1000 nblocksize=1;
1001 maxdrift=4;
1002 minsync2=0.12;
1003 }
1004 if(ipass == 2 ) {
1005 nblocksize=4; // try 3 blocksizes plus bitbybit normalization
1006 maxdrift=0; // no drift for smaller frequency estimator variance
1007 minsync2=0.10;
1008 }
1009 ndecodes_pass=0; // still needed?
1010
1011 for (i=0; i<nffts; i++) {
1012 for(j=0; j<512; j++ ) {
1013 k=i*128+j;
1014 fftin[j][0]=idat[k] * w[j];
1015 fftin[j][1]=qdat[k] * w[j];
1016 }
1017 fftwf_execute(PLAN3);
1018 for (j=0; j<512; j++ ) {
1019 k=j+256;
1020 if( k>511 )
1021 k=k-512;
1022 ps[j][i]=fftout[k][0]*fftout[k][0]+fftout[k][1]*fftout[k][1];
1023 }
1024 }
1025
1026 // Compute average spectrum
1027 for (i=0; i<512; i++) psavg[i]=0.0;
1028 for (i=0; i<nffts; i++) {
1029 for (j=0; j<512; j++) {
1030 psavg[j]=psavg[j]+ps[j][i];
1031 }
1032 }
1033
1034 // Smooth with 7-point window and limit spectrum to +/-150 Hz
1035 int window[7]={1,1,1,1,1,1,1};
1036 float smspec[411];
1037 for (i=0; i<411; i++) {
1038 smspec[i]=0.0;
1039 for(j=-3; j<=3; j++) {
1040 k=256-205+i+j;
1041 smspec[i]=smspec[i]+window[j+3]*psavg[k];
1042 }
1043 }
1044
1045 // Sort spectrum values, then pick off noise level as a percentile
1046 float tmpsort[411];
1047 for (j=0; j<411; j++) {
1048 tmpsort[j]=smspec[j];
1049 }
1050 qsort(tmpsort, 411, sizeof(float), floatcomp);
1051
1052 // Noise level of spectrum is estimated as 123/411= 30'th percentile
1053 float noise_level = tmpsort[122];
1054
1055 /* Renormalize spectrum so that (large) peaks represent an estimate of snr.
1056 * We know from experience that threshold snr is near -7dB in wspr bandwidth,
1057 * corresponding to -7-26.3=-33.3dB in 2500 Hz bandwidth.
1058 * The corresponding threshold is -42.3 dB in 2500 Hz bandwidth for WSPR-15. */
1059
1060 float min_snr, snr_scaling_factor;
1061 min_snr = pow(10.0,-8.0/10.0); //this is min snr in wspr bw
1062 if( wspr_type == 2 ) {
1063 snr_scaling_factor=26.3;
1064 } else {
1065 snr_scaling_factor=35.3;
1066 }
1067 for (j=0; j<411; j++) {
1068 smspec[j]=smspec[j]/noise_level - 1.0;
1069 if( smspec[j] < min_snr) smspec[j]=0.1*min_snr;
1070 continue;
1071 }
1072
1073 // Find all local maxima in smoothed spectrum.
1074 for (i=0; i<200; i++) {
1075 candidates[i].freq=0.0;
1076 candidates[i].snr=0.0;
1077 candidates[i].drift=0.0;
1078 candidates[i].shift=0;
1079 candidates[i].sync=0.0;
1080 }
1081
1082 int npk=0;
1083 unsigned char candidate;
1084 for(j=1; j<410; j++) {
1085 candidate = (smspec[j]>smspec[j-1]) &&
1086 (smspec[j]>smspec[j+1]) &&
1087 (npk<200);
1088 if ( candidate ) {
1089 candidates[npk].freq = (j-205)*df;
1090 candidates[npk].snr = 10*log10(smspec[j])-snr_scaling_factor;
1091 npk++;
1092 }
1093 }
1094 if( more_candidates ) {
1095 for(j=0; j<411; j=j+3) {
1096 candidate = (smspec[j]>min_snr) && (npk<200);
1097 if ( candidate ) {
1098 candidates[npk].freq = (j-205)*df;
1099 candidates[npk].snr = 10*log10(smspec[j])-snr_scaling_factor;
1100 npk++;
1101 }
1102 }
1103 }
1104
1105 // Compute corrected fmin, fmax, accounting for dial frequency error
1106 fmin += dialfreq_error; // dialfreq_error is in units of Hz
1107 fmax += dialfreq_error;
1108
1109 // Don't waste time on signals outside of the range [fmin,fmax].
1110 i=0;
1111 for( j=0; j<npk; j++) {
1112 if( candidates[j].freq >= fmin && candidates[j].freq <= fmax ) {
1113 candidates[i]=candidates[j];
1114 i++;
1115 }
1116 }
1117 npk=i;
1118
1119 // bubble sort on snr
1120 int pass;
1121 struct cand tmp;
1122 for (pass = 1; pass <= npk - 1; pass++) {
1123 for (k = 0; k < npk - pass ; k++) {
1124 if (candidates[k].snr < candidates[k+1].snr) {
1125 tmp = candidates[k];
1126 candidates[k]=candidates[k+1];
1127 candidates[k+1] = tmp;
1128 }
1129 }
1130 }
1131
1132 t0=clock();
1133
1134 /* Make coarse estimates of shift (DT), freq, and drift
1135
1136 * Look for time offsets up to +/- 8 symbols (about +/- 5.4 s) relative
1137 to nominal start time, which is 2 seconds into the file
1138
1139 * Calculates shift relative to the beginning of the file
1140
1141 * Negative shifts mean that signal started before start of file
1142
1143 * The program prints DT = shift-2 s
1144
1145 * Shifts that cause sync vector to fall off of either end of the data
1146 vector are accommodated by "partial decoding", such that missing
1147 symbols produce a soft-decision symbol value of 128
1148
1149 * The frequency drift model is linear, deviation of +/- drift/2 over the
1150 span of 162 symbols, with deviation equal to 0 at the center of the
1151 signal vector.
1152 */
1153
1154 int idrift,ifr,if0,ifd,k0;
1155 int kindex;
1156 float smax,ss,pow,p0,p1,p2,p3;
1157 for(j=0; j<npk; j++) { //For each candidate...
1158 smax=-1e30;
1159 if0=candidates[j].freq/df+256;
1160 for (ifr=if0-2; ifr<=if0+2; ifr++) { //Freq search
1161 for( k0=-10; k0<22; k0++) { //Time search
1162 for (idrift=-maxdrift; idrift<=maxdrift; idrift++) { //Drift search
1163 ss=0.0;
1164 pow=0.0;
1165 for (k=0; k<162; k++) { //Sum over symbols
1166 ifd=ifr+((float)k-81.0)/81.0*( (float)idrift )/(2.0*df);
1167 kindex=k0+2*k;
1168 if( kindex < nffts ) {
1169 p0=ps[ifd-3][kindex];
1170 p1=ps[ifd-1][kindex];
1171 p2=ps[ifd+1][kindex];
1172 p3=ps[ifd+3][kindex];
1173
1174 p0=sqrt(p0);
1175 p1=sqrt(p1);
1176 p2=sqrt(p2);
1177 p3=sqrt(p3);
1178
1179 ss=ss+(2*pr3[k]-1)*((p1+p3)-(p0+p2));
1180 pow=pow+p0+p1+p2+p3;
1181 }
1182 }
1183 sync1=ss/pow;
1184 if( sync1 > smax ) { //Save coarse parameters
1185 smax=sync1;
1186 candidates[j].shift=128*(k0+1);
1187 candidates[j].drift=idrift;
1188 candidates[j].freq=(ifr-256)*df;
1189 candidates[j].sync=sync1;
1190 }
1191 }
1192 }
1193 }
1194 }
1195 tcandidates += (float)(clock()-t0)/CLOCKS_PER_SEC;
1196
1197 /*
1198 Refine the estimates of freq, shift using sync as a metric.
1199 Sync is calculated such that it is a float taking values in the range
1200 [0.0,1.0].
1201
1202 Function sync_and_demodulate has three modes of operation
1203 mode is the last argument:
1204
1205 0 = no frequency or drift search. find best time lag.
1206 1 = no time lag or drift search. find best frequency.
1207 2 = no frequency or time lag search. Calculate soft-decision
1208 symbols using passed frequency and shift.
1209
1210 NB: best possibility for OpenMP may be here: several worker threads
1211 could each work on one candidate at a time.
1212 */
1213 for (j=0; j<npk; j++) {
1214
1215 f1=candidates[j].freq;
1216 drift1=candidates[j].drift;
1217 shift1=candidates[j].shift;
1218 sync1=candidates[j].sync;
1219
1220 // coarse-grid lag and freq search, then if sync>minsync1 continue
1221 fstep=0.0; ifmin=0; ifmax=0;
1222 lagmin=shift1-128;
1223 lagmax=shift1+128;
1224 lagstep=64;
1225 t0 = clock();
1226 sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
1227 lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0);
1228 tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC;
1229
1230 fstep=0.25; ifmin=-2; ifmax=2;
1231 t0 = clock();
1232 sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
1233 lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1);
1234
1235 if(ipass < 2) {
1236 // refine drift estimate
1237 fstep=0.0; ifmin=0; ifmax=0;
1238 float driftp,driftm,syncp,syncm;
1239 driftp=drift1+0.5;
1240 sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
1241 lagmin, lagmax, lagstep, &driftp, symfac, &syncp, 1);
1242
1243 driftm=drift1-0.5;
1244 sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
1245 lagmin, lagmax, lagstep, &driftm, symfac, &syncm, 1);
1246
1247 if(syncp>sync1) {
1248 drift1=driftp;
1249 sync1=syncp;
1250 } else if (syncm>sync1) {
1251 drift1=driftm;
1252 sync1=syncm;
1253 }
1254 }
1255 tsync1 += (float)(clock()-t0)/CLOCKS_PER_SEC;
1256
1257 // fine-grid lag and freq search
1258 if( sync1 > minsync1 ) {
1259
1260 lagmin=shift1-32; lagmax=shift1+32; lagstep=16;
1261 t0 = clock();
1262 sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
1263 lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0);
1264 tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC;
1265
1266 // fine search over frequency
1267 fstep=0.05; ifmin=-2; ifmax=2;
1268 t0 = clock();
1269 sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1,
1270 lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1);
1271 tsync1 += (float)(clock()-t0)/CLOCKS_PER_SEC;
1272
1273 candidates[j].freq=f1;
1274 candidates[j].shift=shift1;
1275 candidates[j].drift=drift1;
1276 candidates[j].sync=sync1;
1277 }
1278 }
1279
1280 int nwat=0;
1281 int idupe;
1282 for ( j=0; j<npk; j++) {
1283 idupe=0;
1284 for (k=0;k<nwat;k++) {
1285 if( fabsf(candidates[j].freq - candidates[k].freq) < 0.05 &&
1286 abs(candidates[j].shift - candidates[k].shift) < 16 ) {
1287 idupe=1;
1288 break;
1289 }
1290 }
1291 if( idupe == 1 ) {
1292 if(candidates[j].sync > candidates[k].sync) candidates[k]=candidates[j];
1293 } else if ( candidates[j].sync > minsync2 ) {
1294 candidates[nwat]=candidates[j];
1295 nwat++;
1296 }
1297 }
1298
1299 int idt, ii, jittered_shift;
1300 float y,sq,rms;
1301 int ib, blocksize, bitmetric;
1302 int n1,n2,n3,nadd,nu,ntype;
1303 int osd_decode;
1304 for (j=0; j<nwat; j++) {
1305 memset(symbols,0,sizeof(char)*nbits*2);
1306 memset(callsign,0,sizeof(char)*13);
1307 memset(grid,0,sizeof(char)*5);
1308 memset(call_loc_pow,0,sizeof(char)*23);
1309 f1=candidates[j].freq;
1310 shift1=candidates[j].shift;
1311 drift1=candidates[j].drift;
1312 not_decoded=1;
1313 osd_decode=0;
1314
1315 ib=1;
1316 while( ib <= nblocksize && not_decoded ) {
1317 if (ib < 4) { blocksize=ib; bitmetric=0; }
1318 if (ib == 4) { blocksize=1; bitmetric=1; }
1319
1320 idt=0; ii=0;
1321 while ( not_decoded && idt<=(128/iifac)) {
1322 ii=(idt+1)/2;
1323 if( idt%2 == 1 ) ii=-ii;
1324 ii=iifac*ii;
1325 jittered_shift=shift1+ii;
1326 nhardmin=0; dmin=0.0;
1327
1328 // Get soft-decision symbols
1329 t0 = clock();
1330 noncoherent_sequence_detection(idat, qdat, npoints, symbols, &f1,
1331 &jittered_shift, &drift1, symfac, &blocksize, &bitmetric);
1332 tsync2 += (float)(clock()-t0)/CLOCKS_PER_SEC;
1333
1334 sq=0.0;
1335 for(i=0; i<162; i++) {
1336 y=(float)symbols[i] - 128.0;
1337 sq += y*y;
1338 }
1339 rms=sqrt(sq/162.0);
1340
1341 if(rms > minrms) {
1342 deinterleave(symbols);
1343 t0 = clock();
1344
1345 if ( stack ) {
1346 not_decoded = jelinek(&metric, &cycles, decdata, symbols, nbits,
1347 stacksize, stack, mettab,maxcycles);
1348 } else {
1349 not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
1350 mettab,delta,maxcycles);
1351 }
1352
1353 tfano += (float)(clock()-t0)/CLOCKS_PER_SEC;
1354
1355 if( (ndepth >= 0) && not_decoded ) {
1356 for(i=0; i<162; i++) {
1357 fsymbs[i]=symbols[i]-128.0;
1358 }
1359 t0 = clock();
1360 osdwspr_(fsymbs,apmask,&ndepth,cw,&nhardmin,&dmin);
1361 tosd += (float)(clock()-t0)/CLOCKS_PER_SEC;
1362
1363 for(i=0; i<162; i++) {
1364 symbols[i]=255*cw[i];
1365 }
1366 fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
1367 mettab,delta,maxcycles);
1368 for(i=0; i<11; i++) {
1369 if( decdata[i]>127 ) {
1370 message[i]=decdata[i]-256;
1371 } else {
1372 message[i]=decdata[i];
1373 }
1374 }
1375 unpack50(message,&n1,&n2);
1376 if( !unpackcall(n1,callsign) ) break;
1377 callsign[12]=0;
1378 if( !unpackgrid(n2, grid) ) break;
1379 grid[4]=0;
1380 ntype = (n2&127) - 64;
1381 int itype;
1382 if( (ntype >= 0) && (ntype <= 62) ) {
1383 nu = ntype%10;
1384 itype=1;
1385 if( !(nu == 0 || nu == 3 || nu == 7) ) {
1386 nadd=nu;
1387 if( nu > 3 ) nadd=nu-3;
1388 if( nu > 7 ) nadd=nu-7;
1389 n3=n2/128+32768*(nadd-1);
1390 if( !unpackpfx(n3,callsign) ) {
1391 break;
1392 }
1393 itype=2;
1394 }
1395 ihash=nhash(callsign,strlen(callsign),(uint32_t)146);
1396 if(strncmp(hashtab+ihash*13,callsign,13)==0) {
1397 if( (itype==1 && strncmp(loctab+ihash*5,grid,5)==0) ||
1398 (itype==2) ) {
1399 not_decoded=0;
1400 osd_decode =1;
1401 }
1402 }
1403 }
1404 }
1405
1406 }
1407 idt++;
1408 if( quickmode ) break;
1409 }
1410 ib++;
1411 }
1412
1413 if( !not_decoded ) {
1414 ndecodes_pass++;
1415
1416 for(i=0; i<11; i++) {
1417
1418 if( decdata[i]>127 ) {
1419 message[i]=decdata[i]-256;
1420 } else {
1421 message[i]=decdata[i];
1422 }
1423
1424 }
1425
1426 // Unpack the decoded message, update the hashtable, apply
1427 // sanity checks on grid and power, and return
1428 // call_loc_pow string and also callsign (for de-duping).
1429 noprint=unpk_(message,hashtab,loctab,call_loc_pow,callsign);
1430 if( subtraction && !noprint ) {
1431 if( get_wspr_channel_symbols(call_loc_pow, hashtab, loctab, channel_symbols) ) {
1432 subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
1433 if(!osd_decode) nhardmin=count_hard_errors(symbols,channel_symbols);
1434 } else {
1435 break;
1436 }
1437 }
1438
1439 // Remove dupes (same callsign and freq within 4 Hz)
1440 int dupe=0;
1441 for (i=0; i<uniques; i++) {
1442 if(!strcmp(callsign,allcalls[i]) &&
1443 (fabs(f1-allfreqs[i]) <4.0)) dupe=1;
1444 }
1445 if( (verbose || !dupe) && !noprint) {
1446 strcpy(allcalls[uniques],callsign);
1447 allfreqs[uniques]=f1;
1448 uniques++;
1449
1450 // Add an extra space at the end of each line so that wspr-x doesn't
1451 // truncate the power (TNX to DL8FCL!)
1452
1453 if( wspr_type == 15 ) {
1454 freq_print=dialfreq+(1500+112.5+f1/8.0)/1e6;
1455 dt_print=shift1*8*dt-1.0;
1456 } else {
1457 freq_print=dialfreq+(1500+f1)/1e6;
1458 dt_print=shift1*dt-1.0;
1459 }
1460
1461 strcpy(decodes[uniques-1].date,date);
1462 strcpy(decodes[uniques-1].time,uttime);
1463 decodes[uniques-1].sync=candidates[j].sync;
1464 decodes[uniques-1].snr=candidates[j].snr;
1465 decodes[uniques-1].dt=dt_print;
1466 decodes[uniques-1].freq=freq_print;
1467 strcpy(decodes[uniques-1].message,call_loc_pow);
1468 decodes[uniques-1].drift=drift1;
1469 decodes[uniques-1].cycles=cycles;
1470 decodes[uniques-1].jitter=ii;
1471 decodes[uniques-1].blocksize=blocksize+3*bitmetric;
1472 decodes[uniques-1].metric=metric;
1473 decodes[uniques-1].nhardmin=nhardmin;
1474 decodes[uniques-1].ipass=ipass;
1475 decodes[uniques-1].decodetype=osd_decode;
1476 }
1477 }
1478 }
1479
1480 if( ipass == 0 && writec2 ) {
1481 char c2filename[15];
1482 double carrierfreq=dialfreq;
1483 int wsprtype=2;
1484 strcpy(c2filename,"000000_0001.c2");
1485 printf("Writing %s\n",c2filename);
1486 writec2file(c2filename, wsprtype, carrierfreq, idat, qdat);
1487 }
1488 }
1489
1490 // sort the result in order of increasing frequency
1491 struct result temp;
1492 for (j = 1; j <= uniques - 1; j++) {
1493 for (k = 0; k < uniques - j ; k++) {
1494 if (decodes[k].freq > decodes[k+1].freq) {
1495 temp = decodes[k];
1496 decodes[k]=decodes[k+1];;
1497 decodes[k+1] = temp;
1498 }
1499 }
1500 }
1501
1502 for (i=0; i<uniques; i++) {
1503 printf("%4s %3.0f %4.1f %10.6f %2d %-s \n",
1504 decodes[i].time, decodes[i].snr,decodes[i].dt, decodes[i].freq,
1505 (int)decodes[i].drift, decodes[i].message);
1506 fprintf(fall_wspr,
1507 "%6s %4s %3.0f %5.2f %11.7f %-22s %2d %5.2f %2d %2d %4d %2d %3d %5u %5d\n",
1508 decodes[i].date, decodes[i].time, decodes[i].snr,
1509 decodes[i].dt, decodes[i].freq, decodes[i].message,
1510 (int)decodes[i].drift, decodes[i].sync,
1511 decodes[i].ipass+1,decodes[i].blocksize,decodes[i].jitter,
1512 decodes[i].decodetype,decodes[i].nhardmin,decodes[i].cycles/81,
1513 decodes[i].metric);
1514 fprintf(fwsprd,
1515 "%6s %4s %3d %3.0f %4.1f %10.6f %-22s %2d %5u %4d\n",
1516 decodes[i].date, decodes[i].time, (int)(10*decodes[i].sync),
1517 decodes[i].snr, decodes[i].dt, decodes[i].freq,
1518 decodes[i].message, (int)decodes[i].drift, decodes[i].cycles/81,
1519 decodes[i].jitter);
1520
1521 }
1522 printf("<DecodeFinished>\n");
1523
1524 fftwf_free(fftin);
1525 fftwf_free(fftout);
1526
1527 if ((fp_fftwf_wisdom_file = fopen(wisdom_fname, "w"))) {
1528 fftwf_export_wisdom_to_file(fp_fftwf_wisdom_file);
1529 fclose(fp_fftwf_wisdom_file);
1530 }
1531
1532 ttotal += (float)(clock()-t00)/CLOCKS_PER_SEC;
1533
1534 fprintf(ftimer,"%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n\n",
1535 treadwav,tcandidates,tsync0,tsync1,tsync2,tfano,tosd,ttotal);
1536
1537 fprintf(ftimer,"Code segment Seconds Frac\n");
1538 fprintf(ftimer,"-----------------------------------\n");
1539 fprintf(ftimer,"readwavfile %7.2f %7.2f\n",treadwav,treadwav/ttotal);
1540 fprintf(ftimer,"Coarse DT f0 f1 %7.2f %7.2f\n",tcandidates,
1541 tcandidates/ttotal);
1542 fprintf(ftimer,"sync_and_demod(0) %7.2f %7.2f\n",tsync0,tsync0/ttotal);
1543 fprintf(ftimer,"sync_and_demod(1) %7.2f %7.2f\n",tsync1,tsync1/ttotal);
1544 fprintf(ftimer,"sync_and_demod(2) %7.2f %7.2f\n",tsync2,tsync2/ttotal);
1545 fprintf(ftimer,"Stack/Fano decoder %7.2f %7.2f\n",tfano,tfano/ttotal);
1546 fprintf(ftimer,"OSD decoder %7.2f %7.2f\n",tosd,tosd/ttotal);
1547 fprintf(ftimer,"-----------------------------------\n");
1548 fprintf(ftimer,"Total %7.2f %7.2f\n",ttotal,1.0);
1549
1550 fclose(fall_wspr);
1551 fclose(fwsprd);
1552 // fclose(fdiag);
1553 fclose(ftimer);
1554 fftwf_destroy_plan(PLAN1);
1555 fftwf_destroy_plan(PLAN2);
1556 fftwf_destroy_plan(PLAN3);
1557
1558 if( usehashtable ) {
1559 fhash=fopen(hash_fname,"w");
1560 for (i=0; i<32768; i++) {
1561 if( strncmp(hashtab+i*13,"\0",1) != 0 ) {
1562 fprintf(fhash,"%5d %s %s\n",i,hashtab+i*13,loctab+i*5);
1563 }
1564 }
1565 fclose(fhash);
1566 }
1567
1568 free(hashtab);
1569 free(loctab);
1570 free(symbols);
1571 free(decdata);
1572 free(channel_symbols);
1573 free(callsign);
1574 free(call_loc_pow);
1575 free(idat);
1576 free(qdat);
1577 free(stack);
1578
1579 return 0;
1580 }
1581