mirror
/
rtlsdr-wsprd
kopia lustrzana https://github.com/Guenael/rtlsdr-wsprd
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

memory leak fix & float approx fix

pull/5/head
Guenael 2016-06-24 20:58:12 -04:00
rodzic 6b42675dd5
commit 81edbc091c
4 zmienionych plików z 169 dodań i 202 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

2
metric_tables.c → metric_tables.h
Wyświetl plik

@ -8,7 +8,7 @@
#pragma GCC diagnostic ignored "-Wmissing-braces" #pragma GCC diagnostic ignored "-Wmissing-braces"
//float symbol_scale[4]={42.6, 53.3, 72.7, 100.2}; //float symbol_scale[4]={42.6, 53.3, 72.7, 100.2};
float metric_tables[4][256]= { static const float metric_tables[4][256]= {
0.9782, 0.9695, 0.9689, 0.9669, 0.9666, 0.9653, 0.9638, 0.9618, 0.9599, 0.9601, 0.9782, 0.9695, 0.9689, 0.9669, 0.9666, 0.9653, 0.9638, 0.9618, 0.9599, 0.9601,
0.9592, 0.9570, 0.9556, 0.9540, 0.9525, 0.9527, 0.9486, 0.9477, 0.9450, 0.9436, 0.9592, 0.9570, 0.9556, 0.9540, 0.9525, 0.9527, 0.9486, 0.9477, 0.9450, 0.9436,
0.9424, 0.9400, 0.9381, 0.9360, 0.9340, 0.9316, 0.9301, 0.9272, 0.9254, 0.9224, 0.9424, 0.9400, 0.9381, 0.9360, 0.9340, 0.9316, 0.9301, 0.9272, 0.9254, 0.9224,

23
rtlsdr_wsprd.c
Wyświetl plik

@ -40,7 +40,7 @@
#include <sys/time.h> #include <sys/time.h>
#include <pthread.h> #include <pthread.h>
#include <curl/curl.h> #include <curl/curl.h>
#include <pthread.h>
#include <rtl-sdr.h> #include <rtl-sdr.h>
#include "rtlsdr_wsprd.h" #include "rtlsdr_wsprd.h"
@ -49,7 +49,6 @@
/* TODO /* TODO
- multi device selection option - multi device selection option
- multispot report in one post - multispot report in one post
- type fix (uint32_t etc..)
- verbose option - verbose option
*/ */
@ -73,6 +72,7 @@ static rtlsdr_dev_t *rtl_device = NULL;
/* Thread stuff for separate decoding */ /* Thread stuff for separate decoding */
struct decoder_state { struct decoder_state {
pthread_t thread; pthread_t thread;
pthread_attr_t tattr;
pthread_rwlock_t rw; pthread_rwlock_t rw;
pthread_cond_t ready_cond; pthread_cond_t ready_cond;
@ -141,7 +141,7 @@ static void rtlsdr_callback(unsigned char *samples, uint32_t samples_count, void
out_Q = in_Q * cos(x) + in_I * sin(x) out_Q = in_Q * cos(x) + in_I * sin(x)
(Weaver technique, keep the upper band, IQ inverted on RTL devices) (Weaver technique, keep the upper band, IQ inverted on RTL devices)
*/ */
int16_t tmp; int8_t tmp;
for (uint32_t i=0; i<sigLenght; i+=8) { for (uint32_t i=0; i<sigLenght; i+=8) {
tmp = -sigIn[i+3]; tmp = -sigIn[i+3];
sigIn[i+3] = sigIn[i+2]; sigIn[i+3] = sigIn[i+2];
@ -425,7 +425,7 @@ void usage(void) {
"\t-Q quick mode, doesn't dig deep for weak signals\n" "\t-Q quick mode, doesn't dig deep for weak signals\n"
"\t-S single pass mode, no subtraction (same as original wsprd)\n" "\t-S single pass mode, no subtraction (same as original wsprd)\n"
"Example:\n" "Example:\n"
"\trtlsdr_wsprd -f 144.489M -c A1XYZ -l AB12cd -g 29 -o -4080\n"); "\trtlsdr_wsprd -f 144.489M -c A1XYZ -l AB12cd -g 29 -o -4200\n");
exit(1); exit(1);
} }
@ -480,7 +480,7 @@ int main(int argc, char** argv) {
case 'p': case 'p':
rx_options.ppm = atoi(optarg); rx_options.ppm = atoi(optarg);
break; break;
case 'u': // Frequency case 'u': // Upconverter frequency
rx_options.upconverter = (uint32_t)atofs(optarg); rx_options.upconverter = (uint32_t)atofs(optarg);
break; break;
case 'H': // Decoder option, use a hastable case 'H': // Decoder option, use a hastable
@ -616,7 +616,7 @@ int main(int argc, char** argv) {
time_t rawtime; time_t rawtime;
time ( &rawtime ); time ( &rawtime );
struct tm *gtm = gmtime(&rawtime); struct tm *gtm = gmtime(&rawtime);
printf("\nStarting rtlsdr-wsprd (%04d-%02d-%02d, %02d:%02dz) -- Version 0.1\n", printf("\nStarting rtlsdr-wsprd (%04d-%02d-%02d, %02d:%02dz) -- Version 0.2\n",
gtm->tm_year + 1900, gtm->tm_mon + 1, gtm->tm_mday, gtm->tm_hour, gtm->tm_min); gtm->tm_year + 1900, gtm->tm_mon + 1, gtm->tm_mday, gtm->tm_hour, gtm->tm_min);
printf(" Callsign : %s\n", dec_options.rcall); printf(" Callsign : %s\n", dec_options.rcall);
printf(" Locator : %s\n", dec_options.rloc); printf(" Locator : %s\n", dec_options.rloc);
@ -637,6 +637,13 @@ int main(int argc, char** argv) {
uint32_t uwait = 120000000 - usec; uint32_t uwait = 120000000 - usec;
printf("Wait for time sync (start in %d sec)\n\n", uwait/1000000); printf("Wait for time sync (start in %d sec)\n\n", uwait/1000000);
/* Prepare a low priority param for the decoder thread */
struct sched_param p;
pthread_attr_init(&dec.tattr);
pthread_attr_setschedpolicy(&dec.tattr, SCHED_RR);
pthread_attr_getschedparam(&dec.tattr, &p);
p.sched_priority = 40;
pthread_attr_setschedparam(&dec.tattr, &p);
/* Create a thread and stuff for separate decoding /* Create a thread and stuff for separate decoding
Info : https://computing.llnl.gov/tutorials/pthreads/ Info : https://computing.llnl.gov/tutorials/pthreads/
@ -644,9 +651,8 @@ int main(int argc, char** argv) {
pthread_rwlock_init(&dec.rw, NULL); pthread_rwlock_init(&dec.rw, NULL);
pthread_cond_init(&dec.ready_cond, NULL); pthread_cond_init(&dec.ready_cond, NULL);
pthread_mutex_init(&dec.ready_mutex, NULL); pthread_mutex_init(&dec.ready_mutex, NULL);
pthread_create(&dec.thread, NULL, wsprDecoder, NULL);
pthread_create(&dongle.thread, NULL, rtlsdr_rx, NULL); pthread_create(&dongle.thread, NULL, rtlsdr_rx, NULL);
pthread_create(&dec.thread, &dec.tattr, wsprDecoder, NULL);
/* Main loop : Wait, read, decode */ /* Main loop : Wait, read, decode */
while (!rx_state.exit_flag) { while (!rx_state.exit_flag) {
@ -669,7 +675,6 @@ int main(int argc, char** argv) {
while( (rx_state.exit_flag == false) && while( (rx_state.exit_flag == false) &&
(rx_state.iqIndex < (SIGNAL_LENGHT * SIGNAL_SAMPLE_RATE) ) ) { (rx_state.iqIndex < (SIGNAL_LENGHT * SIGNAL_SAMPLE_RATE) ) ) {
//sleep(1);
usleep(250000); usleep(250000);
} }
} }

329
wsprd.c
Wyświetl plik

@ -11,6 +11,10 @@
Copyright 2014-2015, Steven Franke, K9AN Copyright 2014-2015, Steven Franke, K9AN
Minor modifications
Copyright 2016, Guenael Jouchet, VA2GKA
License: GNU GPL v3 License: GNU GPL v3
This program is free software: you can redistribute it and/or modify This program is free software: you can redistribute it and/or modify
@ -41,19 +45,24 @@
#include "nhash.h" #include "nhash.h"
#include "wsprd_utils.h" #include "wsprd_utils.h"
#include "wsprsim_utils.h" #include "wsprsim_utils.h"
#include "metric_tables.h"
#define DT 1.0/375.0
#define DF 375.0/256.0 #define DF 375.0/256.0
#define TWOPIDT 2*M_PI*DT #define DT 1.0/375.0
#define TWOPIDT 2 * M_PI * DT
#define NIQ 45000
#define NBITS 81
#define NSYM 162
#define NSPERSYM 256
#define NFILT 256
#define NSIG NSYM * NSPERSYM
#define max(x,y) ((x) > (y) ? (x) : (y)) /* Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT, FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE */
// Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT,
// FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
#define PATIENCE FFTW_ESTIMATE #define PATIENCE FFTW_ESTIMATE
fftwf_plan PLAN1,PLAN2,PLAN3;
unsigned char pr3[162]= { uint8_t pr3[NSYM]= {
1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0, 1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0,
0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1, 0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1,
0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1, 0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1,
@ -65,16 +74,17 @@ unsigned char pr3[162]= {
0,0 0,0
}; };
unsigned long nr; fftwf_plan PLAN1,
PLAN2,
int printdata=0; PLAN3;
int32_t printdata=0;
//*************************************************************************** //***************************************************************************
void sync_and_demodulate(float *id, float *qd, long np, void sync_and_demodulate(float *id, float *qd, long np,
unsigned char *symbols, float *f1, float fstep, uint8_t *symbols, float *f1, float fstep,
int *shift1, int lagmin, int lagmax, int lagstep, int32_t *shift1, int32_t lagmin, int32_t lagmax, int32_t lagstep,
float *drift1, int symfac, float *sync, int mode) { float *drift1, int32_t symfac, float *sync, int32_t mode) {
/*********************************************************************** /***********************************************************************
* mode = 0: no frequency or drift search. find best time lag. * * mode = 0: no frequency or drift search. find best time lag. *
* 1: no time lag or drift search. find best frequency. * * 1: no time lag or drift search. find best frequency. *
@ -84,18 +94,24 @@ void sync_and_demodulate(float *id, float *qd, long np,
float fbest=0.0; float fbest=0.0;
float f0=0.0,fp,ss; float f0=0.0,fp,ss;
int lag; int32_t lag;
static float fplast=-10000.0; static float fplast=-10000.0;
float i0[162],q0[162],i1[162],q1[162],i2[162],q2[162],i3[162],q3[162]; float i0[NSYM],q0[NSYM],
i1[NSYM],q1[NSYM],
i2[NSYM],q2[NSYM],
i3[NSYM],q3[NSYM];
float p0,p1,p2,p3,cmet,totp,syncmax,fac; float p0,p1,p2,p3,cmet,totp,syncmax,fac;
float c0[256],s0[256],c1[256],s1[256],c2[256],s2[256],c3[256],s3[256]; float c0[NSPERSYM],s0[NSPERSYM],
c1[NSPERSYM],s1[NSPERSYM],
c2[NSPERSYM],s2[NSPERSYM],
c3[NSPERSYM],s3[NSPERSYM];
float dphi0, cdphi0, sdphi0, float dphi0, cdphi0, sdphi0,
dphi1, cdphi1, sdphi1, dphi1, cdphi1, sdphi1,
dphi2, cdphi2, sdphi2, dphi2, cdphi2, sdphi2,
dphi3, cdphi3, sdphi3; dphi3, cdphi3, sdphi3;
float fsum=0.0, f2sum=0.0, fsymb[162]; float fsum=0.0, f2sum=0.0, fsymb[NSYM];
int best_shift = 0; int32_t best_shift = 0;
int ifmin=0, ifmax=0; int32_t ifmin=0, ifmax=0;
syncmax=-1e30; syncmax=-1e30;
if( mode == 0 ) { if( mode == 0 ) {
@ -119,12 +135,12 @@ void sync_and_demodulate(float *id, float *qd, long np,
f0=*f1; f0=*f1;
} }
for(int ifreq=ifmin; ifreq<=ifmax; ifreq++) { for(int32_t ifreq=ifmin; ifreq<=ifmax; ifreq++) {
f0=*f1+ifreq*fstep; f0=*f1+ifreq*fstep;
for(lag=lagmin; lag<=lagmax; lag=lag+lagstep) { for(lag=lagmin; lag<=lagmax; lag=lag+lagstep) {
ss=0.0; ss=0.0;
totp=0.0; totp=0.0;
for (int i=0; i<162; i++) { for (int32_t i=0; i<NSYM; i++) {
fp = f0 + ((float)*drift1/2.0)*((float)i-81.0)/81.0; fp = f0 + ((float)*drift1/2.0)*((float)i-81.0)/81.0;
if( i==0 || (fp != fplast) ) { // only calculate sin/cos if necessary if( i==0 || (fp != fplast) ) { // only calculate sin/cos if necessary
dphi0=TWOPIDT*(fp-1.5*DF); dphi0=TWOPIDT*(fp-1.5*DF);
@ -152,7 +168,7 @@ void sync_and_demodulate(float *id, float *qd, long np,
c3[0]=1; c3[0]=1;
s3[0]=0; s3[0]=0;
for (int j=1; j<256; j++) { for (int32_t j=1; j<NSPERSYM; j++) {
c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0; c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0;
s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0; s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0;
c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1; c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1;
@ -174,8 +190,8 @@ void sync_and_demodulate(float *id, float *qd, long np,
i3[i]=0.0; i3[i]=0.0;
q3[i]=0.0; q3[i]=0.0;
for (int j=0; j<256; j++) { for (int32_t j=0; j<NSPERSYM; j++) {
int k=lag+i*256+j; int32_t k=lag+i*NSPERSYM+j;
if( (k>0) & (k<np) ) { if( (k>0) & (k<np) ) {
i0[i]=i0[i] + id[k]*c0[j] + qd[k]*s0[j]; i0[i]=i0[i] + id[k]*c0[j] + qd[k]*s0[j];
q0[i]=q0[i] - id[k]*s0[j] + qd[k]*c0[j]; q0[i]=q0[i] - id[k]*s0[j] + qd[k]*c0[j];
@ -226,12 +242,12 @@ void sync_and_demodulate(float *id, float *qd, long np,
if( mode == 2 ) { if( mode == 2 ) {
*sync=syncmax; *sync=syncmax;
for (int i=0; i<162; i++) { //Normalize the soft symbols for (int32_t i=0; i<NSYM; i++) { //Normalize the soft symbols
fsum=fsum+fsymb[i]/162.0; fsum=fsum+fsymb[i]/NSYM;
f2sum=f2sum+fsymb[i]*fsymb[i]/162.0; f2sum=f2sum+fsymb[i]*fsymb[i]/NSYM;
} }
fac=sqrtf(f2sum-fsum*fsum); fac=sqrtf(f2sum-fsum*fsum);
for (int i=0; i<162; i++) { for (int32_t i=0; i<NSYM; i++) {
fsymb[i]=symfac*fsymb[i]/fac; fsymb[i]=symfac*fsymb[i]/fac;
if( fsymb[i] > 127) fsymb[i]=127.0; if( fsymb[i] > 127) fsymb[i]=127.0;
if( fsymb[i] < -128 ) fsymb[i]=-128.0; if( fsymb[i] < -128 ) fsymb[i]=-128.0;
@ -247,13 +263,13 @@ void sync_and_demodulate(float *id, float *qd, long np,
symbol-by-symbol signal subtraction symbol-by-symbol signal subtraction
****************************************************************************/ ****************************************************************************/
void subtract_signal(float *id, float *qd, long np, void subtract_signal(float *id, float *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols) { float f0, int32_t shift0, float drift0, uint8_t* channel_symbols) {
float i0,q0; float i0,q0;
float c0[256],s0[256]; float c0[NSPERSYM],s0[NSPERSYM];
float dphi, cdphi, sdphi; float dphi, cdphi, sdphi;
for (int i=0; i<162; i++) { for (int32_t i=0; i<NSYM; i++) {
float fp = f0 + ((float)drift0/2.0)*((float)i-81.0)/81.0; float fp = f0 + ((float)drift0/2.0)*((float)i-81.0)/81.0;
dphi=TWOPIDT*(fp+((float)channel_symbols[i]-1.5)*DF); dphi=TWOPIDT*(fp+((float)channel_symbols[i]-1.5)*DF);
@ -263,7 +279,7 @@ void subtract_signal(float *id, float *qd, long np,
c0[0]=1; c0[0]=1;
s0[0]=0; s0[0]=0;
for (int j=1; j<256; j++) { for (int32_t j=1; j<NSPERSYM; j++) {
c0[j]=c0[j-1]*cdphi - s0[j-1]*sdphi; c0[j]=c0[j-1]*cdphi - s0[j-1]*sdphi;
s0[j]=c0[j-1]*sdphi + s0[j-1]*cdphi; s0[j]=c0[j-1]*sdphi + s0[j-1]*cdphi;
} }
@ -271,8 +287,8 @@ void subtract_signal(float *id, float *qd, long np,
i0=0.0; i0=0.0;
q0=0.0; q0=0.0;
for (int j=0; j<256; j++) { for (int32_t j=0; j<NSPERSYM; j++) {
int k=shift0+i*256+j; int32_t k=shift0+i*NSPERSYM+j;
if( (k>0) & (k<np) ) { if( (k>0) & (k<np) ) {
i0=i0 + id[k]*c0[j] + qd[k]*s0[j]; i0=i0 + id[k]*c0[j] + qd[k]*s0[j];
q0=q0 - id[k]*s0[j] + qd[k]*c0[j]; q0=q0 - id[k]*s0[j] + qd[k]*c0[j];
@ -282,11 +298,11 @@ void subtract_signal(float *id, float *qd, long np,
// subtract the signal here. // subtract the signal here.
i0=i0/256.0; //will be wrong for partial symbols at the edges... i0=i0/NSPERSYM; //will be wrong for partial symbols at the edges...
q0=q0/256.0; q0=q0/NSPERSYM;
for (int j=0; j<256; j++) { for (int32_t j=0; j<NSPERSYM; j++) {
int k=shift0+i*256+j; int32_t k=shift0+i*NSPERSYM+j;
if( (k>0) & (k<np) ) { if( (k>0) & (k<np) ) {
id[k]=id[k]- (i0*c0[j] - q0*s0[j]); id[k]=id[k]- (i0*c0[j] - q0*s0[j]);
qd[k]=qd[k]- (q0*c0[j] + i0*s0[j]); qd[k]=qd[k]- (q0*c0[j] + i0*s0[j]);
@ -301,29 +317,12 @@ void subtract_signal(float *id, float *qd, long np,
Fully coherent signal subtraction Fully coherent signal subtraction
*******************************************************************************/ *******************************************************************************/
void subtract_signal2(float *id, float *qd, long np, void subtract_signal2(float *id, float *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols) { float f0, int32_t shift0, float drift0, uint8_t* channel_symbols) {
float phi=0, dphi, cs; float phi=0, dphi, cs;
int nsym=162, nspersym=256, nfilt=256; //nfilt must be even number. float refi[NIQ]= {0}, refq[NIQ]= {0},
int nsig=nsym*nspersym; ci[NIQ]= {0}, cq[NIQ]= {0},
int nc2=45000; cfi[NIQ]= {0}, cfq[NIQ]= {0};
float *refi, *refq, *ci, *cq, *cfi, *cfq;
refi=malloc(sizeof(float)*nc2);
refq=malloc(sizeof(float)*nc2);
ci=malloc(sizeof(float)*nc2);
cq=malloc(sizeof(float)*nc2);
cfi=malloc(sizeof(float)*nc2);
cfq=malloc(sizeof(float)*nc2);
memset(refi,0,sizeof(float)*nc2);
memset(refq,0,sizeof(float)*nc2);
memset(ci,0,sizeof(float)*nc2);
memset(cq,0,sizeof(float)*nc2);
memset(cfi,0,sizeof(float)*nc2);
memset(cfq,0,sizeof(float)*nc2);
/****************************************************************************** /******************************************************************************
Measured signal: s(t)=a(t)*exp( j*theta(t) ) Measured signal: s(t)=a(t)*exp( j*theta(t) )
@ -335,16 +334,16 @@ void subtract_signal2(float *id, float *qd, long np,
// create reference wspr signal vector, centered on f0. // create reference wspr signal vector, centered on f0.
// //
for (int i=0; i<nsym; i++) { for (int32_t i=0; i<NSYM; i++) {
cs=(float)channel_symbols[i]; cs=(float)channel_symbols[i];
dphi=TWOPIDT* ( f0 + dphi=TWOPIDT * ( f0 +
((float)drift0/2.0)*((float)i-(float)nsym/2.0)/((float)nsym/2.0) + ((float)drift0/2.0)*((float)i-(float)NSYM/2.0)/((float)NSYM/2.0) +
(cs-1.5)*DF ); (cs-1.5)*DF );
for (int j=0; j<nspersym; j++ ) { for (int32_t j=0; j<NSPERSYM; j++ ) {
int ii=nspersym*i+j; int32_t ii=NSPERSYM*i+j;
refi[ii]=refi[ii]+cosf(phi); //cannot precompute sin/cos because dphi is changing refi[ii]=refi[ii]+cosf(phi); //cannot precompute sin/cos because dphi is changing
refq[ii]=refq[ii]+sinf(phi); refq[ii]=refq[ii]+sinf(phi);
phi=phi+dphi; phi=phi+dphi;
@ -356,133 +355,111 @@ void subtract_signal2(float *id, float *qd, long np,
// beginning of first symbol in received data is at shift0. // beginning of first symbol in received data is at shift0.
// filter transient lasts nfilt samples // filter transient lasts nfilt samples
// leave nfilt zeros as a pad at the beginning of the unfiltered reference signal // leave nfilt zeros as a pad at the beginning of the unfiltered reference signal
for (int i=0; i<nsym*nspersym; i++) { for (int32_t i=0; i<NSYM*NSPERSYM; i++) {
int k=shift0+i; int32_t k=shift0+i;
if( (k>0) & (k<np) ) { if( (k>0) & (k<np) ) {
ci[i+nfilt] = id[k]*refi[i] + qd[k]*refq[i]; ci[i+NFILT] = id[k]*refi[i] + qd[k]*refq[i];
cq[i+nfilt] = qd[k]*refi[i] - id[k]*refq[i]; cq[i+NFILT] = qd[k]*refi[i] - id[k]*refq[i];
} }
} }
//quick and dirty filter - may want to do better //quick and dirty filter - may want to do better
float w[nfilt], norm=0, partialsum[nfilt]; float w[NFILT]= {0}, norm=0, partialsum[NFILT]= {0};
memset(partialsum,0,sizeof(float)*nfilt); for (int32_t i=0; i<NFILT; i++) {
for (int i=0; i<nfilt; i++) { w[i]=sinf(M_PI*(float)i/(float)(NFILT-1));
w[i]=sinf(M_PI*(float)i/(float)(nfilt-1));
norm=norm+w[i]; norm=norm+w[i];
} }
for (int i=0; i<nfilt; i++) { for (int32_t i=0; i<NFILT; i++) {
w[i]=w[i]/norm; w[i]=w[i]/norm;
} }
for (int i=1; i<nfilt; i++) { for (int32_t i=1; i<NFILT; i++) {
partialsum[i]=partialsum[i-1]+w[i]; partialsum[i]=partialsum[i-1]+w[i];
} }
// LPF // LPF
for (int i=nfilt/2; i<45000-nfilt/2; i++) { for (int32_t i=NFILT/2; i<NIQ-NFILT/2; i++) {
cfi[i]=0.0; cfi[i]=0.0;
cfq[i]=0.0; cfq[i]=0.0;
for (int j=0; j<nfilt; j++) { for (int32_t j=0; j<NFILT; j++) {
cfi[i]=cfi[i]+w[j]*ci[i-nfilt/2+j]; cfi[i]=cfi[i]+w[j]*ci[i-NFILT/2+j];
cfq[i]=cfq[i]+w[j]*cq[i-nfilt/2+j]; cfq[i]=cfq[i]+w[j]*cq[i-NFILT/2+j];
} }
} }
// subtract c(t)*r(t) here // subtract c(t)*r(t) here
// (ci+j*cq)(refi+j*refq)=(ci*refi-cq*refq)+j(ci*refq)+cq*refi) // (ci+j*cq)(refi+j*refq)=(ci*refi-cq*refq)+j(ci*refq)+cq*refi)
// beginning of first symbol in reference signal is at i=nfilt // beginning of first symbol in reference signal is at i=NFILT
// beginning of first symbol in received data is at shift0. // beginning of first symbol in received data is at shift0.
for (int i=0; i<nsig; i++) { for (int32_t i=0; i<NSIG; i++) {
if( i<nfilt/2 ) { // take care of the end effect (LPF step response) here if( i<NFILT/2 ) { // take care of the end effect (LPF step response) here
norm=partialsum[nfilt/2+i]; norm=partialsum[NFILT/2+i];
} else if( i>(nsig-1-nfilt/2) ) { } else if( i>(NSIG-1-NFILT/2) ) {
norm=partialsum[nfilt/2+nsig-1-i]; norm=partialsum[NFILT/2+NSIG-1-i];
} else { } else {
norm=1.0; norm=1.0;
} }
int k=shift0+i; int32_t k=shift0+i;
int j=i+nfilt; int32_t j=i+NFILT;
if( (k>0) & (k<np) ) { if( (k>0) & (k<np) ) {
id[k]=id[k] - (cfi[j]*refi[i]-cfq[j]*refq[i])/norm; id[k]=id[k] - (cfi[j]*refi[i]-cfq[j]*refq[i])/norm;
qd[k]=qd[k] - (cfi[j]*refq[i]+cfq[j]*refi[i])/norm; qd[k]=qd[k] - (cfi[j]*refq[i]+cfq[j]*refi[i])/norm;
} }
} }
free(refi);
free(refq);
free(ci);
free(cq);
free(cfi);
free(cfq);
return; return;
} }
//*************************************************************************** //***************************************************************************
int wspr_decode(float *idat, float *qdat, unsigned int npoints, int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints,
struct decoder_options options, struct decoder_results *decodes, int *n_results) { struct decoder_options options, struct decoder_results *decodes,
int32_t *n_results) {
int i,j,k; int32_t i,j,k;
unsigned char *symbols, *decdata; uint32_t metric, maxcycles, cycles, maxnp;
signed char message[]= {-9,13,-35,123,57,-39,64,0,0,0,0}; uint8_t symbols[NBITS*2]= {0};
char *callsign, *call_loc_pow, *call, *loc, *pwr; uint8_t decdata[(NBITS+7)/8]= {0};
char *data_dir=NULL; int8_t message[]= {-9,13,-35,123,57,-39,64,0,0,0,0};
char wisdom_fname[200],all_fname[200],spots_fname[200];
char timer_fname[200],hash_fname[200]; char callsign[13]= {0};
int delta,verbose=0; char call_loc_pow[23]= {0};
int writenoise=0,wspr_type=2, ipass; char call[13]= {0};
int shift1, lagmin, lagmax, lagstep, worth_a_try, not_decoded; char loc[7]= {0};
unsigned int nbits=81; char pwr[3]= {0};
unsigned int metric, maxcycles, cycles, maxnp;
int32_t delta,verbose=0;
int32_t writenoise=0,wspr_type=2, ipass;
int32_t shift1, lagmin, lagmax, lagstep, worth_a_try, not_decoded;
float freq0[200],snr0[200],drift0[200],sync0[200]; float freq0[200],snr0[200],drift0[200],sync0[200];
int shift0[200]; int32_t shift0[200];
float dt_print; float dt_print;
float freq_print; double freq_print;
float dialfreq= (float)options.freq / 1e6; // check double dialfreq= (double)options.freq / 1e6; // check
float dialfreq_error=0.0; float dialfreq_error=0.0;
float f1, fstep, sync1=0.0, drift1; float f1, fstep, sync1=0.0, drift1;
int noprint=0; int32_t noprint=0;
int uniques=0; int32_t uniques=0;
float fmin=-110.0; float fmin=-110.0;
float fmax=110.0; float fmax=110.0;
char hashtab[32768*13]= {0};
int32_t nh;
char *hashtab; float allfreqs[100]= {0};
hashtab=malloc(sizeof(char)*32768*13); char allcalls[100][13]= {0};
memset(hashtab,0,sizeof(char)*32768*13);
int nh;
symbols=malloc(sizeof(char)*nbits*2);
decdata=malloc((nbits+7)/8);
callsign=malloc(sizeof(char)*13);
call_loc_pow=malloc(sizeof(char)*23);
call=malloc(sizeof(char)*13);
loc=malloc(sizeof(char)*7);
pwr=malloc(sizeof(char)*3);
float allfreqs[100];
char allcalls[100][13];
memset(allfreqs,0,sizeof(float)*100);
memset(allcalls,0,sizeof(char)*100*13);
// Parameters used for performance-tuning: // Parameters used for performance-tuning:
maxcycles=10000; //Fano timeout limit maxcycles=10000; //Fano timeout limit
double minsync1=0.10; //First sync limit double minsync1=0.10; //First sync limit
double minsync2=0.12; //Second sync limit double minsync2=0.12; //Second sync limit
int iifac=3; //Step size in final DT peakup int32_t iifac=3; //Step size in final DT peakup
int symfac=50; //Soft-symbol normalizing factor int32_t symfac=50; //Soft-symbol normalizing factor
int maxdrift=4; //Maximum (+/-) drift int32_t maxdrift=4; //Maximum (+/-) drift
double minrms=52.0 * (symfac/64.0); //Final test for plausible decoding double minrms=52.0 * (symfac/64.0); //Final test for plausible decoding
delta=60; //Fano threshold step delta=60; //Fano threshold step
fftwf_complex *fftin, *fftout; fftwf_complex *fftin, *fftout;
#include "./metric_tables.c" int32_t mettab[2][256];
int mettab[2][256];
float bias=0.42; float bias=0.42;
// setup metric table // setup metric table
@ -492,30 +469,14 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
} }
FILE *fp_fftw_wisdom_file, *fhash; FILE *fp_fftw_wisdom_file, *fhash;
strcpy(wisdom_fname,".");
strcpy(all_fname,"."); if ((fp_fftw_wisdom_file = fopen("wspr_wisdom.dat", "r"))) { //Open FFTW wisdom
strcpy(spots_fname,".");
strcpy(timer_fname,".");
strcpy(hash_fname,".");
if(data_dir != NULL) {
strcpy(wisdom_fname,data_dir);
strcpy(all_fname,data_dir);
strcpy(spots_fname,data_dir);
strcpy(timer_fname,data_dir);
strcpy(hash_fname,data_dir);
}
strncat(wisdom_fname,"/wspr_wisdom.dat",20);
strncat(all_fname,"/ALL_WSPR.TXT",20);
strncat(spots_fname,"/wspr_spots.txt",20);
strncat(timer_fname,"/wspr_timer.out",20);
strncat(hash_fname,"/hashtable.txt",20);
if ((fp_fftw_wisdom_file = fopen(wisdom_fname, "r"))) { //Open FFTW wisdom
fftwf_import_wisdom_from_file(fp_fftw_wisdom_file); fftwf_import_wisdom_from_file(fp_fftw_wisdom_file);
fclose(fp_fftw_wisdom_file); fclose(fp_fftw_wisdom_file);
} }
// Do windowed ffts over 2 symbols, stepped by half symbols // Do windowed ffts over 2 symbols, stepped by half symbols
int nffts=4*floor(npoints/512)-1; int32_t nffts=4*floor(npoints/512)-1;
fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512); fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512);
fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512); fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512);
PLAN3 = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE); PLAN3 = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE);
@ -528,13 +489,13 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
if( options.usehashtable ) { if( options.usehashtable ) {
char line[80], hcall[12]; char line[80], hcall[12];
if( (fhash=fopen(hash_fname,"r+")) ) { if( (fhash=fopen("hashtable.txt","r+")) ) {
while (fgets(line, sizeof(line), fhash) != NULL) { while (fgets(line, sizeof(line), fhash) != NULL) {
sscanf(line,"%d %s",&nh,hcall); sscanf(line,"%d %s",&nh,hcall);
strcpy(hashtab+nh*13,hcall); strcpy(hashtab+nh*13,hcall);
} }
} else { } else {
fhash=fopen(hash_fname,"w+"); fhash=fopen("hashtable.txt","w+");
} }
fclose(fhash); fclose(fhash);
} }
@ -565,7 +526,6 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
// Compute average spectrum // Compute average spectrum
float psavg[512]= {0}; float psavg[512]= {0};
//memset(psavg,0.0, sizeof(float)*512);
for (i=0; i<nffts; i++) { for (i=0; i<nffts; i++) {
for (j=0; j<512; j++) { for (j=0; j<512; j++) {
psavg[j]=psavg[j]+ps[j][i]; psavg[j]=psavg[j]+ps[j][i];
@ -573,7 +533,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
} }
// Smooth with 7-point window and limit spectrum to +/-150 Hz // Smooth with 7-point window and limit spectrum to +/-150 Hz
int window[7]= {1,1,1,1,1,1,1}; int32_t window[7]= {1,1,1,1,1,1,1};
float smspec[411]; float smspec[411];
for (i=0; i<411; i++) { for (i=0; i<411; i++) {
smspec[i]=0.0; smspec[i]=0.0;
@ -620,7 +580,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
sync0[i]=0.0; sync0[i]=0.0;
} }
int npk=0; int32_t npk=0;
for(j=1; j<410; j++) { for(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); freq0[npk]=(j-205)*(DF/2.0);
@ -645,7 +605,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
npk=i; npk=i;
// bubble sort on snr, bringing freq along for the ride // bubble sort on snr, bringing freq along for the ride
int pass; int32_t pass;
float tmp; float tmp;
for (pass = 1; pass <= npk - 1; pass++) { for (pass = 1; pass <= npk - 1; pass++) {
for (k = 0; k < npk - pass ; k++) { for (k = 0; k < npk - pass ; k++) {
@ -674,18 +634,18 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
signal vector. signal vector.
*/ */
int idrift,ifr,if0,ifd,k0; int32_t idrift,ifr,if0,ifd,k0;
int kindex; int32_t kindex;
float smax,ss,pow,p0,p1,p2,p3; float smax,ss,pow,p0,p1,p2,p3;
for(j=0; j<npk; j++) { //For each candidate... for(j=0; j<npk; j++) { //For each candidate...
smax=-1e30; smax=-1e30;
if0=freq0[j]/(DF/2.0)+256; if0=freq0[j]/(DF/2.0)+NSPERSYM;
for (ifr=if0-1; ifr<=if0+1; ifr++) { //Freq search for (ifr=if0-1; ifr<=if0+1; ifr++) { //Freq search
for( k0=-10; k0<22; k0++) { //Time search for( k0=-10; k0<22; k0++) { //Time search
for (idrift=-maxdrift; idrift<=maxdrift; idrift++) { //Drift search for (idrift=-maxdrift; idrift<=maxdrift; idrift++) { //Drift search
ss=0.0; ss=0.0;
pow=0.0; pow=0.0;
for (k=0; k<162; k++) { //Sum over symbols for (k=0; k<NSYM; k++) { //Sum over symbols
ifd=ifr+((float)k-81.0)/81.0*( (float)idrift )/DF; ifd=ifr+((float)k-81.0)/81.0*( (float)idrift )/DF;
kindex=k0+2*k; kindex=k0+2*k;
if( kindex < nffts ) { if( kindex < nffts ) {
@ -694,10 +654,10 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
p2=ps[ifd+1][kindex]; p2=ps[ifd+1][kindex];
p3=ps[ifd+3][kindex]; p3=ps[ifd+3][kindex];
p0=sqrt(p0); p0=sqrtf(p0);
p1=sqrt(p1); p1=sqrtf(p1);
p2=sqrt(p2); p2=sqrtf(p2);
p3=sqrt(p3); p3=sqrtf(p3);
ss=ss+(2*pr3[k]-1)*((p1+p3)-(p0+p2)); ss=ss+(2*pr3[k]-1)*((p1+p3)-(p0+p2));
pow=pow+p0+p1+p2+p3; pow=pow+p0+p1+p2+p3;
@ -708,7 +668,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
smax=sync1; smax=sync1;
shift0[j]=128*(k0+1); shift0[j]=128*(k0+1);
drift0[j]=idrift; drift0[j]=idrift;
freq0[j]=(ifr-256)*(DF/2.0); freq0[j]=(ifr-NSPERSYM)*(DF/2.0);
sync0[j]=sync1; sync0[j]=sync1;
} }
} }
@ -734,7 +694,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
*/ */
for (j=0; j<npk; j++) { for (j=0; j<npk; j++) {
memset(symbols,0,sizeof(char)*nbits*2); memset(symbols,0,sizeof(char)*NBITS*2);
memset(callsign,0,sizeof(char)*13); memset(callsign,0,sizeof(char)*13);
memset(call_loc_pow,0,sizeof(char)*23); memset(call_loc_pow,0,sizeof(char)*23);
memset(call,0,sizeof(char)*13); memset(call,0,sizeof(char)*13);
@ -751,7 +711,9 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
lagmin=shift1-144; lagmin=shift1-144;
lagmax=shift1+144; lagmax=shift1+144;
lagstep=8; lagstep=8;
if(options.quickmode) lagstep=16;
if(options.quickmode)
lagstep=16;
sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep, &shift1, sync_and_demodulate(idat, qdat, npoints, symbols, &f1, fstep, &shift1,
lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0); lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0);
@ -767,7 +729,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
worth_a_try = 0; worth_a_try = 0;
} }
int idt=0, ii=0, jiggered_shift; int32_t idt=0, ii=0, jiggered_shift;
double y,sq,rms; double y,sq,rms;
not_decoded=1; not_decoded=1;
@ -783,16 +745,15 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
&sync1, 2); &sync1, 2);
sq=0.0; sq=0.0;
for(i=0; i<162; i++) { for(i=0; i<NSYM; i++) {
y=(double)symbols[i] - 128.0; y=(double)symbols[i] - 128.0;
sq += y*y; sq += y*y;
} }
rms=sqrt(sq/162.0); rms=sqrt(sq/NSYM);
if((sync1 > minsync2) && (rms > minrms)) { if((sync1 > minsync2) && (rms > minrms)) {
deinterleave(symbols); deinterleave(symbols);
not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,NBITS,
not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
mettab,delta,maxcycles); mettab,delta,maxcycles);
} }
idt++; idt++;
@ -814,7 +775,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
noprint=unpk_(message,hashtab,call_loc_pow,call,loc,pwr,callsign); noprint=unpk_(message,hashtab,call_loc_pow,call,loc,pwr,callsign);
if( options.subtraction && (ipass == 0) && !noprint ) { if( options.subtraction && (ipass == 0) && !noprint ) {
unsigned char channel_symbols[162]; unsigned char channel_symbols[NSYM];
if( get_wspr_channel_symbols(call_loc_pow, hashtab, channel_symbols) ) { if( get_wspr_channel_symbols(call_loc_pow, hashtab, channel_symbols) ) {
subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, channel_symbols); subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
@ -825,7 +786,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
} }
// Remove dupes (same callsign and freq within 3 Hz) // Remove dupes (same callsign and freq within 3 Hz)
int dupe=0; int32_t dupe=0;
for (i=0; i<uniques; i++) { for (i=0; i<uniques; i++) {
if(!strcmp(callsign,allcalls[i]) && (fabs(f1-allfreqs[i]) <3.0)) if(!strcmp(callsign,allcalls[i]) && (fabs(f1-allfreqs[i]) <3.0))
dupe=1; dupe=1;
@ -880,7 +841,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
fftwf_free(fftin); fftwf_free(fftin);
fftwf_free(fftout); fftwf_free(fftout);
if ((fp_fftw_wisdom_file = fopen(wisdom_fname, "w"))) { if ((fp_fftw_wisdom_file = fopen("wspr_wisdom.dat", "w"))) {
fftwf_export_wisdom_to_file(fp_fftw_wisdom_file); fftwf_export_wisdom_to_file(fp_fftw_wisdom_file);
fclose(fp_fftw_wisdom_file); fclose(fp_fftw_wisdom_file);
} }
@ -890,7 +851,7 @@ int wspr_decode(float *idat, float *qdat, unsigned int npoints,
fftwf_destroy_plan(PLAN3); fftwf_destroy_plan(PLAN3);
if( options.usehashtable ) { if( options.usehashtable ) {
fhash=fopen(hash_fname,"w"); fhash=fopen("hashtable.txt","w");
for (i=0; i<32768; i++) { for (i=0; i<32768; i++) {
if( strncmp(hashtab+i*13,"\0",1) != 0 ) { if( strncmp(hashtab+i*13,"\0",1) != 0 ) {
fprintf(fhash,"%5d %s\n",i,hashtab+i*13); fprintf(fhash,"%5d %s\n",i,hashtab+i*13);

17
wsprd.h
Wyświetl plik

@ -46,7 +46,7 @@ struct decoder_options {
struct decoder_results { struct decoder_results {
float freq; double freq;
float sync; float sync;
float snr; float snr;
float dt; float dt;
@ -61,12 +61,13 @@ struct decoder_results {
void sync_and_demodulate(float *id, float *qd, long np, void sync_and_demodulate(float *id, float *qd, long np,
unsigned char *symbols, float *f1, float fstep, uint8_t *symbols, float *f1, float fstep,
int *shift1, int lagmin, int lagmax, int lagstep, int32_t *shift1, int32_t lagmin, int32_t lagmax, int32_t lagstep,
float *drift1, int symfac, float *sync, int mode); float *drift1, int32_t symfac, float *sync, int32_t mode);
void subtract_signal(float *id, float *qd, long np, void subtract_signal(float *id, float *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols); float f0, int32_t shift0, float drift0, uint8_t* channel_symbols);
void subtract_signal2(float *id, float *qd, long np, void subtract_signal2(float *id, float *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols); float f0, int32_t shift0, float drift0, uint8_t* channel_symbols);
int wspr_decode(float *idat, float *qdat, unsigned int npoints, int32_t wspr_decode(float *idat, float *qdat, uint32_t npoints,
struct decoder_options options, struct decoder_results *decodes, int *n_results); struct decoder_options options, struct decoder_results *decodes,
int32_t *n_results);