kopia lustrzana https://github.com/Guenael/rtlsdr-wsprd
WIP chore: wsprd lib update part2
Guenael
rodzic
b248ceb121
commit
1a167c87bf
4
Makefile
4
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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
602
wsprd/wsprd.c
602
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue