From 1a167c87bf4c93dd585fcd15eeb358486dc60be9 Mon Sep 17 00:00:00 2001 From: Guenael Date: Wed, 22 Dec 2021 21:22:06 -0500 Subject: [PATCH] WIP chore: wsprd lib update part2 --- Makefile | 4 +- wsprd/metric_tables.h | 5 - wsprd/wsprd.c | 602 +++++++++++++++--------------------------- wsprd/wsprd.h | 30 --- 4 files changed, 216 insertions(+), 425 deletions(-) diff --git a/Makefile b/Makefile index 14b689b..2672352 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ CC = clang -CFLAGS = -O3 -std=gnu17 -LIBS = -lusb-1.0 -lrtlsdr -lpthread -lfftw3f -lcurl -lm +CFLAGS = -O3 -std=gnu17 -Wall +LIBS = -lusb-1.0 -lrtlsdr -lpthread -lfftw3f -lcurl -lm -Wall # Note # gcc is a bit faster that clang on this app diff --git a/wsprd/metric_tables.h b/wsprd/metric_tables.h index a86a956..7db6e27 100644 --- a/wsprd/metric_tables.h +++ b/wsprd/metric_tables.h @@ -4,9 +4,6 @@ * should be normalized to have rms amplitude equal to "symbol_scale". ********************************************************************************/ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wmissing-braces" - //float symbol_scale[5]={42.6, 53.3, 72.7, 100.2, 125.4}; float metric_tables[5][256]={ {0.9782, 0.9695, 0.9689, 0.9669, 0.9666, 0.9653, 0.9638, 0.9618, 0.9599, 0.9601, @@ -140,5 +137,3 @@ float metric_tables[5][256]={ -17.5973, -17.0403, -17.7039, -18.0073, -18.1840, -18.3848, -18.6286, -20.7063, 1.43370769e-019, 2.64031087e-006, 6.6908396e+031, 1.77537994e+028, 2.79322819e+020, 1.94326e-019, 0.00019371575, 2.80722121e-041}}; - -#pragma GCC diagnostic pop \ No newline at end of file diff --git a/wsprd/wsprd.c b/wsprd/wsprd.c index f4ea1b5..1843a04 100644 --- a/wsprd/wsprd.c +++ b/wsprd/wsprd.c @@ -56,11 +56,14 @@ #define DF15 DF * 1.5 #define TWOPIDT 2.0 * M_PI * DT -#define NFILT 256 -#define NSIG NSYM * NSPERSYM - -/* Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT, FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE */ +/* Possible PATIENCE options: F + FTW_ESTIMATE, + FFTW_ESTIMATE_PATIENT, + FFTW_MEASURE, + FFTW_PATIENT, + FFTW_EXHAUSTIVE +*/ #define PATIENCE FFTW_ESTIMATE fftwf_plan PLAN; @@ -78,91 +81,88 @@ uint8_t pr3vector[NSYM] = { 0, 0}; -//*************************************************************************** -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. * - ************************************************************************/ +/* 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. + */ +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) { - static float fplast = -10000.0; - float f0 = 0.0, fp, ss, - fbest = 0.0, - fsum = 0.0, - f2sum = 0.0, - fsymb[NSYM]; - int best_shift = 0; float i0[NSYM], q0[NSYM], - i1[NSYM], q1[NSYM], - i2[NSYM], q2[NSYM], + i1[NSYM], q1[NSYM], + i2[NSYM], q2[NSYM], i3[NSYM], q3[NSYM]; - float p0, p1, p2, p3, cmet, totp, syncmax, fac; - float c0[NSPERSYM], s0[NSPERSYM], - c1[NSPERSYM], s1[NSPERSYM], - c2[NSPERSYM], s2[NSPERSYM], + float c0[NSPERSYM], s0[NSPERSYM], + c1[NSPERSYM], s1[NSPERSYM], + c2[NSPERSYM], s2[NSPERSYM], c3[NSPERSYM], s3[NSPERSYM]; - float dphi0, cdphi0, sdphi0, - dphi1, cdphi1, sdphi1, - dphi2, cdphi2, sdphi2, - dphi3, cdphi3, sdphi3; + float fsymb[NSYM]; + + float fbest = 0.0, + fsum = 0.0, + f2sum = 0.0; + + int best_shift = 0; + static float fplast = -10000.0; + float syncmax = -1e30; - syncmax = -1e30; if (mode == 0) { ifmin = 0; ifmax = 0; fstep = 0.0; - f0 = *f1; - } - if (mode == 1) { + } else if (mode == 1) { lagmin = *shift1; lagmax = *shift1; - f0 = *f1; - } - if (mode == 2) { + } else if (mode == 2) { lagmin = *shift1; lagmax = *shift1; - ifmin = 0; - ifmax = 0; - f0 = *f1; + ifmin = 0; + ifmax = 0; } for (int ifreq = ifmin; ifreq <= ifmax; ifreq++) { - f0 = *f1 + ifreq * fstep; + float f0 = *f1 + ifreq * fstep; for (int lag = lagmin; lag <= lagmax; lag = lag + lagstep) { - ss = 0.0; - totp = 0.0; + float ss = 0.0; + float totp = 0.0; for (int i = 0; i < NSYM; i++) { - fp = f0 + (*drift1 / 2.0) * ((float)i - (float)NBITS) / (float)NBITS; + float fp = f0 + (*drift1 / 2.0) * ((float)i - (float)NBITS) / (float)NBITS; if (i == 0 || (fp != fplast)) { // only calculate sin/cos if necessary - dphi0 = TWOPIDT * (fp - DF15); - cdphi0 = cos(dphi0); - sdphi0 = sin(dphi0); + float dphi0 = TWOPIDT * (fp - DF15); + float cdphi0 = cosf(dphi0); + float sdphi0 = sinf(dphi0); - dphi1 = TWOPIDT * (fp - DF05); - cdphi1 = cos(dphi1); - sdphi1 = sin(dphi1); + float dphi1 = TWOPIDT * (fp - DF05); + float cdphi1 = cosf(dphi1); + float sdphi1 = sinf(dphi1); - dphi2 = TWOPIDT * (fp + DF05); - cdphi2 = cos(dphi2); - sdphi2 = sin(dphi2); + float dphi2 = TWOPIDT * (fp + DF05); + float cdphi2 = cosf(dphi2); + float sdphi2 = sinf(dphi2); - dphi3 = TWOPIDT * (fp + DF15); - cdphi3 = cos(dphi3); - sdphi3 = sin(dphi3); + float dphi3 = TWOPIDT * (fp + DF15); + float cdphi3 = cosf(dphi3); + float sdphi3 = sinf(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; + 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 (int j = 1; j < NSPERSYM; j++) { c0[j] = c0[j - 1] * cdphi0 - s0[j - 1] * sdphi0; @@ -177,14 +177,10 @@ void sync_and_demodulate(float *id, float *qd, long np, 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; + 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 (int j = 0; j < NSPERSYM; j++) { int k = lag + i * NSPERSYM + j; @@ -199,18 +195,14 @@ void sync_and_demodulate(float *id, float *qd, long np, q3[i] = q3[i] - id[k] * s3[j] + qd[k] * c3[j]; } } - p0 = i0[i] * i0[i] + q0[i] * q0[i]; - p1 = i1[i] * i1[i] + q1[i] * q1[i]; - p2 = i2[i] * i2[i] + q2[i] * q2[i]; - p3 = i3[i] * i3[i] + q3[i] * q3[i]; - p0 = sqrt(p0); - p1 = sqrt(p1); - p2 = sqrt(p2); - p3 = sqrt(p3); + float p0 = sqrt(i0[i] * i0[i] + q0[i] * q0[i]); + float p1 = sqrt(i1[i] * i1[i] + q1[i] * q1[i]); + float p2 = sqrt(i2[i] * i2[i] + q2[i] * q2[i]); + float p3 = sqrt(i3[i] * i3[i] + q3[i] * q3[i]); totp = totp + p0 + p1 + p2 + p3; - cmet = (p1 + p3) - (p0 + p2); + float cmet = (p1 + p3) - (p0 + p2); ss = (pr3vector[i] == 1) ? ss + cmet : ss - cmet; if (mode == 2) { // Compute soft symbols if (pr3vector[i] == 1) { @@ -227,7 +219,7 @@ void sync_and_demodulate(float *id, float *qd, long np, fbest = f0; } } // lag loop - } // freq loop + } // freq loop if (mode <= 1) { // Send best params back to caller *sync = syncmax; @@ -239,10 +231,10 @@ void sync_and_demodulate(float *id, float *qd, long np, if (mode == 2) { *sync = syncmax; for (int i = 0; i < NSYM; i++) { // Normalize the soft symbols - fsum = fsum + fsymb[i] / NSYM; - f2sum = f2sum + fsymb[i] * fsymb[i] / NSYM; + fsum += fsymb[i] / NSYM; + f2sum += fsymb[i] * fsymb[i] / NSYM; } - fac = sqrt(f2sum - fsum * fsum); + float fac = sqrt(f2sum - fsum * fsum); for (int i = 0; i < NSYM; i++) { fsymb[i] = symfac * fsymb[i] / fac; if (fsymb[i] > 127) fsymb[i] = 127.0; @@ -254,182 +246,24 @@ void sync_and_demodulate(float *id, float *qd, long np, return; } -void noncoherent_sequence_detection(float *id, float *qd, long np, - unsigned char *symbols, float *f1, int *shift1, - float *drift1, int symfac, int *nblocksize, int *bitmetric) { - /************************************************************************ - * 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; - int i, j, k, lag, itone, ib, b, nblock, nseq, imask; - float xi[512], xq[512]; - float is[4][NSYM], qs[4][NSYM], - cf[4][NSYM], sf[4][NSYM], - cm, sm, cmp, smp; - float p[512], fac, xm1, xm0; - float c0[NSPERSYM+1], s0[NSPERSYM+1], - c1[NSPERSYM+1], s1[NSPERSYM+1], - c2[NSPERSYM+1], s2[NSPERSYM+1], - c3[NSPERSYM+1], s3[NSPERSYM+1]; - float dphi0, cdphi0, sdphi0, - dphi1, cdphi1, sdphi1, - dphi2, cdphi2, sdphi2, - dphi3, cdphi3, sdphi3; - float f0, fp, fsum = 0.0, f2sum = 0.0, fsymb[NSYM]; - - f0 = *f1; - lag = *shift1; - nblock = *nblocksize; - nseq = 1 << nblock; - int bitbybit = *bitmetric; - - for (i = 0; i < NSYM; i++) { - fp = f0 + (*drift1 / 2.0) * ((float)i - (float)NBITS) / (float)NBITS; - 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 < (NSPERSYM+1); 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; - } - - cf[0][i] = c0[NSPERSYM]; - sf[0][i] = s0[NSPERSYM]; - cf[1][i] = c1[NSPERSYM]; - sf[1][i] = s1[NSPERSYM]; - cf[2][i] = c2[NSPERSYM]; - sf[2][i] = s2[NSPERSYM]; - cf[3][i] = c3[NSPERSYM]; - sf[3][i] = s3[NSPERSYM]; - - is[0][i] = 0.0; - qs[0][i] = 0.0; - is[1][i] = 0.0; - qs[1][i] = 0.0; - is[2][i] = 0.0; - qs[2][i] = 0.0; - is[3][i] = 0.0; - qs[3][i] = 0.0; - - for (j = 0; j < NSPERSYM; j++) { - k = lag + i * NSPERSYM + j; - if ((k > 0) && (k < np)) { - is[0][i] = is[0][i] + id[k] * c0[j] + qd[k] * s0[j]; - qs[0][i] = qs[0][i] - id[k] * s0[j] + qd[k] * c0[j]; - is[1][i] = is[1][i] + id[k] * c1[j] + qd[k] * s1[j]; - qs[1][i] = qs[1][i] - id[k] * s1[j] + qd[k] * c1[j]; - is[2][i] = is[2][i] + id[k] * c2[j] + qd[k] * s2[j]; - qs[2][i] = qs[2][i] - id[k] * s2[j] + qd[k] * c2[j]; - is[3][i] = is[3][i] + id[k] * c3[j] + qd[k] * s3[j]; - qs[3][i] = qs[3][i] - id[k] * s3[j] + qd[k] * c3[j]; - } - } - } - - for (i = 0; i < NSYM; i = i + nblock) { - for (j = 0; j < nseq; j++) { - xi[j] = 0.0; - xq[j] = 0.0; - cm = 1; - sm = 0; - for (ib = 0; ib < nblock; ib++) { - b = (j & (1 << (nblock - 1 - ib))) >> (nblock - 1 - ib); - itone = pr3vector[i + ib] + 2 * b; - xi[j] = xi[j] + is[itone][i + ib] * cm + qs[itone][i + ib] * sm; - xq[j] = xq[j] + qs[itone][i + ib] * cm - is[itone][i + ib] * sm; - cmp = cf[itone][i + ib] * cm - sf[itone][i + ib] * sm; - smp = sf[itone][i + ib] * cm + cf[itone][i + ib] * sm; - cm = cmp; - sm = smp; - } - p[j] = xi[j] * xi[j] + xq[j] * xq[j]; - p[j] = sqrt(p[j]); - } - for (ib = 0; ib < nblock; ib++) { - imask = 1 << (nblock - 1 - ib); - xm1 = 0.0; - xm0 = 0.0; - for (j = 0; j < nseq; j++) { - if ((j & imask) != 0) { - if (p[j] > xm1) xm1 = p[j]; - } - if ((j & imask) == 0) { - if (p[j] > xm0) xm0 = p[j]; - } - } - fsymb[i + ib] = xm1 - xm0; - if (bitbybit == 1) { - fsymb[i + ib] = fsymb[i + ib] / (xm1 > xm0 ? xm1 : xm0); - } - } - } - for (i = 0; i < NSYM; i++) { // Normalize the soft symbols - fsum = fsum + fsymb[i] / NSYM; - f2sum = f2sum + fsymb[i] * fsymb[i] / NSYM; - } - fac = sqrt(f2sum - fsum * fsum); - for (i = 0; i < NSYM; 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) { +/* 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 c0[NSPERSYM], s0[NSPERSYM]; - float dphi, cdphi, sdphi; for (int i = 0; i < NSYM; i++) { float fp = f0 + ((float)drift0 / 2.0) * ((float)i - (float)NBITS) / (float)NBITS; - dphi = TWOPIDT * (fp + ((float)channel_symbols[i] - 1.5) * DF); - cdphi = cos(dphi); - sdphi = sin(dphi); + float dphi = TWOPIDT * (fp + ((float)channel_symbols[i] - 1.5) * DF); + float cdphi = cosf(dphi); + float sdphi = sinf(dphi); c0[0] = 1; s0[0] = 0; @@ -451,7 +285,6 @@ void subtract_signal(float *id, float *qd, long np, } // subtract the signal here. - i0 = i0 / (float)NSPERSYM; // will be wrong for partial symbols at the edges... q0 = q0 / (float)NSPERSYM; @@ -467,9 +300,7 @@ void subtract_signal(float *id, float *qd, long np, } -/****************************************************************************** - Subtract the coherent component of a signal - *******************************************************************************/ +/* Subtract the coherent component of a signal */ void subtract_signal2(float *id, float *qd, long np, @@ -478,7 +309,7 @@ void subtract_signal2(float *id, float drift0, unsigned char *channel_symbols) { - float phi = 0; + float phi = 0.0; const int nfilt = 360; // nfilt must be even number. float refi[SIGNAL_SAMPLES] = {0}, refq[SIGNAL_SAMPLES] = {0}, @@ -493,8 +324,7 @@ void subtract_signal2(float *id, Multiply r(t) by c(t) and subtract from s(t), i.e. s'(t)=s(t)-c(t)r(t) *******************************************************************************/ - // create reference wspr signal vector, centered on f0. - // + /* create reference wspr signal vector, centered on f0. */ for (int i = 0; i < NSYM; i++) { float cs = (float)channel_symbols[i]; @@ -502,20 +332,20 @@ void subtract_signal2(float *id, for (int j = 0; j < NSPERSYM; j++) { int ii = NSPERSYM * i + j; - refi[ii] = cos(phi); // cannot precompute sin/cos because dphi is changing - refq[ii] = sin(phi); + refi[ii] = cosf(phi); // cannot precompute sin/cos because dphi is changing + refq[ii] = sinf(phi); phi = phi + dphi; } } float w[nfilt], norm = 0, partialsum[nfilt]; - // lowpass filter and remove startup transient + /* lowpass filter and remove startup transient */ for (int i = 0; i < nfilt; i++) { partialsum[i] = 0.0; } for (int i = 0; i < nfilt; i++) { - w[i] = sin(M_PI * (float)i / (float)(nfilt - 1)); + w[i] = sinf(M_PI * (float)i / (float)(nfilt - 1)); norm = norm + w[i]; } for (int i = 0; i < nfilt; i++) { @@ -571,72 +401,54 @@ void subtract_signal2(float *id, } -unsigned int count_hard_errors(unsigned char *symbols, unsigned char *channel_symbols) { - int i, is; - unsigned char cw[162]; - unsigned int nerrors; - for (i = 0; i < 162; i++) { - cw[i] = channel_symbols[i] >= 2 ? 1 : 0; - } - deinterleave(cw); - nerrors = 0; - for (i = 0; i < 162; i++) { - is = symbols[i] > 127 ? 1 : 0; - nerrors = nerrors + (is == cw[i] ? 0 : 1); - } - return nerrors; -} - - int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, struct decoder_options options, struct decoder_results *decodes, int32_t *n_results) { - int32_t i, j, k; + uint8_t symbols[NBITS * 2] = {0}; uint8_t decdata[(NBITS + 7) / 8] = {0}; - int8_t message[12] = {0}; + int8_t message[12] = {0}; - float freq0[200], freq1 = 0.0; - float drift0[200], drift1 = 0.0; - float sync0[200], sync1 = 0.0; - float snr0[200]; + float freq0[200], freq1 = 0.0; + float drift0[200], drift1 = 0.0; + float sync0[200], sync1 = 0.0; + float snr0[200]; int32_t shift0[200], shift1 = 0; - int32_t ifmin, ifmax; - - char callsign[13] = {0}; - char call_loc_pow[23] = {0}; - char call[13] = {0}; - char loc[7] = {0}; - char pwr[3] = {0}; uint32_t metric, cycles, maxnp; int32_t worth_a_try; int32_t uniques = 0; + + // Search tuning parameters + float fstep; + int lagmin; + int lagmax; + int lagstep; + int ifmin; + int ifmax; float fmin = -110.0; float fmax = 110.0; - // Hash table - char hashtab[32768 * 13] = {0}; - char loctab[32768 * 5] = {0}; - int32_t nh; - // Parameters used for performance-tuning: - uint32_t maxcycles = 10000; // Fano timeout limit - double minsync1 = 0.10; // First sync limit - double minsync2 = 0.12; // Second sync limit - int32_t iifac = 3; // Step size in final DT peakup - int32_t symfac = 50; // Soft-symbol normalizing factor - int32_t maxdrift = 4; // Maximum (+/-) drift - double minrms = 52.0 * (symfac / 64.0); // Final test for plausible decoding - int32_t delta = 60; // Fano threshold step + uint32_t maxcycles = 10000; // Fano timeout limit + double minsync1 = 0.10; // First sync limit + double minsync2 = 0.12; // Second sync limit + int32_t iifac = 3; // Step size in final DT peakup + int32_t symfac = 50; // Soft-symbol normalizing factor + int32_t maxdrift = 4; // Maximum (+/-) drift + double minrms = 52.0 * (symfac / 64.0); // Final test for plausible decoding + int32_t delta = 60; // Fano threshold step // Results - float allfreqs[100]; - char allcalls[100][13]; - memset(allfreqs, 0, sizeof(float) * 100); - memset(allcalls, 0, sizeof(char) * 100 * 13); + float allfreqs[100] = {0}; + char allcalls[100][13] = {0}; + char callsign[13] = {0}; + char call_loc_pow[23] = {0}; + char call[13] = {0}; + char loc[7] = {0}; + char pwr[3] = {0}; - // Setup metric table + /* Setup metric table */ int32_t mettab[2][256]; float bias = 0.42; for (int i = 0; i < 256; i++) { @@ -644,57 +456,78 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, mettab[1][i] = round(10 * (metric_tables[2][255 - i] - bias)); } - // FFT buffer - fftwf_complex *fftin, *fftout; - - FILE *fp_fftw_wisdom_file, *fhash; - if ((fp_fftw_wisdom_file = fopen("fftw_wisdom.dat", "r"))) { // Open FFTW wisdom - fftwf_import_wisdom_from_file(fp_fftw_wisdom_file); - fclose(fp_fftw_wisdom_file); - } - - // Do windowed ffts over 2 symbols, stepped by half symbols - int32_t nffts = 4 * floor(npoints / 512) - 1; - fftin = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * 512); - fftout = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * 512); - PLAN = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE); - - float ps[512][nffts]; - float w[512]; - for (int i = 0; i < 512; i++) { - w[i] = sinf(0.006147931 * i); - } + /* Setup/Load hash tables */ + FILE *fhash; + int nh; + char hashtab[32768 * 13] = {0}; + char loctab[32768 * 5] = {0}; if (options.usehashtable) { - char line[80], hcall[12]; + char line[80], hcall[12], hgrid[5];; if ((fhash = fopen("hashtable.txt", "r+"))) { while (fgets(line, sizeof(line), fhash) != NULL) { - sscanf(line, "%d %s", &nh, hcall); + hgrid[0] = '\0'; + sscanf(line, "%d %s %s", &nh, hcall, hgrid); strcpy(hashtab + nh * 13, hcall); + if (strlen(hgrid) > 0) strcpy(loctab + nh * 5, hgrid); } fclose(fhash); } } - // Main loop starts here + /* FFT buffer (512 bins) */ + fftwf_complex *fftin, *fftout; + fftin = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * 512); + fftout = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * 512); + PLAN = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE); + + /* Recover FFTW optimization settings */ + FILE *fp_fftw_wisdom_file; + if ((fp_fftw_wisdom_file = fopen("fftw_wisdom.dat", "r"))) { // Open FFTW wisdom + fftwf_import_wisdom_from_file(fp_fftw_wisdom_file); + fclose(fp_fftw_wisdom_file); + } + + /* Hann function */ + float hann[512]; + for (int i = 0; i < 512; i++) { + hann[i] = sinf(0.006147931 * i); + } + + /* FFT output alloc */ + const int blocks = 4 * floor(npoints / 512) - 1; + float ps[512][blocks]; + memset(ps, 0.0, sizeof(float) * 512 * blocks); + + /* Main loop starts here */ for (int ipass = 0; ipass < options.npasses; ipass++) { if (ipass == 1 && uniques == 0) break; - if (ipass == 1) // otherwise we bog down on the second pass - options.quickmode = 1; + if (ipass < 2) { + maxdrift = 4; + minsync2 = 0.12; + } + if (ipass == 2) { + maxdrift = 0; // no drift for smaller frequency estimator variance + minsync2 = 0.10; + } - memset(ps, 0.0, sizeof(float) * 512 * nffts); - for (int i = 0; i < nffts; i++) { - for (j = 0; j < 512; j++) { - k = i * 128 + j; - fftin[j][0] = idat[k] * w[j]; - fftin[j][1] = qdat[k] * w[j]; + /* Compute FFT + * FFT over 2 symbols, stepped by half symbols + */ + for (int i = 0; i < blocks; i++) { + /* Load samples */ + for (int j = 0; j < 512; j++) { + int k = i * 128 + j; + fftin[j][0] = idat[k] * hann[j]; + fftin[j][1] = qdat[k] * hann[j]; } fftwf_execute(PLAN); + /* Recover frequencies */ for (int j = 0; j < 512; j++) { - k = j + 256; + int k = j + 256; if (k > 511) k = k - 512; ps[j][i] = fftout[k][0] * fftout[k][0] + fftout[k][1] * fftout[k][1]; @@ -703,20 +536,21 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, // Compute average spectrum float psavg[512] = {0}; - for (int i = 0; i < nffts; i++) { - for (j = 0; j < 512; j++) { - psavg[j] = psavg[j] + ps[j][i]; + for (int i = 0; i < blocks; i++) { + for (int j = 0; j < 512; j++) { + psavg[j] += ps[j][i]; } } + // Already restricted by previous FIR // Smooth with 7-point window and limit spectrum to +/-150 Hz int32_t window[7] = {1, 1, 1, 1, 1, 1, 1}; float smspec[411]; for (int i = 0; i < 411; i++) { smspec[i] = 0.0; - for (j = -3; j <= 3; j++) { - k = 256 - 205 + i + j; - smspec[i] = smspec[i] + window[j + 3] * psavg[k]; + for (int j = -3; j <= 3; j++) { + int k = 256 - 205 + i + j; + smspec[i] += window[j + 3] * psavg[k]; } } @@ -746,30 +580,26 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, // Find all local maxima in smoothed spectrum. for (int i = 0; i < 200; i++) { - freq0[i] = 0.0; - snr0[i] = 0.0; + freq0[i] = 0.0; + snr0[i] = 0.0; drift0[i] = 0.0; shift0[i] = 0; - sync0[i] = 0.0; + sync0[i] = 0.0; } int32_t npk = 0; for (int j = 1; j < 410; j++) { - if ((smspec[j] > smspec[j - 1]) && (smspec[j] > smspec[j + 1]) && (npk < 200)) { + if ((smspec[j] > smspec[j - 1]) && + (smspec[j] > smspec[j + 1]) && + (npk < 200)) { freq0[npk] = (j - 205) * (DF / 2.0); snr0[npk] = 10.0 * log10f(smspec[j]) - snr_scaling_factor; npk++; } } - /* Compute corrected fmin, fmax, accounting for dial frequency error - float dialfreq_error = 0.0; // dialfreq_error is in units of Hz - fmin += dialfreq_error; - fmax += dialfreq_error; - */ - // Don't waste time on signals outside of the range [fmin,fmax]. - i = 0; + int i = 0; for (int j = 0; j < npk; j++) { if (freq0[j] >= fmin && freq0[j] <= fmax) { freq0[i] = freq0[j]; @@ -782,7 +612,7 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, // bubble sort on snr, bringing freq along for the ride float tmp; for (int pass = 1; pass <= npk - 1; pass++) { - for (k = 0; k < npk - pass; k++) { + for (int k = 0; k < npk - pass; k++) { if (snr0[k] < snr0[k + 1]) { tmp = snr0[k]; snr0[k] = snr0[k + 1]; @@ -815,10 +645,10 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, for (int idrift = -maxdrift; idrift <= maxdrift; idrift++) { // Drift search float ss = 0.0; float pow = 0.0; - for (k = 0; k < NSYM; k++) { // Sum over symbols + for (int k = 0; k < NSYM; k++) { // Sum over symbols int ifd = ifr + ((float)k - (float)NBITS) / (float)NBITS * ((float)idrift) / DF; int kindex = k0 + 2 * k; - if (kindex < nffts) { + if (kindex < blocks) { float p0 = sqrtf(ps[ifd - 3][kindex]); float p1 = sqrtf(ps[ifd - 1][kindex]); float p2 = sqrtf(ps[ifd + 1][kindex]); @@ -871,23 +701,20 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, shift1 = shift0[j]; sync1 = sync0[j]; - // Fine search for best sync lag (mode 0) - float fstep = 0.0; - int32_t lagmin = shift1 - 144; - int32_t lagmax = shift1 + 144; - int32_t lagstep = 8; + // Coarse-grid search for best sync lag (mode 0) + fstep = 0.0; ifmin = 0; ifmax = 0; - + lagmin = shift1 - 128; + lagmax = shift1 + 128; + lagstep = 8; if (options.quickmode) lagstep = 16; - sync_and_demodulate(idat, qdat, npoints, symbols, &freq1, ifmin, ifmax, fstep, &shift1, lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0); - // Fine search for frequency peak (mode 1) + // Coarse-grid search for frequency peak (mode 1) fstep = 0.1; - ifmin = 0; ifmin = -2; ifmax = 2; sync_and_demodulate(idat, qdat, npoints, symbols, &freq1, ifmin, ifmax, fstep, &shift1, @@ -899,26 +726,24 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, worth_a_try = 0; } - int32_t idt = 0, ii = 0, jiggered_shift; - float y, sq, rms; - int32_t not_decoded = 1; - + int idt = 0, ii = 0; + int 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; + int jiggered_shift = shift1 + ii; // Use mode 2 to get soft-decision symbols sync_and_demodulate(idat, qdat, npoints, symbols, &freq1, ifmin, ifmax, fstep, &jiggered_shift, lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 2); - sq = 0.0; + float sq = 0.0; for (i = 0; i < NSYM; i++) { - y = (float)symbols[i] - 128.0; + float y = (float)symbols[i] - 128.0; sq += y * y; } - rms = sqrtf(sq / (float)NSYM); + float rms = sqrtf(sq / (float)NSYM); if ((sync1 > minsync2) && (rms > minrms)) { deinterleave(symbols); @@ -987,8 +812,8 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, // sort the result struct decoder_results temp; - for (j = 1; j <= uniques - 1; j++) { - for (k = 0; k < uniques - j; k++) { + for (int j = 1; j <= uniques - 1; j++) { + for (int k = 0; k < uniques - j; k++) { if (decodes[k].snr < decodes[k + 1].snr) { temp = decodes[k]; decodes[k] = decodes[k + 1]; @@ -1013,12 +838,13 @@ int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, if (options.usehashtable) { fhash = fopen("hashtable.txt", "w"); - for (i = 0; i < 32768; i++) { + for (int i = 0; i < 32768; i++) { if (strncmp(hashtab + i * 13, "\0", 1) != 0) { - fprintf(fhash, "%5d %s\n", i, hashtab + i * 13); + fprintf(fhash, "%5d %s %s\n", i, hashtab + i * 13, loctab + i * 5); } } fclose(fhash); } + return 0; } diff --git a/wsprd/wsprd.h b/wsprd/wsprd.h index ec3b7c5..0a44d89 100644 --- a/wsprd/wsprd.h +++ b/wsprd/wsprd.h @@ -66,44 +66,14 @@ struct decoder_results { uint32_t cycles; }; -// 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; -// int blocksize; -// unsigned int metric; -// int nhardmin; -// int ipass; -// int decodetype; -// }; - 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); -void noncoherent_sequence_detection(float *id, float *qd, long np, - unsigned char *symbols, float *f1, int *shift1, - float *drift1, int symfac, int *nblocksize, int *bitmetric); void subtract_signal(float *id, float *qd, long np, float f0, int shift0, float drift0, unsigned char *channel_symbols); void subtract_signal2(float *id, float *qd, long np, float f0, int shift0, float drift0, unsigned char *channel_symbols); -unsigned int count_hard_errors(unsigned char *symbols, unsigned char *channel_symbols); int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints, struct decoder_options options, struct decoder_results *decodes, int32_t *n_results); - - - - - - - -