diff --git a/Makefile b/Makefile index 4db95f1d..348bd823 100644 --- a/Makefile +++ b/Makefile @@ -28,23 +28,48 @@ LIBSOURCES = fft_fftw.c libcsdr_wrapper.c #SOURCES = csdr.c $(LIBSOURCES) -cpufeature = $(if $(findstring $(1),$(shell cat /proc/cpuinfo)),$(2)) + +WSPR_SOURCES = $(wildcard wspr/*.c) + +# OSX see: +# github.com/sgentle/csdr/blob/osx/Makefile +# http://llvm.org/docs/Vectorizers.html +ifeq ($(shell uname -s),Darwin) + CPUFEATURES = $(shell sysctl -a | grep machdep.cpu.features | tr [A-Z] [a-z]) + PARAMS_CC = -DOSX -D_DARWIN_C_SOURCE -I/usr/local/include + PARAMS_LOOPVECT = -O3 -ffast-math +# uncomment to get messages about vectorization results +# PARAMS_LOOPVECT += -Rpass=loop-vectorize + PARAMS_LIBS += -L/usr/local/lib + PARAMS_SHLIB = -fpic -shared -o libcsdr.so.$(SOVERSION) + PARSEVECT ?= no + LDCONFIG = + PREFIX ?= /usr/local +else + CPUFEATURES = $(shell cat /proc/cpuinfo) + PARAMS_CC = + PARAMS_LOOPVECT = -O3 -ffast-math -fdump-tree-vect-details -dumpbase dumpvect + PARAMS_LIBS += -lrt + PARAMS_SHLIB = -fpic -shared -Wl,-soname,libcsdr.so.$(SOVERSION) -o libcsdr.so.$(SOVERSION) + PARSEVECT ?= yes + LDCONFIG = ldconfig || echo please run ldconfig + PREFIX ?= /usr +endif +cpufeature = $(if $(findstring $(1),$(CPUFEATURES)),$(2)) + PARAMS_SSE = $(call cpufeature,sse,-msse) $(call cpufeature,sse2,-msse2) $(call cpufeature,sse3,-msse3) $(call cpufeature,sse4a,-msse4a) $(call cpufeature,sse4_1,-msse4.1) $(call cpufeature,sse4_2,-msse4.2 -msse4) -mfpmath=sse PARAMS_NEON = -mfloat-abi=hard -march=armv7-a -mtune=cortex-a8 -mfpu=neon -mvectorize-with-neon-quad -funsafe-math-optimizations -Wformat=0 -DNEON_OPTS #tnx Jan Szumiec for the Raspberry Pi support PARAMS_RASPI = -mfloat-abi=hard -mcpu=arm1176jzf-s -mfpu=vfp -funsafe-math-optimizations -Wformat=0 PARAMS_ARM = $(if $(call cpufeature,BCM2708,dummy-text),$(PARAMS_RASPI),$(PARAMS_NEON)) PARAMS_SIMD = $(if $(call cpufeature,sse,dummy-text),$(PARAMS_SSE),$(PARAMS_ARM)) -PARAMS_LOOPVECT = -O3 -ffast-math -fdump-tree-vect-details -dumpbase dumpvect -PARAMS_LIBS = -g -lm -lrt -lfftw3f -DUSE_FFTW -DLIBCSDR_GPL -DUSE_IMA_ADPCM +PARAMS_LIBS += -g -lm -lfftw3f -DUSE_FFTW -DLIBCSDR_GPL -DUSE_IMA_ADPCM PARAMS_SO = -fpic PARAMS_MISC = -Wno-unused-result #DEBUG_ON = 0 #debug is always on by now (anyway it could be compiled with `make DEBUG_ON=1`) #PARAMS_DEBUG = $(if $(DEBUG_ON),-g,) FFTW_PACKAGE = fftw-3.3.3 -PREFIX ?= /usr SOVERSION = 0.15 -PARSEVECT ?= yes .PHONY: clean-vect clean codequality checkdocs v all: codequality csdr nmux @@ -53,17 +78,17 @@ libcsdr.so: fft_fftw.c fft_rpi.c libcsdr_wrapper.c libcsdr.c libcsdr_gpl.c fastd @echo Auto-detected optimization parameters: $(PARAMS_SIMD) @echo rm -f dumpvect*.vect - gcc -std=gnu99 $(PARAMS_LOOPVECT) $(PARAMS_SIMD) $(LIBSOURCES) $(PARAMS_LIBS) $(PARAMS_MISC) -fpic -shared -Wl,-soname,libcsdr.so.$(SOVERSION) -o libcsdr.so.$(SOVERSION) + gcc $(PARAMS_CC) -std=gnu99 $(PARAMS_LOOPVECT) $(PARAMS_SIMD) $(LIBSOURCES) $(PARAMS_LIBS) $(PARAMS_MISC) $(PARAMS_SHLIB) @ln -fs libcsdr.so.$(SOVERSION) libcsdr.so ifeq ($(PARSEVECT),yes) -./parsevect dumpvect*.vect endif csdr: csdr.c libcsdr.so - gcc -std=gnu99 $(PARAMS_LOOPVECT) $(PARAMS_SIMD) csdr.c $(PARAMS_LIBS) -L. -lcsdr $(PARAMS_MISC) -o csdr + gcc $(PARAMS_CC) -std=gnu99 $(PARAMS_LOOPVECT) $(PARAMS_SIMD) csdr.c $(WSPR_SOURCES) $(PARAMS_LIBS) -L. -lcsdr $(PARAMS_MISC) -o csdr ddcd: ddcd.cpp libcsdr.so ddcd.h - g++ $(PARAMS_LOOPVECT) $(PARAMS_SIMD) ddcd.cpp $(PARAMS_LIBS) -L. -lcsdr -lpthread $(PARAMS_MISC) -o ddcd + g++ $(PARAMS_CC) $(PARAMS_LOOPVECT) $(PARAMS_SIMD) ddcd.cpp $(PARAMS_LIBS) -L. -lcsdr -lpthread $(PARAMS_MISC) -o ddcd nmux: nmux.cpp libcsdr.so nmux.h tsmpool.cpp tsmpool.h - g++ $(PARAMS_LOOPVECT) $(PARAMS_SIMD) nmux.cpp tsmpool.cpp $(PARAMS_LIBS) -L. -lcsdr -lpthread $(PARAMS_MISC) -o nmux + g++ $(PARAMS_CC) $(PARAMS_LOOPVECT) $(PARAMS_SIMD) nmux.cpp tsmpool.cpp $(PARAMS_LIBS) -L. -lcsdr -lpthread $(PARAMS_MISC) -o nmux arm-cross: clean-vect #note: this doesn't work since having added FFTW arm-linux-gnueabihf-gcc -std=gnu99 -O3 -fshort-double -ffast-math -dumpbase dumpvect-arm -fdump-tree-vect-details -mfloat-abi=softfp -march=armv7-a -mtune=cortex-a9 -mfpu=neon -mvectorize-with-neon-quad -Wno-unused-result -Wformat=0 $(SOURCES) -lm -o ./csdr @@ -71,15 +96,16 @@ clean-vect: rm -f dumpvect*.vect clean: clean-vect rm -f libcsdr.so.$(SOVERSION) csdr ddcd nmux *.o *.so + rm -rf *.dSYM install: all install -m 0755 libcsdr.so.$(SOVERSION) $(PREFIX)/lib install -m 0755 csdr $(PREFIX)/bin install -m 0755 csdr-fm $(PREFIX)/bin install -m 0755 nmux $(PREFIX)/bin #-install -m 0755 ddcd $(PREFIX)/bin - @ldconfig || echo please run ldconfig + @$(LDCONFIG) uninstall: - rm $(PREFIX)/lib/libcsdr.so.$(SOVERSION) $(PREFIX)/bin/csdr $(PREFIX)/bin/csdr-fm + rm $(PREFIX)/lib/libcsdr.so.$(SOVERSION) $(PREFIX)/bin/csdr $(PREFIX)/bin/csdr-fm $(PREFIX)/bin/nmux ldconfig disasm: objdump -S libcsdr.so.$(SOVERSION) > libcsdr.disasm diff --git a/README.md b/README.md index 3ffb1c6b..6a12d4e8 100755 --- a/README.md +++ b/README.md @@ -952,7 +952,7 @@ Syntax: csdr add_n_zero_samples_at_beginning_f -When the function is executed, it furst writes `` 32-bit floating point zeros at the output, after that it just clones the input at the output. +When the function is executed, it first writes `` 32-bit floating point zeros at the output, after that it just clones the input at the output. ---- diff --git a/csdr.c b/csdr.c old mode 100755 new mode 100644 index 27dbf6dc..b827110e --- a/csdr.c +++ b/csdr.c @@ -53,6 +53,15 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "fastddc.h" #include +#ifdef OSX + #ifndef F_LINUX_SPECIFIC_BASE + #define F_LINUX_SPECIFIC_BASE 1024 + #endif + #ifndef F_SETPIPE_SZ + #define F_SETPIPE_SZ (F_LINUX_SPECIFIC_BASE + 7) + #endif +#endif + char usage[]= "csdr - a simple commandline tool for Software Defined Radio receiver DSP.\n\n" "usage: \n\n" @@ -62,10 +71,11 @@ char usage[]= " convert_f_u8\n" " convert_s8_f\n" " convert_f_s8\n" -" convert_f_s16\n" -" convert_s16_f\n" +" convert_f_s16 [--bigendian]\n" +" convert_s16_f [--bigendian]\n" " convert_f_s24 [--bigendian]\n" " convert_s24_f [--bigendian]\n" +" iq_swap_ff\n" " realpart_cf\n" " clipdetect_ff\n" " limit_ff [max_amplitude]\n" @@ -114,7 +124,6 @@ char usage[]= " flowcontrol \n" " through\n" " dsb_fc [q_value]\n" -" convert_f_samperf \n" " fmmod_fc\n" " fixed_amplitude_cc \n" " mono2stereo_s16\n" @@ -160,7 +169,7 @@ char usage[]= " bfsk_demod_cf \n" " normalized_timing_variance_u32_f [--debug]\n" " ?\n" -" ??\n" +" \?\?\n" " =\n" " shift_addfast_cc #only if system supports NEON \n" " shift_unroll_cc \n" @@ -201,7 +210,7 @@ int bigbufs = 0; char **argv_global; int argc_global; -int errhead() +void errhead() { fprintf(stderr, "%s%s%s: ", argv_global[0], ((argc_global>=2)?" ":""), ((argc_global>=2)?argv_global[1]:"")); } @@ -391,7 +400,7 @@ int sendbufsize(int size) return size; } -int parse_env() +void parse_env() { char* envtmp; envtmp=getenv("CSDR_DYNAMIC_BUFSIZE_ON"); @@ -530,6 +539,99 @@ int main(int argc, char *argv[]) } + #include "wspr/wspr.h" + if(!strcmp(argv[1], "wspr")) + { + int fft = (argc>2) && !strcmp(argv[2], "--fft"); + int demod = (argc>2) && !strcmp(argv[2], "--demod"); + + // csdr wspr --file filename data_rate reads_per_second + // essentially same scheme as csdr "flowcontrol" with addition of WSPR time sync + int file = (argc>3) && !strcmp(argv[2], "--file"); + time_t t; time(&t); struct tm tm; gmtime_r(&t, &tm); + bool sync = false; + int last_min = tm.tm_min, last_sec = tm.tm_sec; + + FILE *fp; + int flowcontrol_bufsize, flowcontrol_sleep; + unsigned char* flowcontrol_buffer; + + if (file) { + if(argc<6) return badsyntax("WSPR need required parameters (data_rate, reads_per_seconds)"); + int data_rate; + sscanf(argv[4],"%d",&data_rate); + int reads_per_second; + sscanf(argv[5],"%d",&reads_per_second); + flowcontrol_bufsize=ceil(1.*(double)data_rate/reads_per_second); + //if(!getbufsize()) return -2; + sendbufsize(flowcontrol_bufsize); + flowcontrol_buffer = (unsigned char*)malloc(sizeof(unsigned char)*flowcontrol_bufsize); + flowcontrol_sleep=floor(1000000./reads_per_second); + errhead(); fprintf(stderr, "WSPR --file flowcontrol_bufsize = %d, flowcontrol_sleep = %d\n", flowcontrol_bufsize, flowcontrol_sleep); + fp = fopen(argv[3], "r"); + } else { + if(!sendbufsize(initialize_buffers())) return -2; + } + + // csdr wspr --fft demod_pipe_name decimate + if (fft) { + if(argc<5) return badsyntax("WSPR need required parameter (demod_pipe_name, decimate)"); + int decimate=0; + sscanf(argv[4],"%d",&decimate); + SET_NONBLOCK(STDIN_FILENO); + wspr_init(argv[3], decimate, the_bufsize); + // doesn't return + } + + for(;;) + { + FEOF_CHECK; + + if (file) { + + // Allow samples to pass before the time sync so they appear on the FFT display. + // Rewind the file on time sync to position the samples correctly. + if (!sync) { + time(&t); gmtime_r(&t, &tm); + if (last_sec != tm.tm_sec) { + //fprintf(stderr, "WSPR --file %02d:%02d\n", tm.tm_min, tm.tm_sec); fflush(stderr); + + if ((last_min&1) && !(tm.tm_min&1)) { // wait for WSPR even 2 min boundary + //fprintf(stderr, "WSPR --file %02d:%02d START\n", tm.tm_min, tm.tm_sec); fflush(stderr); + rewind(fp); + sync = true; + } + + last_min = tm.tm_min; + last_sec = tm.tm_sec; + } + } + + //fprintf(stderr, "WSPR --file read fp=%p\n", fp); fflush(stderr); + int n = fread(flowcontrol_buffer, sizeof(unsigned char), flowcontrol_bufsize, fp); + //fprintf(stderr, "WSPR --file n=%d\n", n); fflush(stderr); + if (n && n != flowcontrol_bufsize) { + //fprintf(stderr, "WSPR --file n=%d != fcbs=%d\n", n, flowcontrol_bufsize); fflush(stderr); + } + if (n == 0) { + rewind(fp); + continue; + } + fwrite(flowcontrol_buffer, sizeof(unsigned char), n, stdout); + usleep(flowcontrol_sleep); + } else + + if (demod) { + // for now just passes the text output from the WSPR decoder on the FFT chain via the "iqtee2_pipe" + // (not used currently) + FREAD_C; + FWRITE_C; + } else + + TRY_YIELD; + } + return 0; + } if(!strcmp(argv[1],"convert_u8_f")) { @@ -581,26 +683,50 @@ int main(int argc, char *argv[]) } if((!strcmp(argv[1],"convert_f_i16")) || (!strcmp(argv[1],"convert_f_s16"))) { + int bigendian = (argc>2) && (!strcmp(argv[2],"--bigendian")); if(!sendbufsize(initialize_buffers())) return -2; - for(;;) - { - FEOF_CHECK; - FREAD_R; - convert_f_i16(input_buffer, buffer_i16, the_bufsize); - fwrite(buffer_i16, sizeof(short), the_bufsize, stdout); - TRY_YIELD; + if (bigendian) { + for(;;) + { + FEOF_CHECK; + FREAD_R; + convert_f_s16_big_endian(input_buffer, buffer_i16, the_bufsize); + fwrite(buffer_i16, sizeof(short), the_bufsize, stdout); + TRY_YIELD; + } + } else { + for(;;) + { + FEOF_CHECK; + FREAD_R; + convert_f_i16(input_buffer, buffer_i16, the_bufsize); + fwrite(buffer_i16, sizeof(short), the_bufsize, stdout); + TRY_YIELD; + } } } if((!strcmp(argv[1],"convert_i16_f")) || (!strcmp(argv[1],"convert_s16_f"))) { + int bigendian = (argc>2) && (!strcmp(argv[2],"--bigendian")); if(!sendbufsize(initialize_buffers())) return -2; - for(;;) - { - FEOF_CHECK; - fread(buffer_i16, sizeof(short), the_bufsize, stdin); - convert_i16_f(buffer_i16, output_buffer, the_bufsize); - FWRITE_R; - TRY_YIELD; + if (bigendian) { + for(;;) + { + FEOF_CHECK; + fread(buffer_i16, sizeof(short), the_bufsize, stdin); + convert_s16_big_endian_f(buffer_i16, output_buffer, the_bufsize); + FWRITE_R; + TRY_YIELD; + } + } else { + for(;;) + { + FEOF_CHECK; + fread(buffer_i16, sizeof(short), the_bufsize, stdin); + convert_i16_f(buffer_i16, output_buffer, the_bufsize); + FWRITE_R; + TRY_YIELD; + } } } if(!strcmp(argv[1],"convert_f_s24")) @@ -631,6 +757,19 @@ int main(int argc, char *argv[]) TRY_YIELD; } } + if(!strcmp(argv[1],"iq_swap_ff")) + { + double t; + if(!sendbufsize(initialize_buffers())) return -2; + for(;;) + { + FEOF_CHECK; + FREAD_C; + for(int i=0;i)"); errhead(); fprintf(stderr, "initial squelch level is %g\n", squelch_level); if((argc<=5)||((argc>5)&&(strcmp(argv[4],"--outfifo")))) return badsyntax("need required parameter (--outfifo )"); @@ -2305,7 +2444,7 @@ int main(int argc, char *argv[]) int plusarg=0; int fd; - if(fd=init_fifo(argc,argv)) + if((fd=init_fifo(argc,argv))) { while(!read_fifo_ctl(fd,"%g\n",&shift_rate)) usleep(10000); plusarg=1; @@ -2362,7 +2501,7 @@ int main(int argc, char *argv[]) complexf* output = (complexf*)fft_malloc(sizeof(complexf)*ddc.post_input_size); decimating_shift_addition_status_t shift_stat; - bzero(&shift_stat, sizeof(shift_stat)); + memset(&shift_stat, 0, sizeof(shift_stat)); for(;;) { FEOF_CHECK; @@ -2766,7 +2905,7 @@ int main(int argc, char *argv[]) { if(!initialize_buffers()) return -2; sendbufsize(1); - char local_input_buffer[8]; + unsigned char local_input_buffer[8]; for(;;) { FEOF_CHECK; @@ -3346,7 +3485,7 @@ int main(int argc, char *argv[]) while(current_buffer_write_cntr!=current_buffer_read_cntr) { int result = fwrite(async_tee_buffers+(the_bufsize*current_buffer_write_cntr)+current_byte_write_cntr, sizeof(unsigned char), the_bufsize-current_byte_write_cntr, teefile); - if(!result) { errhead(); fprintf(stderr, "\t fwrite tee zero, next turn\n"); break; } + if(!result) { errhead(); fprintf(stderr, "\t fwrite %s tee zero, next turn\n", argv[2]); break; } current_byte_write_cntr += result; //errhead(); fprintf(stderr, "\tfwrite tee, current_byte_write_cntr = %d, current_buffer_write_cntr = %d, current_buffer_read_cntr = %d\n", // current_byte_write_cntr, current_buffer_write_cntr, current_buffer_read_cntr); @@ -3370,7 +3509,7 @@ int main(int argc, char *argv[]) float rate; int fd; - if(fd=init_fifo(argc,argv)) + if((fd=init_fifo(argc,argv))) { while(!read_fifo_ctl(fd,"%g\n",&rate)) usleep(10000); } diff --git a/libcsdr.c b/libcsdr.c old mode 100755 new mode 100644 index 1eb91b38..67ca02a2 --- a/libcsdr.c +++ b/libcsdr.c @@ -1537,7 +1537,7 @@ char psk31_varicode_decoder_push(unsigned long long* status_shr, unsigned char s { *status_shr=((*status_shr)<<1)|(!!symbol); //shift new bit in shift register //fprintf(stderr,"*status_shr = %llx\n", *status_shr); - if((*status_shr)&0xFFF==0) return 0; + if(((*status_shr)&0xFFF)==0) return 0; for(int i=0;isamples_per_symbol/2) sinearest++; unsigned socorrect = initial_sample_offset+(sinearest*samples_per_symbol); //the sample offset which input[i] should have been, in order to sample at the maximum effect point - int sodiff = abs(socorrect-input[i]); + // FIXME: is this fix correct? + // A compiler correctly complains the abs() of the difference of two unsigned quantities has no effect. + // But I'm not sure if simply casting the difference to an int gives the correct result. (jks) + //int sodiff = abs(socorrect-input[i]); + int sodiff = abs((int) (socorrect-input[i])); float ndiff = (float)sodiff/samples_per_symbol; ndiff_rad[i] = ndiff*PI; @@ -2375,6 +2379,15 @@ void convert_s16_f(short* input, float* output, int input_size) for(int i=0;i> 8) & 0xff); + output[i]=(float)sint/SHRT_MAX; //@convert_s16_big_endian_f + } +} + void convert_f_u8(float* input, unsigned char* output, int input_size) { for(int i=0;i> 8) & 0xff); + output[i]=sint; //@convert_f_s16_big_endian + } +} + void convert_i16_f(short* input, float* output, int input_size) { convert_s16_f(input, output, input_size); } void convert_f_i16(float* input, short* output, int input_size) { convert_f_s16(input, output, input_size); } @@ -2470,7 +2492,7 @@ int deinit_get_random_samples_f(FILE* status) return fclose(status); } -int firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol) +void firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol) { //needs a taps_length 2 × samples_per_symbol + 1 int middle_i=taps_length/2; @@ -2479,7 +2501,7 @@ int firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol) normalize_fir_f(taps, taps, taps_length); } -int firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta) +void firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta) { //needs an odd taps_length int middle_i=taps_length/2; @@ -2515,12 +2537,12 @@ matched_filter_type_t matched_filter_get_type_from_string(char* input) return MATCHED_FILTER_DEFAULT; } -float* add_ff(float* input1, float* input2, float* output, int input_size) +void add_ff(float* input1, float* input2, float* output, int input_size) { for(int i=0;i +#include + +int csdr_clock_gettime(clockid_t clock_is, struct timespec *tp) +{ + clock_serv_t cclock; + mach_timespec_t mts; + host_get_clock_service(mach_host_self(), CALENDAR_CLOCK, &cclock); + clock_get_time(cclock, &mts); + mach_port_deallocate(mach_task_self(), cclock); + tp->tv_sec = mts.tv_sec; + tp->tv_nsec = mts.tv_nsec; + return 0; +} +#endif diff --git a/libcsdr.h b/libcsdr.h index 7e472056..67756f6c 100644 --- a/libcsdr.h +++ b/libcsdr.h @@ -64,6 +64,13 @@ typedef struct complexf_s { float i; float q; } complexf; //they dropped M_PI in C99, so we define it: #define PI ((float)3.14159265358979323846) +#ifdef OSX + //clock_gettime is supposed to be included starting with OSX 10.12 (Sierra) + #include + #define clock_gettime(clock_id, tp) csdr_clock_gettime(clock_id, tp) + int csdr_clock_gettime(clockid_t clock_id, struct timespec *tp); +#endif + #define TIME_TAKEN(start,end) ((end.tv_sec-start.tv_sec)+(end.tv_nsec-start.tv_nsec)/1e9) //window @@ -222,7 +229,9 @@ void convert_f_u8(float* input, unsigned char* output, int input_size); void convert_s8_f(signed char* input, float* output, int input_size); void convert_f_s8(float* input, signed char* output, int input_size); void convert_f_s16(float* input, short* output, int input_size); +void convert_f_s16_big_endian(float* input, short* output, int input_size); void convert_s16_f(short* input, float* output, int input_size); +void convert_s16_big_endian_f(short* input, float* output, int input_size); void convert_f_i16(float* input, short* output, int input_size); void convert_i16_f(short* input, float* output, int input_size); void convert_f_s24(float* input, unsigned char* output, int input_size, int bigendian); @@ -386,7 +395,7 @@ FILE* init_get_random_samples_f(); void get_random_samples_f(float* output, int output_size, FILE* status); void get_random_gaussian_samples_c(complexf* output, int output_size, FILE* status); int deinit_get_random_samples_f(FILE* status); -float* add_ff(float* input1, float* input2, float* output, int input_size); +void add_ff(float* input1, float* input2, float* output, int input_size); float total_logpower_cf(complexf* input, int input_size); float normalized_timing_variance_u32_f(unsigned* input, float* temp, int input_size, int samples_per_symbol, int initial_sample_offset, int debug_print); @@ -398,14 +407,14 @@ typedef enum matched_filter_type_e #define MATCHED_FILTER_DEFAULT MATCHED_FILTER_RRC -int firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol); -int firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta); +void firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol); +void firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta); matched_filter_type_t matched_filter_get_type_from_string(char* input); int apply_real_fir_cc(complexf* input, complexf* output, int input_size, float* taps, int taps_length); void generic_slicer_f_u8(float* input, unsigned char* output, int input_size, int n_symbols); void plain_interpolate_cc(complexf* input, complexf* output, int input_size, int interpolation);; void normalize_fir_f(float* input, float* output, int length); -float* add_const_cc(complexf* input, complexf* output, int input_size, complexf x); +void add_const_cc(complexf* input, complexf* output, int input_size, complexf x); void pack_bits_1to8_u8_u8(unsigned char* input, unsigned char* output, int input_size); unsigned char pack_bits_8to1_u8_u8(unsigned char* input); void dbpsk_decoder_c_u8(complexf* input, unsigned char* output, int input_size); diff --git a/nmux.cpp b/nmux.cpp index 7d750124..e8c685e6 100644 --- a/nmux.cpp +++ b/nmux.cpp @@ -30,6 +30,10 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "nmux.h" +#ifdef OSX + #define MSG_NOSIGNAL 0 // done via SO_NOSIGPIPE +#endif + char help_text[]="nmux is a TCP stream multiplexer. It reads data from the standard input, and sends it to each client connected through TCP sockets. Available command line options are:\n" "\t--port (-p), --address (-a): TCP port and address to listen.\n" "\t--bufsize (-b), --bufcnt (-n): Internal buffer size and count.\n" @@ -128,6 +132,11 @@ int main(int argc, char* argv[]) if( setsockopt(listen_socket, SOL_SOCKET, SO_REUSEADDR, (char *)&sockopt, sizeof(sockopt)) == -1 ) error_exit(MSG_START "cannot set SO_REUSEADDR"); //the best description on SO_REUSEADDR ever: http://stackoverflow.com/a/14388707/3182453 +#ifdef OSX + if( setsockopt(listen_socket, SOL_SOCKET, SO_NOSIGPIPE, (char *)&sockopt, sizeof(sockopt)) == -1 ) + error_exit(MSG_START "cannot set SO_NOSIGPIPE"); +#endif + memset(&addr_host,'0',sizeof(addr_host)); addr_host.sin_family = AF_INET; addr_host.sin_port = htons(host_port); @@ -192,7 +201,7 @@ int main(int argc, char* argv[]) if(NMUX_DEBUG) { fprintf(stderr, "\x1b[1m\x1b[33mmainfor: clients before closing: "); - for(int i=0;ithread, NULL, client_thread, (void*)new_client)==0) { clients.push_back(new_client); - fprintf(stderr, MSG_START "pthread_create() done, clients now: %d\n", clients.size()); + fprintf(stderr, MSG_START "pthread_create() done, clients now: %ld\n", clients.size()); } else { @@ -278,14 +287,14 @@ int main(int argc, char* argv[]) void* client_thread (void* param) { - fprintf(stderr, "client 0x%x: started!\n", (intptr_t)param); + fprintf(stderr, "client 0x%lx: started!\n", (intptr_t)param); client_t* this_client = (client_t*)param; this_client->status = CS_THREAD_RUNNING; int retval; tsmpool* lpool = this_client->lpool; - if(NMUX_DEBUG) fprintf(stderr, "client 0x%x: socket = %d!\n", (intptr_t)param, this_client->socket); + if(NMUX_DEBUG) fprintf(stderr, "client 0x%lx: socket = %d!\n", (intptr_t)param, this_client->socket); - if(NMUX_DEBUG) fprintf(stderr, "client 0x%x: poll init...", (intptr_t)param); + if(NMUX_DEBUG) fprintf(stderr, "client 0x%lx: poll init...", (intptr_t)param); struct pollfd pollfds[1]; pollfds[0].fd = this_client->socket; pollfds[0].events = POLLOUT; @@ -307,10 +316,10 @@ void* client_thread (void* param) // (Wait for the server process to wake me up.) while(!pool_read_buffer || client_buffer_index >= lpool->size) { - if(NMUX_DEBUG) fprintf(stderr, "client 0x%x: trying to grb\n", (intptr_t)param); + if(NMUX_DEBUG) fprintf(stderr, "client 0x%lx: trying to grb\n", (intptr_t)param); pool_read_buffer = (char*)lpool->get_read_buffer(this_client->tsmthread); if(pool_read_buffer) { client_buffer_index = 0; break; } - if(NMUX_DEBUG) fprintf(stderr, "client 0x%x: cond_waiting for more data\n", (intptr_t)param); + if(NMUX_DEBUG) fprintf(stderr, "client 0x%lx: cond_waiting for more data\n", (intptr_t)param); pthread_mutex_lock(&wait_mutex); this_client->sleeping = 1; pthread_cond_wait(&wait_condition, &wait_mutex); @@ -318,14 +327,14 @@ void* client_thread (void* param) } //Wait for the socket to be available for write. - if(NMUX_DEBUG) fprintf(stderr, "client 0x%x: polling for socket write...", (intptr_t)param); + if(NMUX_DEBUG) fprintf(stderr, "client 0x%lx: polling for socket write...", (intptr_t)param); int ret = poll(pollfds, 1, -1); if(NMUX_DEBUG) fprintf(stderr, "client polled for socket write.\n"); if(ret == 0) continue; else if (ret == -1) { client_goto_source = 1; goto client_thread_exit; } //Read data from global tsmpool and write it to client socket - if(NMUX_DEBUG) fprintf(stderr, "client 0x%x: sending...", (intptr_t)param); + if(NMUX_DEBUG) fprintf(stderr, "client 0x%lx: sending...", (intptr_t)param); ret = send(this_client->socket, pool_read_buffer + client_buffer_index, lpool->size - client_buffer_index, MSG_NOSIGNAL); if(NMUX_DEBUG) fprintf(stderr, "client sent.\n"); if(ret == -1) @@ -340,7 +349,7 @@ void* client_thread (void* param) } client_thread_exit: - fprintf(stderr, "client 0x%x: CS_THREAD_FINISHED, client_goto_source = %d, errno = %d", (intptr_t)param, client_goto_source, errno); + fprintf(stderr, "client 0x%lx: CS_THREAD_FINISHED, client_goto_source = %d, errno = %d", (intptr_t)param, client_goto_source, errno); this_client->status = CS_THREAD_FINISHED; pthread_exit(NULL); return NULL; diff --git a/tsmpool.cpp b/tsmpool.cpp index d70e8bae..7c8f88f2 100644 --- a/tsmpool.cpp +++ b/tsmpool.cpp @@ -42,7 +42,7 @@ tsmthread_t* tsmpool::register_thread() return thread; } -int tsmpool::remove_thread(tsmthread_t* thread) +void tsmpool::remove_thread(tsmthread_t* thread) { pthread_mutex_lock(&this->mutex); for(int i=0;i +#include +#include +#include "fano.h" + +struct node { + unsigned long encstate; // Encoder state of next node + long gamma; // Cumulative metric to this node + int metrics[4]; // Metrics indexed by all possible tx syms + int tm[2]; // Sorted metrics for current hypotheses + int i; // Current branch being tested +}; + +// Convolutional coding polynomials. All are rate 1/2, K=32 +#ifdef NASA_STANDARD +/* "NASA standard" code by Massey & Costello + * Nonsystematic, quick look-in, dmin=11, dfree=23 + * used on Pioneer 10-12, Helios A,B + */ +#define POLY1 0xbbef6bb7 +#define POLY2 0xbbef6bb5 +#endif + +#ifdef MJ +/* Massey-Johannesson code + * Nonsystematic, quick look-in, dmin=13, dfree>=23 + * Purported to be more computationally efficient than Massey-Costello + */ +#define POLY1 0xb840a20f +#define POLY2 0xb840a20d +#endif + +#ifdef LL +/* Layland-Lushbaugh code + * Nonsystematic, non-quick look-in, dmin=?, dfree=? + */ +#define POLY1 0xf2d05351 +#define POLY2 0xe4613c47 +#endif + +/* Convolutionally encode a packet. The input data bytes are read + * high bit first and the encoded packet is written into 'symbols', + * one symbol per byte. The first symbol is generated from POLY1, + * the second from POLY2. + * + * Storing only one symbol per byte uses more space, but it is faster + * and easier than trying to pack them more compactly. + */ + #if 0 +int encode( + unsigned char *symbols, // Output buffer, 2*nbytes*8 + unsigned char *data, // Input buffer, nbytes + unsigned int nbytes) // Number of bytes in data +{ + unsigned long encstate; + int sym; + int i; + + encstate = 0; + while(nbytes-- != 0) { + for(i=7;i>=0;i--) { + encstate = (encstate << 1) | ((*data >> i) & 1); + ENCODE(sym,encstate); + *symbols++ = sym >> 1; + *symbols++ = sym & 1; + } + data++; + } + return 0; +} +#endif + +/* Decode packet with the Fano algorithm. + * Return 1 on success, 0 on timeout + */ +int fano( + unsigned int *metric, // Final path metric (returned value) + unsigned int *cycles, // Cycle count (returned value) + unsigned int *maxnp, // Progress before timeout (returned value) + unsigned char *data, // Decoded output data + unsigned char *symbols, // Raw deinterleaved input symbols + unsigned int nbits, // Number of output bits + int mettab[2][256], // Metric table, [sent sym][rx symbol] + int delta, // Threshold adjust parameter + unsigned int maxcycles) // Decoding timeout in cycles per bit +{ + struct node *nodes; // First node + struct node *np; // Current node + struct node *lastnode; // Last node + struct node *tail; // First node of tail + int t; // Threshold + int m0,m1; + int ngamma; + unsigned int lsym; + unsigned int i; + + if((nodes = (struct node *)malloc((nbits+1)*sizeof(struct node))) == NULL) { + printf("malloc failed\n"); + return 0; + } + lastnode = &nodes[nbits-1]; + tail = &nodes[nbits-31]; + *maxnp = 0; + +/* Compute all possible branch metrics for each symbol pair + * This is the only place we actually look at the raw input symbols + */ + for(np=nodes;np <= lastnode;np++) { + np->metrics[0] = mettab[0][symbols[0]] + mettab[0][symbols[1]]; + np->metrics[1] = mettab[0][symbols[0]] + mettab[1][symbols[1]]; + np->metrics[2] = mettab[1][symbols[0]] + mettab[0][symbols[1]]; + np->metrics[3] = mettab[1][symbols[0]] + mettab[1][symbols[1]]; + symbols += 2; + } + np = nodes; + np->encstate = 0; + +// Compute and sort branch metrics from root node */ + ENCODE(lsym,np->encstate); // 0-branch (LSB is 0) + m0 = np->metrics[lsym]; + +/* Now do the 1-branch. To save another ENCODE call here and + * inside the loop, we assume that both polynomials are odd, + * providing complementary pairs of branch symbols. + + * This code should be modified if a systematic code were used. + */ + + m1 = np->metrics[3^lsym]; + if(m0 > m1) { + np->tm[0] = m0; // 0-branch has better metric + np->tm[1] = m1; + } else { + np->tm[0] = m1; // 1-branch is better + np->tm[1] = m0; + np->encstate++; // Set low bit + } + np->i = 0; // Start with best branch + maxcycles *= nbits; + np->gamma = t = 0; + + // Start the Fano decoder + for(i=1;i <= maxcycles;i++) { + TRY_YIELD_DECODE; + if((int)(np-nodes) > (int)*maxnp) *maxnp=(int)(np-nodes); +#ifdef debug + printf("k=%ld, g=%ld, t=%d, m[%d]=%d, maxnp=%d, encstate=%lx\n", + np-nodes,np->gamma,t,np->i,np->tm[np->i],*maxnp,np->encstate); +#endif +// Look forward */ + ngamma = np->gamma + np->tm[np->i]; + if(ngamma >= t) { + if(np->gamma < t + delta) { // Node is acceptable + /* First time we've visited this node; + * Tighten threshold. + * + * This loop could be replaced with + * t += delta * ((ngamma - t)/delta); + * but the multiply and divide are slower. + */ + while(ngamma >= t + delta) t += delta; + } + np[1].gamma = ngamma; // Move forward + np[1].encstate = np->encstate << 1; + if( ++np == (lastnode+1) ) { + break; // Done! + } + + /* Compute and sort metrics, starting with the + * zero branch + */ + ENCODE(lsym,np->encstate); + if(np >= tail) { + /* The tail must be all zeroes, so don't + * bother computing the 1-branches here. + */ + np->tm[0] = np->metrics[lsym]; + } else { + m0 = np->metrics[lsym]; + m1 = np->metrics[3^lsym]; + if(m0 > m1) { + np->tm[0] = m0; // 0-branch is better + np->tm[1] = m1; + } else { + np->tm[0] = m1; // 1-branch is better + np->tm[1] = m0; + np->encstate++; // Set low bit + } + } + np->i = 0; // Start with best branch + continue; + } + // Threshold violated, can't go forward + for(;;) { // Look backward + if(np == nodes || np[-1].gamma < t) { + /* Can't back up either. + * Relax threshold and and look + * forward again to better branch. + */ + t -= delta; + if(np->i != 0) { + np->i = 0; + np->encstate ^= 1; + } + break; + } + // Back up + if(--np < tail && np->i != 1) { + np->i++; // Search next best branch + np->encstate ^= 1; + break; + } // else keep looking back + } + } + *metric = np->gamma; // Return the final path metric + + // Copy decoded data to user's buffer + nbits >>= 3; + np = &nodes[7]; + while(nbits-- != 0) { + *data++ = np->encstate; + np += 8; + } + *cycles = i+1; + + free(nodes); + if(i >= maxcycles) return 0; // Decoder timed out + return 1; // Successful completion +} diff --git a/wspr/fano.h b/wspr/fano.h new file mode 100644 index 00000000..3290a094 --- /dev/null +++ b/wspr/fano.h @@ -0,0 +1,39 @@ +/* + This file is part of wsprd. + + File name: fano.h + + Description: Header file for sequential Fano decoder. + + Copyright 1994, Phil Karn, KA9Q + Minor modifications by Joe Taylor, K1JT +*/ + +#ifndef FANO_H +#define FANO_H + +int fano(unsigned int *metric, unsigned int *cycles, unsigned int *maxnp, + unsigned char *data,unsigned char *symbols, unsigned int nbits, + int mettab[2][256],int delta,unsigned int maxcycles); + +int encode(unsigned char *symbols,unsigned char *data,unsigned int nbytes); + +extern unsigned char Partab[]; + +/* Convolutional encoder macro. Takes the encoder state, generates + * a rate 1/2 symbol pair and stores it in 'sym'. The symbol generated from + * POLY1 goes into the 2-bit of sym, and the symbol generated from POLY2 + * goes into the 1-bit. + */ +#define ENCODE(sym,encstate){\ +unsigned long _tmp;\ +\ +_tmp = (encstate) & POLY1;\ +_tmp ^= _tmp >> 16;\ +(sym) = Partab[(_tmp ^ (_tmp >> 8)) & 0xff] << 1;\ +_tmp = (encstate) & POLY2;\ +_tmp ^= _tmp >> 16;\ +(sym) |= Partab[(_tmp ^ (_tmp >> 8)) & 0xff];\ +} + +#endif diff --git a/wspr/jelinek.c b/wspr/jelinek.c new file mode 100644 index 00000000..c9dd30ed --- /dev/null +++ b/wspr/jelinek.c @@ -0,0 +1,164 @@ +/* + Soft-decision stack-based sequential decoder for K=32 r=1/2 + convolutional code. This code implements the "stack-bucket" algorithm + described in: + "Fast Sequential Decoding Algorithm Using a Stack", F. Jelinek + + The ENCODE macro from Phil Karn's (KA9Q) Fano decoder is used. + + Written by Steve Franke, K9AN for WSJT-X (July 2015) + */ + +#include "jelinek.h" + +#include +#include +#include +#include /* memset */ + +#include "fano.h" + +/* WSPR uses the Layland-Lushbaugh code + * Nonsystematic, non-quick look-in, dmin=?, dfree=? + */ +#define POLY1 0xf2d05351 +#define POLY2 0xe4613c47 + +//Decoder - returns 1 on success, 0 on timeout +int jelinek( + unsigned int *metric, /* Final path metric (returned value) */ + unsigned int *cycles, /* Cycle count (returned value) */ + unsigned char *data, /* Decoded output data */ + unsigned char *symbols, /* Raw deinterleaved input symbols */ + unsigned int nbits, /* Number of output bits */ + unsigned int stacksize, + struct snode *stack, + int mettab[2][256], /* Metric table, [sent sym][rx symbol] */ + unsigned int maxcycles)/* Decoding timeout in cycles per bit */ +{ + + // Compute branch metrics for each symbol pair + // The sequential decoding algorithm only uses the metrics, not the + // symbol values. + unsigned int i; + long int metrics[81][4]; + for(i=0; i>5)+200; //fast, but not particularly safe - totmet can be negative + if( bucket > high_bucket ) high_bucket=bucket; + if( bucket < low_bucket ) low_bucket=bucket; + + // place the 0 node on the stack, overwriting the parent (current) node + stack[ptr].encstate=encstate; + stack[ptr].gamma=totmet0; + stack[ptr].depth=depth; + stack[ptr].jpointer=buckets[bucket]; + buckets[bucket]=ptr; + + // if in the tail, only need to evaluate the "0" branch. + // Otherwise, enter this "if" and place the 1 node on the stack, + if( depth <= nbits_minus_ntail ) { + if( stackptr < stacksize_minus_1 ) { + stackptr++; + ptr=stackptr; + } else { // stack full + while( buckets[low_bucket] == 0 ) { //write latest to where the top of the lowest bucket points + low_bucket++; + } + ptr=buckets[low_bucket]; + buckets[low_bucket]=stack[ptr].jpointer; //make bucket point to next older entry + } + + bucket=(totmet1>>5)+200; //this may not be safe on all compilers + if( bucket > high_bucket ) high_bucket=bucket; + if( bucket < low_bucket ) low_bucket=bucket; + + stack[ptr].encstate=encstate+1; + stack[ptr].gamma=totmet1; + stack[ptr].depth=depth; + stack[ptr].jpointer=buckets[bucket]; + buckets[bucket]=ptr; + } + + // pick off the latest entry from the high bucket + while( buckets[high_bucket] == 0 ) { + high_bucket--; + } + ptr=buckets[high_bucket]; + buckets[high_bucket]=stack[ptr].jpointer; + depth=stack[ptr].depth; + gamma=stack[ptr].gamma; + encstate=stack[ptr].encstate; + + // we are done if the top entry on the stack is at depth nbits + if (depth == nbits) { + break; + } + } + + *cycles = i+1; + *metric = gamma; /* Return final path metric */ + + // printf("cycles %d stackptr=%d, depth=%d, gamma=%d, encstate=%lx\n", + // *cycles, stackptr, depth, *metric, encstate); + + for (i=0; i<7; i++) { + data[i]=(encstate>>(48-i*8))&(0x00000000000000ff); + } + for (i=7; i<11; i++) { + data[i]=0; + } + + if(*cycles/nbits >= maxcycles) //timed out + { + return 0; + } + return 1; //success +} diff --git a/wspr/jelinek.h b/wspr/jelinek.h new file mode 100644 index 00000000..9c22cf08 --- /dev/null +++ b/wspr/jelinek.h @@ -0,0 +1,34 @@ +/* + This file is part of wsprd. + + File name: jelinek.h + + Description: Header file for sequential Jelinek decoder. + + Written by Steve Franke, K9AN for WSJT-X (July 2015) +*/ + +#ifndef JELINEK_H +#define JELINEK_H + +#include + +struct snode { + uint64_t encstate; // Encoder state + int gamma; // Cumulative metric to this node + unsigned int depth; // depth of this node + unsigned int jpointer; +}; + +int jelinek(unsigned int *metric, + unsigned int *cycles, + unsigned char *data, + unsigned char *symbols, + unsigned int nbits, + unsigned int stacksize, + struct snode *stack, + int mettab[2][256], + unsigned int maxcycles); + +#endif + diff --git a/wspr/metric_tables.h b/wspr/metric_tables.h new file mode 100644 index 00000000..973651d7 --- /dev/null +++ b/wspr/metric_tables.h @@ -0,0 +1,138 @@ +/******************************************************************************* +* 4 metric tables calculated via simulation for 2-FSK with Es/No=0,3,6,9 dB +* tables were calculated for constant rms noise level of 50. The symbol vector +* should be normalized to have rms amplitude equal to "symbol_scale". +********************************************************************************/ +//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, + 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.9196, 0.9171, 0.9154, 0.9123, 0.9076, 0.9061, 0.9030, 0.9000, 0.8965, 0.8934, + 0.8903, 0.8874, 0.8834, 0.8792, 0.8760, 0.8726, 0.8685, 0.8639, 0.8599, 0.8550, + 0.8504, 0.8459, 0.8422, 0.8364, 0.8320, 0.8262, 0.8215, 0.8159, 0.8111, 0.8052, + 0.7996, 0.7932, 0.7878, 0.7812, 0.7745, 0.7685, 0.7616, 0.7550, 0.7479, 0.7405, + 0.7336, 0.7255, 0.7184, 0.7102, 0.7016, 0.6946, 0.6860, 0.6769, 0.6687, 0.6598, + 0.6503, 0.6416, 0.6325, 0.6219, 0.6122, 0.6016, 0.5920, 0.5818, 0.5711, 0.5606, + 0.5487, 0.5374, 0.5266, 0.5142, 0.5020, 0.4908, 0.4784, 0.4663, 0.4532, 0.4405, + 0.4271, 0.4144, 0.4006, 0.3865, 0.3731, 0.3594, 0.3455, 0.3304, 0.3158, 0.3009, + 0.2858, 0.2708, 0.2560, 0.2399, 0.2233, 0.2074, 0.1919, 0.1756, 0.1590, 0.1427, + 0.1251, 0.1074, 0.0905, 0.0722, 0.0550, 0.0381, 0.0183, 0.0000, -0.0185, -0.0391, + -0.0571, -0.0760, -0.0966, -0.1160, -0.1370, -0.1584, -0.1787, -0.1999, -0.2214, -0.2423, + -0.2643, -0.2879, -0.3114, -0.3336, -0.3568, -0.3806, -0.4050, -0.4293, -0.4552, -0.4798, + -0.5046, -0.5296, -0.5564, -0.5836, -0.6093, -0.6372, -0.6645, -0.6933, -0.7208, -0.7495, + -0.7763, -0.8065, -0.8378, -0.8660, -0.8964, -0.9293, -0.9592, -0.9907, -1.0214, -1.0509, + -1.0850, -1.1168, -1.1528, -1.1847, -1.2157, -1.2511, -1.2850, -1.3174, -1.3540, -1.3900, + -1.4201, -1.4580, -1.4956, -1.5292, -1.5683, -1.6030, -1.6411, -1.6789, -1.7147, -1.7539, + -1.7887, -1.8289, -1.8699, -1.9043, -1.9469, -1.9849, -2.0267, -2.0610, -2.1028, -2.1391, + -2.1855, -2.2215, -2.2712, -2.3033, -2.3440, -2.3870, -2.4342, -2.4738, -2.5209, -2.5646, + -2.6016, -2.6385, -2.6868, -2.7356, -2.7723, -2.8111, -2.8524, -2.9009, -2.9428, -2.9879, + -3.0103, -3.0832, -3.1340, -3.1628, -3.2049, -3.2557, -3.3101, -3.3453, -3.4025, -3.4317, + -3.4828, -3.5270, -3.5745, -3.6181, -3.6765, -3.7044, -3.7410, -3.8118, -3.8368, -3.9549, + -3.9488, -3.9941, -4.0428, -4.0892, -4.1648, -4.1965, -4.1892, -4.2565, -4.3356, -4.3948, + -4.4481, -4.4607, -4.5533, -4.5809, -4.5927, -5.1047}, + {0.9978, 0.9962, 0.9961, 0.9959, 0.9958, 0.9954, 0.9949, 0.9950, 0.9947, 0.9942, + 0.9940, 0.9939, 0.9933, 0.9931, 0.9928, 0.9924, 0.9921, 0.9916, 0.9911, 0.9909, + 0.9903, 0.9900, 0.9892, 0.9887, 0.9883, 0.9877, 0.9869, 0.9863, 0.9857, 0.9848, + 0.9842, 0.9835, 0.9825, 0.9817, 0.9808, 0.9799, 0.9791, 0.9777, 0.9767, 0.9757, + 0.9744, 0.9729, 0.9716, 0.9704, 0.9690, 0.9674, 0.9656, 0.9641, 0.9625, 0.9609, + 0.9587, 0.9567, 0.9548, 0.9524, 0.9501, 0.9478, 0.9453, 0.9426, 0.9398, 0.9371, + 0.9339, 0.9311, 0.9277, 0.9242, 0.9206, 0.9168, 0.9131, 0.9087, 0.9043, 0.8999, + 0.8953, 0.8907, 0.8857, 0.8803, 0.8747, 0.8690, 0.8632, 0.8572, 0.8507, 0.8439, + 0.8368, 0.8295, 0.8217, 0.8138, 0.8058, 0.7972, 0.7883, 0.7784, 0.7694, 0.7597, + 0.7489, 0.7378, 0.7269, 0.7152, 0.7030, 0.6911, 0.6782, 0.6643, 0.6506, 0.6371, + 0.6211, 0.6054, 0.5897, 0.5740, 0.5565, 0.5393, 0.5214, 0.5027, 0.4838, 0.4643, + 0.4436, 0.4225, 0.4004, 0.3787, 0.3562, 0.3324, 0.3089, 0.2839, 0.2584, 0.2321, + 0.2047, 0.1784, 0.1499, 0.1213, 0.0915, 0.0628, 0.0314, 0.0000, -0.0321, -0.0657, + -0.0977, -0.1324, -0.1673, -0.2036, -0.2387, -0.2768, -0.3150, -0.3538, -0.3936, -0.4327, + -0.4739, -0.5148, -0.5561, -0.6000, -0.6438, -0.6889, -0.7331, -0.7781, -0.8247, -0.8712, + -0.9177, -0.9677, -1.0142, -1.0631, -1.1143, -1.1686, -1.2169, -1.2680, -1.3223, -1.3752, + -1.4261, -1.4806, -1.5356, -1.5890, -1.6462, -1.7041, -1.7591, -1.8124, -1.8735, -1.9311, + -1.9891, -2.0459, -2.1048, -2.1653, -2.2248, -2.2855, -2.3466, -2.4079, -2.4668, -2.5263, + -2.5876, -2.6507, -2.7142, -2.7761, -2.8366, -2.8995, -2.9620, -3.0279, -3.0973, -3.1576, + -3.2238, -3.2890, -3.3554, -3.4215, -3.4805, -3.5518, -3.6133, -3.6812, -3.7473, -3.8140, + -3.8781, -3.9450, -4.0184, -4.0794, -4.1478, -4.2241, -4.2853, -4.3473, -4.4062, -4.4839, + -4.5539, -4.6202, -4.6794, -4.7478, -4.8309, -4.9048, -4.9669, -5.0294, -5.1194, -5.1732, + -5.2378, -5.3094, -5.3742, -5.4573, -5.5190, -5.5728, -5.6637, -5.7259, -5.7843, -5.8854, + -5.9553, -6.0054, -6.0656, -6.1707, -6.2241, -6.3139, -6.3393, -6.4356, -6.5153, -6.5758, + -6.6506, -6.7193, -6.7542, -6.8942, -6.9219, -6.9605, -7.1013, -7.1895, -7.1549, -7.2799, + -7.4119, -7.4608, -7.5256, -7.5879, -7.7598, -8.4120}, + {0.9999, 0.9998, 0.9998, 0.9998, 0.9998, 0.9998, 0.9997, 0.9997, 0.9997, 0.9997, + 0.9997, 0.9996, 0.9996, 0.9996, 0.9995, 0.9995, 0.9994, 0.9994, 0.9994, 0.9993, + 0.9993, 0.9992, 0.9991, 0.9991, 0.9990, 0.9989, 0.9988, 0.9988, 0.9988, 0.9986, + 0.9985, 0.9984, 0.9983, 0.9982, 0.9980, 0.9979, 0.9977, 0.9976, 0.9974, 0.9971, + 0.9969, 0.9968, 0.9965, 0.9962, 0.9960, 0.9957, 0.9953, 0.9950, 0.9947, 0.9941, + 0.9937, 0.9933, 0.9928, 0.9922, 0.9917, 0.9911, 0.9904, 0.9897, 0.9890, 0.9882, + 0.9874, 0.9863, 0.9855, 0.9843, 0.9832, 0.9819, 0.9806, 0.9792, 0.9777, 0.9760, + 0.9743, 0.9724, 0.9704, 0.9683, 0.9659, 0.9634, 0.9609, 0.9581, 0.9550, 0.9516, + 0.9481, 0.9446, 0.9406, 0.9363, 0.9317, 0.9270, 0.9218, 0.9160, 0.9103, 0.9038, + 0.8972, 0.8898, 0.8822, 0.8739, 0.8647, 0.8554, 0.8457, 0.8357, 0.8231, 0.8115, + 0.7984, 0.7854, 0.7704, 0.7556, 0.7391, 0.7210, 0.7038, 0.6840, 0.6633, 0.6408, + 0.6174, 0.5939, 0.5678, 0.5410, 0.5137, 0.4836, 0.4524, 0.4193, 0.3850, 0.3482, + 0.3132, 0.2733, 0.2315, 0.1891, 0.1435, 0.0980, 0.0493, 0.0000, -0.0510, -0.1052, + -0.1593, -0.2177, -0.2759, -0.3374, -0.4005, -0.4599, -0.5266, -0.5935, -0.6626, -0.7328, + -0.8051, -0.8757, -0.9498, -1.0271, -1.1019, -1.1816, -1.2642, -1.3459, -1.4295, -1.5077, + -1.5958, -1.6818, -1.7647, -1.8548, -1.9387, -2.0295, -2.1152, -2.2154, -2.3011, -2.3904, + -2.4820, -2.5786, -2.6730, -2.7652, -2.8616, -2.9546, -3.0526, -3.1445, -3.2445, -3.3416, + -3.4357, -3.5325, -3.6324, -3.7313, -3.8225, -3.9209, -4.0248, -4.1278, -4.2261, -4.3193, + -4.4220, -4.5262, -4.6214, -4.7242, -4.8234, -4.9245, -5.0298, -5.1250, -5.2232, -5.3267, + -5.4332, -5.5342, -5.6431, -5.7270, -5.8401, -5.9350, -6.0407, -6.1418, -6.2363, -6.3384, + -6.4536, -6.5429, -6.6582, -6.7433, -6.8438, -6.9478, -7.0789, -7.1894, -7.2714, -7.3815, + -7.4810, -7.5575, -7.6852, -7.8071, -7.8580, -7.9724, -8.1000, -8.2207, -8.2867, -8.4017, + -8.5287, -8.6347, -8.7082, -8.8319, -8.9448, -9.0355, -9.1885, -9.2095, -9.2863, -9.4186, + -9.5064, -9.6386, -9.7207, -9.8286, -9.9453, -10.0701, -10.1735, -10.3001, -10.2858, -10.5427, + -10.5982, -10.7361, -10.7042, -10.9212, -11.0097, -11.0469, -11.1155, -11.2812, -11.3472, -11.4988, + -11.5327, -11.6692, -11.9376, -11.8606, -12.1372, -13.2539}, + {1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, + 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, + 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, + 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, + 0.9999, 0.9998, 0.9998, 0.9998, 0.9998, 0.9997, 0.9997, 0.9997, 0.9997, 0.9996, + 0.9996, 0.9995, 0.9995, 0.9994, 0.9994, 0.9993, 0.9992, 0.9991, 0.9991, 0.9989, + 0.9988, 0.9986, 0.9985, 0.9983, 0.9981, 0.9980, 0.9977, 0.9974, 0.9971, 0.9968, + 0.9965, 0.9962, 0.9956, 0.9950, 0.9948, 0.9941, 0.9933, 0.9926, 0.9919, 0.9910, + 0.9899, 0.9889, 0.9877, 0.9863, 0.9845, 0.9829, 0.9811, 0.9791, 0.9769, 0.9741, + 0.9716, 0.9684, 0.9645, 0.9611, 0.9563, 0.9519, 0.9463, 0.9406, 0.9344, 0.9272, + 0.9197, 0.9107, 0.9016, 0.8903, 0.8791, 0.8653, 0.8523, 0.8357, 0.8179, 0.7988, + 0.7779, 0.7562, 0.7318, 0.7024, 0.6753, 0.6435, 0.6089, 0.5700, 0.5296, 0.4860, + 0.4366, 0.3855, 0.3301, 0.2735, 0.2114, 0.1443, 0.0682, 0.0000, -0.0715, -0.1604, + -0.2478, -0.3377, -0.4287, -0.5277, -0.6291, -0.7384, -0.8457, -0.9559, -1.0742, -1.1913, + -1.3110, -1.4238, -1.5594, -1.6854, -1.8093, -1.9414, -2.0763, -2.2160, -2.3611, -2.4876, + -2.6374, -2.7710, -2.9225, -3.0591, -3.2077, -3.3452, -3.4916, -3.6316, -3.7735, -3.9296, + -4.0682, -4.2334, -4.3607, -4.5270, -4.6807, -4.8108, -4.9753, -5.1212, -5.2631, -5.4042, + -5.5510, -5.7227, -5.8794, -6.0244, -6.1677, -6.3271, -6.4862, -6.6130, -6.7449, -6.9250, + -7.1232, -7.1736, -7.3628, -7.5596, -7.6906, -7.8129, -7.9817, -8.1440, -8.3016, -8.4797, + -8.5734, -8.7692, -8.9198, -9.0610, -9.1746, -9.3536, -9.5939, -9.6957, -9.8475, -9.9639, + -10.1730, -10.2427, -10.4573, -10.5413, -10.7303, -10.9339, -11.0215, -11.2047, -11.2894, -11.4572, + -11.6256, -11.7794, -11.8801, -12.1717, -12.2354, -12.3686, -12.6195, -12.6527, -12.8247, -12.9560, + -13.3265, -13.1667, -13.4274, -13.6064, -13.5515, -13.9501, -13.9926, -14.4049, -14.1653, -14.4348, + -14.7983, -14.7807, -15.2349, -15.3536, -15.3026, -15.2739, -15.7170, -16.2161, -15.9185, -15.9490, + -16.6258, -16.5568, -16.4318, -16.7999, -16.4101, -17.6393, -17.7643, -17.2644, -17.5973, -17.0403, + -17.7039, -18.0073, -18.1840, -18.3848, -18.6286, -20.7063}, + {1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, + 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, + 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, + 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, 0.9999, + 0.9999, 0.9998, 0.9998, 0.9998, 0.9998, 0.9997, 0.9997, 0.9997, 0.9997, 0.9996, + 0.9996, 0.9995, 0.9995, 0.9994, 0.9994, 0.9993, 0.9992, 0.9991, 0.9991, 0.9989, + 0.9988, 0.9986, 0.9985, 0.9983, 0.9981, 0.9980, 0.9977, 0.9974, 0.9971, 0.9968, + 0.9965, 0.9962, 0.9956, 0.9950, 0.9948, 0.9941, 0.9933, 0.9926, 0.9919, 0.9910, + 0.9899, 0.9889, 0.9877, 0.9863, 0.9845, 0.9829, 0.9811, 0.9791, 0.9769, 0.9741, + 0.9716, 0.9684, 0.9645, 0.9611, 0.9563, 0.9519, 0.9463, 0.9406, 0.9344, 0.9272, + 0.9197, 0.9107, 0.9016, 0.8903, 0.8791, 0.8653, 0.8523, 0.8357, 0.8179, 0.7988, + 0.7779, 0.7562, 0.7318, 0.7024, 0.6753, 0.6435, 0.6089, 0.5700, 0.5296, 0.4860, + 0.4366, 0.3855, 0.3301, 0.2735, 0.2114, 0.1443, 0.0682, 0.0000, -0.0715, -0.1604, + -0.2478, -0.3377, -0.4287, -0.5277, -0.6291, -0.7384, -0.8457, -0.9559, -1.0742, -1.1913, + -1.3110, -1.4238, -1.5594, -1.6854, -1.8093, -1.9414, -2.0763, -2.2160, -2.3611, -2.4876, + -2.6374, -2.7710, -2.9225, -3.0591, -3.2077, -3.3452, -3.4916, -3.6316, -3.7735, -3.9296, + -4.0682, -4.2334, -4.3607, -4.5270, -4.6807, -4.8108, -4.9753, -5.1212, -5.2631, -5.4042, + -5.5510, -5.7227, -5.8794, -6.0244, -6.1677, -6.3271, -6.4862, -6.6130, -6.7449, -6.9250, + -7.1232, -7.1736, -7.3628, -7.5596, -7.6906, -7.8129, -7.9817, -8.1440, -8.3016, -8.4797, + -8.5734, -8.7692, -8.9198, -9.0610, -9.1746, -9.3536, -9.5939, -9.6957, -9.8475, -9.9639, + -10.1730, -10.2427, -10.4573, -11.7794, -11.8801, -12.1717, -12.2354, -12.3686, -12.6195, -12.6527, + -12.8247, -12.9560, -13.3265, -13.1667, -13.4274, -13.6064, -13.5515, -13.9501, -13.9926, -14.4049, + -14.1653, -14.4348, -14.7983, -14.7807, -15.2349, -15.3536, -15.3026, -15.2739, -15.7170, -16.2161, + -15.9185, -15.9490, -16.6258, -16.5568, -16.4318, -16.7999, -16.4101, -17.6393, -17.7643, -17.2644, + -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}}; diff --git a/wspr/misc.c b/wspr/misc.c new file mode 100644 index 00000000..1708cf37 --- /dev/null +++ b/wspr/misc.c @@ -0,0 +1,45 @@ +#include "wspr.h" + +// used by qsort +// NB: assumes first element in struct is float sort field +int qsort_floatcomp(const void *elem1, const void *elem2) +{ + const float f1 = *(const float *) elem1, f2 = *(const float *) elem2; + if (f1 < f2) + return -1; + return f1 > f2; +} + +static bool init = false; +static u4_t epoch_sec; + +static void set_epoch() +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + epoch_sec = ts.tv_sec; + + init = true; +} + +// overflows 136 years after timer epoch +u4_t timer_sec() +{ + struct timespec ts; + + if (!init) set_epoch(); + clock_gettime(CLOCK_MONOTONIC, &ts); + return ts.tv_sec - epoch_sec; +} + +// overflows 49.7 days after timer epoch +u4_t timer_ms() +{ + struct timespec ts; + + if (!init) set_epoch(); + clock_gettime(CLOCK_MONOTONIC, &ts); + int msec = ts.tv_nsec/1000000; + assert(msec >= 0 && msec < 1000); + return (ts.tv_sec - epoch_sec)*1000 + msec; +} diff --git a/wspr/misc.h b/wspr/misc.h new file mode 100644 index 00000000..56c53416 --- /dev/null +++ b/wspr/misc.h @@ -0,0 +1,5 @@ +#pragma once + +int qsort_floatcomp(const void* elem1, const void* elem2); +u4_t timer_sec(); +u4_t timer_ms(); diff --git a/wspr/nhash.c b/wspr/nhash.c new file mode 100644 index 00000000..1498f4f2 --- /dev/null +++ b/wspr/nhash.c @@ -0,0 +1,368 @@ +/* + This file is part of wsprd. + + File name: nhash.c + + *------------------------------------------------------------------------------ + * + * This file is part of the WSPR application, Weak Signal Propogation Reporter + * + * File Name: nhash.c + * Description: Functions to produce 32-bit hashes for hash table lookup + * + * Copyright (C) 2008-2014 Joseph Taylor, K1JT + * License: GNU GPL v3+ + * + * This program is free software; you can redistribute it and/or modify it under + * the terms of the GNU General Public License as published by the Free Software + * Foundation; either version 3 of the License, or (at your option) any later + * version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * this program; if not, write to the Free Software Foundation, Inc., 51 Franklin + * Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * Files: lookup3.c + * Copyright: Copyright (C) 2006 Bob Jenkins + * License: public-domain + * You may use this code any way you wish, private, educational, or commercial. + * It's free. + * + *------------------------------------------------------------------------------- +*/ + +/* +These are functions for producing 32-bit hashes for hash table lookup. +hashword(), hashlittle(), hashlittle2(), hashbig(), mix(), and final() +are externally useful functions. Routines to test the hash are included +if SELF_TEST is defined. You can use this free for any purpose. It's in +the public domain. It has no warranty. + +You probably want to use hashlittle(). hashlittle() and hashbig() +hash byte arrays. hashlittle() is is faster than hashbig() on +little-endian machines. Intel and AMD are little-endian machines. +On second thought, you probably want hashlittle2(), which is identical to +hashlittle() except it returns two 32-bit hashes for the price of one. +You could implement hashbig2() if you wanted but I haven't bothered here. + +If you want to find a hash of, say, exactly 7 integers, do + a = i1; b = i2; c = i3; + mix(a,b,c); + a += i4; b += i5; c += i6; + mix(a,b,c); + a += i7; + final(a,b,c); +then use c as the hash value. If you have a variable length array of +4-byte integers to hash, use hashword(). If you have a byte array (like +a character string), use hashlittle(). If you have several byte arrays, or +a mix of things, see the comments above hashlittle(). + +Why is this so big? I read 12 bytes at a time into 3 4-byte integers, +then mix those integers. This is fast (you can do a lot more thorough +mixing with 12*3 instructions on 3 integers than you can with 3 instructions +on 1 byte), but shoehorning those bytes into integers efficiently is messy. +*/ + +#define SELF_TEST 1 + +#include /* defines printf for tests */ +#include /* defines time_t for timings in the test */ +#include "nhash.h" + +#define HASH_LITTLE_ENDIAN 1 + +#define hashsize(n) ((uint32_t)1<<(n)) +#define hashmask(n) (hashsize(n)-1) +#define rot(x,k) (((x)<<(k)) | ((x)>>(32-(k)))) + +/* +------------------------------------------------------------------------------- +mix -- mix 3 32-bit values reversibly. + +This is reversible, so any information in (a,b,c) before mix() is +still in (a,b,c) after mix(). + +If four pairs of (a,b,c) inputs are run through mix(), or through +mix() in reverse, there are at least 32 bits of the output that +are sometimes the same for one pair and different for another pair. +This was tested for: +* pairs that differed by one bit, by two bits, in any combination + of top bits of (a,b,c), or in any combination of bottom bits of + (a,b,c). +* "differ" is defined as +, -, ^, or ~^. For + and -, I transformed + the output delta to a Gray code (a^(a>>1)) so a string of 1's (as + is commonly produced by subtraction) look like a single 1-bit + difference. +* the base values were pseudorandom, all zero but one bit set, or + all zero plus a counter that starts at zero. + +Some k values for my "a-=c; a^=rot(c,k); c+=b;" arrangement that +satisfy this are + 4 6 8 16 19 4 + 9 15 3 18 27 15 + 14 9 3 7 17 3 +Well, "9 15 3 18 27 15" didn't quite get 32 bits diffing +for "differ" defined as + with a one-bit base and a two-bit delta. I +used http://burtleburtle.net/bob/hash/avalanche.html to choose +the operations, constants, and arrangements of the variables. + +This does not achieve avalanche. There are input bits of (a,b,c) +that fail to affect some output bits of (a,b,c), especially of a. The +most thoroughly mixed value is c, but it doesn't really even achieve +avalanche in c. + +This allows some parallelism. Read-after-writes are good at doubling +the number of bits affected, so the goal of mixing pulls in the opposite +direction as the goal of parallelism. I did what I could. Rotates +seem to cost as much as shifts on every machine I could lay my hands +on, and rotates are much kinder to the top and bottom bits, so I used +rotates. +------------------------------------------------------------------------------- +*/ +#define mix(a,b,c) \ +{ \ + a -= c; a ^= rot(c, 4); c += b; \ + b -= a; b ^= rot(a, 6); a += c; \ + c -= b; c ^= rot(b, 8); b += a; \ + a -= c; a ^= rot(c,16); c += b; \ + b -= a; b ^= rot(a,19); a += c; \ + c -= b; c ^= rot(b, 4); b += a; \ +} + +/* +------------------------------------------------------------------------------- +final -- final mixing of 3 32-bit values (a,b,c) into c + +Pairs of (a,b,c) values differing in only a few bits will usually +produce values of c that look totally different. This was tested for +* pairs that differed by one bit, by two bits, in any combination + of top bits of (a,b,c), or in any combination of bottom bits of + (a,b,c). +* "differ" is defined as +, -, ^, or ~^. For + and -, I transformed + the output delta to a Gray code (a^(a>>1)) so a string of 1's (as + is commonly produced by subtraction) look like a single 1-bit + difference. +* the base values were pseudorandom, all zero but one bit set, or + all zero plus a counter that starts at zero. + +These constants passed: + 14 11 25 16 4 14 24 + 12 14 25 16 4 14 24 +and these came close: + 4 8 15 26 3 22 24 + 10 8 15 26 3 22 24 + 11 8 15 26 3 22 24 +------------------------------------------------------------------------------- +*/ +#define final(a,b,c) \ +{ \ + c ^= b; c -= rot(b,14); \ + a ^= c; a -= rot(c,11); \ + b ^= a; b -= rot(a,25); \ + c ^= b; c -= rot(b,16); \ + a ^= c; a -= rot(c,4); \ + b ^= a; b -= rot(a,14); \ + c ^= b; c -= rot(b,24); \ +} + +/* +------------------------------------------------------------------------------- +hashlittle() -- hash a variable-length key into a 32-bit value + k : the key (the unaligned variable-length array of bytes) + length : the length of the key, counting by bytes + initval : can be any 4-byte value +Returns a 32-bit value. Every bit of the key affects every bit of +the return value. Two keys differing by one or two bits will have +totally different hash values. + +The best hash table sizes are powers of 2. There is no need to do +mod a prime (mod is sooo slow!). If you need less than 32 bits, +use a bitmask. For example, if you need only 10 bits, do + h = (h & hashmask(10)); +In which case, the hash table should have hashsize(10) elements. + +If you are hashing n strings (uint8_t **)k, do it like this: + for (i=0, h=0; i 12) + { + a += k[0]; + b += k[1]; + c += k[2]; + mix(a,b,c); + length -= 12; + k += 3; + } + + /*----------------------------- handle the last (probably partial) block */ + /* + * "k[2]&0xffffff" actually reads beyond the end of the string, but + * then masks off the part it's not allowed to read. Because the + * string is aligned, the masked-off tail is in the same word as the + * rest of the string. Every machine with memory protection I've seen + * does it on word boundaries, so is OK with this. But VALGRIND will + * still catch it and complain. The masking trick does make the hash + * noticably faster for short strings (like English words). + */ +#ifndef VALGRIND + + switch(length) + { + case 12: c+=k[2]; b+=k[1]; a+=k[0]; break; + case 11: c+=k[2]&0xffffff; b+=k[1]; a+=k[0]; break; + case 10: c+=k[2]&0xffff; b+=k[1]; a+=k[0]; break; + case 9 : c+=k[2]&0xff; b+=k[1]; a+=k[0]; break; + case 8 : b+=k[1]; a+=k[0]; break; + case 7 : b+=k[1]&0xffffff; a+=k[0]; break; + case 6 : b+=k[1]&0xffff; a+=k[0]; break; + case 5 : b+=k[1]&0xff; a+=k[0]; break; + case 4 : a+=k[0]; break; + case 3 : a+=k[0]&0xffffff; break; + case 2 : a+=k[0]&0xffff; break; + case 1 : a+=k[0]&0xff; break; + case 0 : return c; /* zero length strings require no mixing */ + } + +#else /* make valgrind happy */ + + k8 = (const uint8_t *)k; + switch(length) + { + case 12: c+=k[2]; b+=k[1]; a+=k[0]; break; + case 11: c+=((uint32_t)k8[10])<<16; /* fall through */ + case 10: c+=((uint32_t)k8[9])<<8; /* fall through */ + case 9 : c+=k8[8]; /* fall through */ + case 8 : b+=k[1]; a+=k[0]; break; + case 7 : b+=((uint32_t)k8[6])<<16; /* fall through */ + case 6 : b+=((uint32_t)k8[5])<<8; /* fall through */ + case 5 : b+=k8[4]; /* fall through */ + case 4 : a+=k[0]; break; + case 3 : a+=((uint32_t)k8[2])<<16; /* fall through */ + case 2 : a+=((uint32_t)k8[1])<<8; /* fall through */ + case 1 : a+=k8[0]; break; + case 0 : return c; + } + +#endif /* !valgrind */ + + } else if (HASH_LITTLE_ENDIAN && ((u.i & 0x1) == 0)) { + const uint16_t *k = (const uint16_t *)key; /* read 16-bit chunks */ + const uint8_t *k8; + + /*--------------- all but last block: aligned reads and different mixing */ + while (length > 12) + { + a += k[0] + (((uint32_t)k[1])<<16); + b += k[2] + (((uint32_t)k[3])<<16); + c += k[4] + (((uint32_t)k[5])<<16); + mix(a,b,c); + length -= 12; + k += 6; + } + + /*----------------------------- handle the last (probably partial) block */ + k8 = (const uint8_t *)k; + switch(length) + { + case 12: c+=k[4]+(((uint32_t)k[5])<<16); + b+=k[2]+(((uint32_t)k[3])<<16); + a+=k[0]+(((uint32_t)k[1])<<16); + break; + case 11: c+=((uint32_t)k8[10])<<16; /* fall through */ + case 10: c+=k[4]; + b+=k[2]+(((uint32_t)k[3])<<16); + a+=k[0]+(((uint32_t)k[1])<<16); + break; + case 9 : c+=k8[8]; /* fall through */ + case 8 : b+=k[2]+(((uint32_t)k[3])<<16); + a+=k[0]+(((uint32_t)k[1])<<16); + break; + case 7 : b+=((uint32_t)k8[6])<<16; /* fall through */ + case 6 : b+=k[2]; + a+=k[0]+(((uint32_t)k[1])<<16); + break; + case 5 : b+=k8[4]; /* fall through */ + case 4 : a+=k[0]+(((uint32_t)k[1])<<16); + break; + case 3 : a+=((uint32_t)k8[2])<<16; /* fall through */ + case 2 : a+=k[0]; + break; + case 1 : a+=k8[0]; + break; + case 0 : return c; /* zero length requires no mixing */ + } + + } else { /* need to read the key one byte at a time */ + const uint8_t *k = (const uint8_t *)key; + + /*--------------- all but the last block: affect some 32 bits of (a,b,c) */ + while (length > 12) + { + a += k[0]; + a += ((uint32_t)k[1])<<8; + a += ((uint32_t)k[2])<<16; + a += ((uint32_t)k[3])<<24; + b += k[4]; + b += ((uint32_t)k[5])<<8; + b += ((uint32_t)k[6])<<16; + b += ((uint32_t)k[7])<<24; + c += k[8]; + c += ((uint32_t)k[9])<<8; + c += ((uint32_t)k[10])<<16; + c += ((uint32_t)k[11])<<24; + mix(a,b,c); + length -= 12; + k += 12; + } + + /*-------------------------------- last block: affect all 32 bits of (c) */ + switch(length) /* all the case statements fall through */ + { + case 12: c+=((uint32_t)k[11])<<24; + case 11: c+=((uint32_t)k[10])<<16; + case 10: c+=((uint32_t)k[9])<<8; + case 9 : c+=k[8]; + case 8 : b+=((uint32_t)k[7])<<24; + case 7 : b+=((uint32_t)k[6])<<16; + case 6 : b+=((uint32_t)k[5])<<8; + case 5 : b+=k[4]; + case 4 : a+=((uint32_t)k[3])<<24; + case 3 : a+=((uint32_t)k[2])<<16; + case 2 : a+=((uint32_t)k[1])<<8; + case 1 : a+=k[0]; + break; + case 0 : return c; + } + } + + final(a,b,c); + c=(32767&c); + + return (u2_t) c; +} diff --git a/wspr/nhash.h b/wspr/nhash.h new file mode 100644 index 00000000..9f9dd55b --- /dev/null +++ b/wspr/nhash.h @@ -0,0 +1,16 @@ +#pragma once + +#include "types.h" + +#include +#include /* defines uint32_t etc */ + +#ifdef __cplusplus +extern "C" { +#endif + +u2_t nhash(const void *key, size_t length, uint32_t initval); + +#ifdef __cplusplus +} +#endif diff --git a/wspr/tab.c b/wspr/tab.c new file mode 100644 index 00000000..96cf177f --- /dev/null +++ b/wspr/tab.c @@ -0,0 +1,42 @@ +#include "wspr.h" + +/* + This file is part of wsprd. + + File name: tab.c + Description: 8-bit parity lookup table. +*/ +unsigned char Partab[] = { + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, +}; diff --git a/wspr/types.h b/wspr/types.h new file mode 100644 index 00000000..b2613ffb --- /dev/null +++ b/wspr/types.h @@ -0,0 +1,88 @@ +#pragma once + +#include + +typedef unsigned long long u64_t; +typedef unsigned int u4_t; +typedef unsigned char u1_t; +typedef unsigned short u2_t; + +typedef signed long long s64_t; +typedef signed int s4_t; +typedef signed short s2_t; +typedef signed char s1_t; + +typedef void (*func_t)(); +typedef void (*funcP_t)(void *); +typedef int (*funcPR_t)(void *); + +#define TO_VOID_PARAM(p) ((void *) (long) (p)) +#define FROM_VOID_PARAM(p) ((long) (p)) + +#define U1(v) ((u1_t) (v)) +#define S1(v) ((s1_t) (v)) +#define U2(v) ((u2_t) (v)) +#define S2(v) ((s2_t) (v)) +#define U4(v) ((u4_t) (v)) +#define S4(v) ((s4_t) (v)) +#define U8(v) ((u64_t) (v)) +#define S8(v) ((s64_t) (v)) + +#define S16x4_S64(a,b,c,d) S8( (U8(a)<<48) | (U8(b)<<32) | (U8(c)<<16) | U8(d) ) +#define S14_16(w) S2( U2(w) | ( ((U2(w)) & 0x2000)? 0xc000:0 ) ) +#define S14_32(w) S4( S2( U2(w) | ((U2(w) & 0x2000)? 0xffffc000:0 ) ) ) +#define S24_8_16(h8,l16) S4( (U1(h8)<<16) | U2(l16) | ((U1(h8) & 0x80)? 0xff000000:0) ) +#define S24_16_8(h16,l8) S4( (U2(h16)<<8) | U1(l8) | ((U2(h16) & 0x8000)? 0xff000000:0) ) + +#define B2I(bytes) (((bytes)+3)/4) +#define I2B(ints) ((ints)*4) +#define B2S(bytes) (((bytes)+1)/2) +#define S2B(shorts) ((shorts)*2) + +#define B3(i) (((i) >> 24) & 0xff) +#define B2(i) (((i) >> 16) & 0xff) +#define B1(i) (((i) >> 8) & 0xff) +#define B0(i) (((i) >> 0) & 0xff) + +#define FLIP32(i) ((B0(i) << 24) | (B1(i) << 16) | (B2(i) << 8) | (B3(i) << 0)) +#define FLIP16(i) ((B0(i) << 8) | (B1(i) << 0)) + +#ifndef TRUE + #define TRUE 1 +#endif + +#ifndef FALSE + #define FALSE 0 +#endif + +#ifndef __cplusplus + typedef unsigned char bool; + #define true TRUE + #define false FALSE +#endif + +#define NOT_FOUND -1 + +#define ARRAY_LEN(x) ((int) (sizeof (x) / sizeof ((x) [0]))) + +#define K 1024 +#define M (K*K) +#define B (M*K) + +#define MHz 1000000 +#define kHz 1000 + +#define MAX(a,b) ((a)>(b)?(a):(b)) +#define max(a,b) MAX(a,b) +#define MIN(a,b) ((a)<(b)?(a):(b)) +#define min(a,b) MIN(a,b) + +#define SI_CLAMP(a,n) ( ((a) > ((n)-1))? ((n)-1) : ( ((a) < -(n))? -(n) : (a) ) ) + +#define STRINGIFY(x) #x +#define STRINGIFY_DEFINE(x) STRINGIFY(x) // indirection needed for a -Dx=y define +#define CAT_STRING(x,y) x y // just a reminder of how this is done: "foo" "bar" +#define CAT_DEFINE_VAR(x,y) x ## y // just a reminder of how this is done: foo ## bar + +// documentation assistance +#define SPACE_FOR_NULL 1 diff --git a/wspr/wspr.c b/wspr/wspr.c new file mode 100644 index 00000000..4ded00c0 --- /dev/null +++ b/wspr/wspr.c @@ -0,0 +1,1076 @@ +/* + This file is part of program wsprd, a detector/demodulator/decoder + for the Weak Signal Propagation Reporter (WSPR) mode. + + Copyright 2001-2015, Joe Taylor, K1JT + + Much of the present code is based on work by Steven Franke, K9AN, + which in turn was based on earlier work by K1JT. + + Copyright 2014-2015, Steven Franke, K9AN + + License: GNU GPL v3 + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + */ + +#include "wspr.h" + +static const unsigned char pr3[NSYM_162]= +{1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0, + 0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1, + 0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1, + 1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1, + 0,0,1,0,1,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0, + 0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1, + 0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1, + 0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,0,1,1,0, + 0,0}; + +void sync_and_demodulate( + WSPR_CPX_t *id, WSPR_CPX_t *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. * + ************************************************************************/ + + static float fplast=-10000.0; + static float dt=1.0/FSRATE, df=FSRATE/FSPS; + static float pi=K_PI; + float twopidt, df15=df*1.5, df05=df*0.5; + + int i, j, k, lag; + float + i0[NSYM_162], q0[NSYM_162], + i1[NSYM_162], q1[NSYM_162], + i2[NSYM_162], q2[NSYM_162], + i3[NSYM_162], q3[NSYM_162]; + double p0, p1, p2, p3, cmet, totp, syncmax, fac; + double + c0[SPS], s0[SPS], + c1[SPS], s1[SPS], + c2[SPS], s2[SPS], + c3[SPS], s3[SPS]; + double + dphi0, cdphi0, sdphi0, + dphi1, cdphi1, sdphi1, + dphi2, cdphi2, sdphi2, + dphi3, cdphi3, sdphi3; + float f0=0.0, fp, ss, fbest=0.0, fsum=0.0, f2sum=0.0, fsymb[NSYM_162]; + int best_shift = 0, ifreq; + + syncmax=-1e30; + + if( mode == 0 ) { ifmin=0; ifmax=0; fstep=0.0; f0=*f1; } + if( mode == 1 ) { lagmin=*shift1; lagmax=*shift1; f0=*f1; } + if( mode == 2 ) { lagmin=*shift1; lagmax=*shift1; ifmin=0; ifmax=0; f0=*f1; } + + twopidt=2*pi*dt; + for(ifreq=ifmin; ifreq<=ifmax; ifreq++) { + f0=*f1+ifreq*fstep; + for(lag=lagmin; lag<=lagmax; lag=lag+lagstep) { + ss=0.0; + totp=0.0; + for (i=0; i0) && (k syncmax ) { //Save best parameters + syncmax=ss; + best_shift=lag; + fbest=f0; + } + } // lag loop + } //freq loop + + if( mode <=1 ) { //Send best params back to caller + *sync=syncmax; + *shift1=best_shift; + *f1=fbest; + return; + } + + if( mode == 2 ) { + *sync=syncmax; + for (i=0; i 127) fsymb[i]=127.0; + if( fsymb[i] < -128 ) fsymb[i]=-128.0; + symbols[i] = fsymb[i] + 128; + } + TRY_YIELD_DECODE; + return; + } + TRY_YIELD_DECODE; + return; +} + +// uses TRY_YIELD instead of TRY_YIELD_DECODE since called from non-decoder FFT code +void renormalize(wspr_t *w, float psavg[], float smspec[]) +{ + int i,j,k; + + // smooth with 7-point window and limit the spectrum to +/-150 Hz + int window[7] = {1,1,1,1,1,1,1}; + + for (i=0; imin_snr = pow(10.0,-7.0/10.0); //this is min snr in wspr bw + w->snr_scaling_factor = (w->wspr_type == WSPR_TYPE_2MIN)? 26.3 : 35.3; + for (j=0; j<411; j++) { + smspec[j]=smspec[j]/noise_level - 1.0; + if( smspec[j] < w->min_snr) smspec[j]=0.1*w->min_snr; + } + TRY_YIELD; +} + +/*************************************************************************** + 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 dt=1.0/FSRATE, df=FSRATE/FSPS; + int i, j, k; + float pi=K_PI, twopidt, fp; + + float i0,q0; + float c0[SPS],s0[SPS]; + float dphi, cdphi, sdphi; + + twopidt=2*pi*dt; + + for (i=0; i0) & (k0) & (k0) && (k(nsig-1-nfilt/2) ) { + norm=partialsum[nfilt/2+nsig-1-i]; + } else { + norm=1.0; + } + k=shift0+i; + j=i+nfilt; + if( (k>0) && (kdecode_ping_pong); + + WSPR_CPX_t *idat = w->i_data[w->decode_ping_pong]; + WSPR_CPX_t *qdat = w->q_data[w->decode_ping_pong]; + + static bool wspr_decode_init; + if (!wspr_decode_init) { + if (w->stackdecoder) + w->stack = (struct snode *) malloc(WSPR_STACKSIZE * sizeof(struct snode)); + + wspr_decode_init = true; + } + + int uniques = 0; + + // multi-pass strategies + //#define SUBTRACT_SIGNAL // FIXME: how to implement spectrum subtraction given our incrementally-computed FFTs? + #define MORE_EFFORT // this scheme repeats work as maxcycles is increased, but it's difficult to eliminate that + + #if defined(SUBTRACT_SIGNAL) + npasses = 2; + #elif defined(MORE_EFFORT) + npasses = 0; // unlimited + #endif + + for (ipass=0; (npasses == 0 || ipass < npasses) && !w->abort_decode; ipass++) { + + #if defined(SUBTRACT_SIGNAL) + if (ipass > 0 && ndecodes_pass == 0) break; + #elif defined(MORE_EFFORT) + if (ipass == 0) { + maxcycles = 200; + iifac = 2; + } else { + if (maxcycles < 10000) maxcycles *= 2; + iifac = 1; + } + #endif + + // only build the peak list on the first pass + if (ipass == 0) { + float smspec[nbins_411]; + renormalize(w, w->pwr_sampavg[w->decode_ping_pong], smspec); + + // Find all local maxima in smoothed spectrum. + memset(pk, 0, sizeof(pk)); + + npk = 0; + bool candidate; + + if (w->more_candidates) { + for (j=0; j w->min_snr) && (npk < NPK); + if (candidate) { + pk[npk].bin0 = j; + pk[npk].freq0 = (j-hbins_205)*df; + pk[npk].snr0 = 10*log10(smspec[j]) - w->snr_scaling_factor; + npk++; + } + } + } else { + for (j=1; j<(nbins_411-1); j++) { + candidate = (smspec[j]>smspec[j-1]) && + (smspec[j]>smspec[j+1]) && + (npksnr_scaling_factor; + npk++; + } + } + } + //wdprintf("initial npk %d/%d\n", npk, NPK); + TRY_YIELD_DECODE; + + // Don't waste time on signals outside of the range [fmin,fmax]. + i=0; + for (pki=0; pki < npk; pki++) { + if (pk[pki].freq0 >= FMIN && pk[pki].freq0 <= FMAX) { + pk[i].ignore = false; + pk[i].bin0 = pk[pki].bin0; + pk[i].freq0 = pk[pki].freq0; + pk[i].snr0 = pk[pki].snr0; + i++; + } + } + npk = i; + //wdprintf("freq range limited npk %d\n", npk); + TRY_YIELD_DECODE; + + // only look at a limited number of strong peaks + qsort(pk, npk, sizeof(pk_t), snr_comp); // sort in decreasing snr order + if (npk > MAX_NPK) npk = MAX_NPK; + //wdprintf("MAX_NPK limited npk %d/%d\n", npk, MAX_NPK); + + // remember our frequency-sorted index + qsort(pk, npk, sizeof(pk_t), freq_comp); + for (pki=0; pki < npk; pki++) { + pk[pki].freq_idx = pki; + } + + // send peak info to client in increasing frequency order + memcpy(pk_freq, pk, sizeof(pk)); + + for (pki=0; pki < npk; pki++) { + sprintf(pk_freq[pki].snr_call, "%d", (int) roundf(pk_freq[pki].snr0)); + } + + wspr_send_peaks(w, pk_freq, npk); + + // keep 'pk' in snr order so strong sigs are processed first in case we run out of decode time + qsort(pk, npk, sizeof(pk_t), snr_comp); // sort in decreasing snr order + } + + int valid_peaks = 0; + for (pki=0; pki < npk; pki++) { + if (pk[pki].ignore) continue; + valid_peaks++; + } + wdprintf("PASS %d npeaks=%d valid_peaks=%d maxcycles=%d jig_range=%+d..%d jig_step=%d ---------------------------------------------\n", + ipass+1, npk, valid_peaks, maxcycles, jig_range/2, -jig_range/2, iifac); + + /* Make coarse estimates of shift (DT), freq, and drift + + * Look for time offsets up to +/- 8 symbols (about +/- 5.4 s) relative + to nominal start time, which is 2 seconds into the file + + * Calculates shift relative to the beginning of the file + + * Negative shifts mean that signal started before start of file + + * The program prints DT = shift-2 s + + * Shifts that cause sync vector to fall off of either end of the data + vector are accommodated by "partial decoding", such that missing + symbols produce a soft-decision symbol value of 128 + + * The frequency drift model is linear, deviation of +/- drift/2 over the + span of 162 symbols, with deviation equal to 0 at the center of the + signal vector. + */ + + int idrift,ifr,if0,ifd,k0; + int kindex; + float smax,ss,power,p0,p1,p2,p3; + + for (pki=0; pki < npk; pki++) { //For each candidate... + pk_t *p = &pk[pki]; + smax = -1e30; + if0 = p->freq0/df+SPS; + + #if defined(MORE_EFFORT) + if (ipass != 0 && p->ignore) + continue; + #endif + + for (ifr=if0-2; ifr<=if0+2; ifr++) { //Freq search + for( k0=-10; k0<22; k0++) { //Time search + for (idrift=-maxdrift; idrift<=maxdrift; idrift++) { //Drift search + ss=0.0; + power=0.0; + for (k=0; kpwr_samp[w->decode_ping_pong][ifd-3][kindex]; + p1=w->pwr_samp[w->decode_ping_pong][ifd-1][kindex]; + p2=w->pwr_samp[w->decode_ping_pong][ifd+1][kindex]; + p3=w->pwr_samp[w->decode_ping_pong][ifd+3][kindex]; + + p0=sqrt(p0); + p1=sqrt(p1); + p2=sqrt(p2); + p3=sqrt(p3); + + ss=ss+(2*pr3[k]-1)*((p1+p3)-(p0+p2)); + power=power+p0+p1+p2+p3; + } + } + sync1=ss/power; + if( sync1 > smax ) { //Save coarse parameters + smax=sync1; + p->shift0=HSPS*(k0+1); + p->drift0=idrift; + p->freq0=(ifr-SPS)*df; + p->sync0=sync1; + } + //wdprintf("drift %d k0 %d sync %f\n",idrift,k0,smax); + } + TRY_YIELD_DECODE; + } + } + wdprintf("npeak #%02d %6.1f snr %9.6f (%7.2f) freq %4.1f drift %5d shift %6.3f sync %3d bin\n", + pki, p->snr0, w->dialfreq_MHz+(bfo+p->freq0)/1e6, w->cf_offset+p->freq0, p->drift0, p->shift0, p->sync0, p->bin0); + } + + /* + Refine the estimates of freq, shift using sync as a metric. + Sync is calculated such that it is a float taking values in the range + [0.0,1.0]. + + Function sync_and_demodulate has three modes of operation + mode is the last argument: + + 0 = no frequency or drift search. find best time lag. + 1 = no time lag or drift search. find best frequency. + 2 = no frequency or time lag search. Calculate soft-decision + symbols using passed frequency and shift. + + NB: best possibility for OpenMP may be here: several worker threads + could each work on one candidate at a time. + */ + int candidates = 0; + ndecodes_pass = 0; + + for (pki=0; pki < npk && !w->abort_decode; pki++) { + bool f_decoded = false, f_delete = false, f_image = false, f_decoding = false; + + pk_t *p = &pk[pki]; + u4_t decode_start = timer_ms(); + + #if defined(MORE_EFFORT) + if (ipass != 0 && p->ignore) + continue; + #endif + + candidates++; + f1 = p->freq0; + snr = p->snr0; + drift1 = p->drift0; + shift1 = p->shift0; + sync1 = p->sync0; + + f_decoding = true; + pk_freq[p->freq_idx].flags |= WSPR_F_DECODING; + wspr_send_peaks(w, pk_freq, npk); + + wdprintf("start #%02d %6.1f snr %9.6f (%7.2f) freq %4.1f drift %5d shift %6.3f sync\n", + pki, snr, w->dialfreq_MHz+(bfo+f1)/1e6, w->cf_offset+f1, drift1, shift1, sync1); + + // coarse-grid lag and freq search, then if sync > minsync1 continue + fstep=0.0; ifmin=0; ifmax=0; + lagmin = shift1-128; + lagmax = shift1+128; + lagstep = 64; + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, &shift1, + lagmin, lagmax, lagstep, drift1, symfac, &sync1, 0); + + fstep = 0.25; ifmin = -2; ifmax = 2; + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, &shift1, + lagmin, lagmax, lagstep, drift1, symfac, &sync1, 1); + + // refine drift estimate + fstep=0.0; ifmin=0; ifmax=0; + float driftp,driftm,syncp,syncm; + driftp = drift1+0.5; + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, &shift1, + lagmin, lagmax, lagstep, driftp, symfac, &syncp, 1); + + driftm = drift1-0.5; + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, &shift1, + lagmin, lagmax, lagstep, driftm, symfac, &syncm, 1); + + if (syncp > sync1) { + drift1 = driftp; + sync1 = syncp; + } else + + if (syncm > sync1) { + drift1 = driftm; + sync1 = syncm; + } + + wdprintf("coarse #%02d %6.1f snr %9.6f (%7.2f) freq %4.1f drift %5d shift %6.3f sync\n", + pki, snr, w->dialfreq_MHz+(bfo+f1)/1e6, w->cf_offset+f1, drift1, shift1, sync1); + + // fine-grid lag and freq search + bool r_minsync1 = (sync1 > minsync1); + + if (r_minsync1) { + lagmin = shift1-32; lagmax = shift1+32; lagstep = 16; + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, &shift1, + lagmin, lagmax, lagstep, drift1, symfac, &sync1, 0); + + // fine search over frequency + fstep = 0.05; ifmin = -2; ifmax = 2; + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, &shift1, + lagmin, lagmax, lagstep, drift1, symfac, &sync1, 1); + } else { + p->ignore = true; + wdprintf("MINSYNC1 #%02d\n", pki); + f_delete = true; + } + + int idt=0, ii=0, jiggered_shift; + float y, sq, rms; + int r_decoded = 0; + bool r_tooWeak = true; + + // ii: 0 +1 -1 +2 -2 +3 -3 ... (*iifac) + // ii always covers jig_range, stepped by iifac resolution + while (!w->abort_decode && r_minsync1 && !r_decoded && idt <= (jig_range/iifac)) { + ii = (idt+1)/2; + if ((idt&1) == 1) ii = -ii; + ii = iifac*ii; + jiggered_shift = shift1+ii; + + // Use mode 2 to get soft-decision symbols + sync_and_demodulate(idat, qdat, npoints, w->symbols, &f1, ifmin, ifmax, fstep, + &jiggered_shift, lagmin, lagmax, lagstep, drift1, symfac, + &sync1, 2); + + sq = 0.0; + for (i=0; isymbols[i] - 128.0; + sq += y*y; + } + rms = sqrt(sq/FNSYM_162); + + bool weak = true; + if ((sync1 > minsync2) && (rms > minrms)) { + deinterleave(w->symbols); + + if (w->stackdecoder) { + r_decoded = jelinek(&metric, &cycles, w->decdata, w->symbols, NBITS, + WSPR_STACKSIZE, w->stack, mettab, maxcycles); + } else { + r_decoded = fano(&metric, &cycles, &maxnp, w->decdata, w->symbols, NBITS, + mettab, delta, maxcycles); + } + + r_tooWeak = weak = false; + } + + wdprintf("jig <>%3d #%02d %6.1f snr %9.6f (%7.2f) freq %4.1f drift %5d(%+4d) shift %6.3f sync %4.1f rms", + idt, pki, snr, w->dialfreq_MHz+(bfo+f1)/1e6, w->cf_offset+f1, drift1, jiggered_shift, ii, sync1, rms); + if (!weak) { + wprintf(" %4u metric %3u cycles\n", metric, cycles); + } else { + if (sync1 <= minsync2) wprintf(" SYNC-WEAK"); + if (rms <= minrms) wprintf(" RMS-WEAK"); + wprintf("\n"); + } + + idt++; + if (w->quickmode) break; + } + + bool r_timeUp = (!w->abort_decode && r_minsync1 && !r_decoded && idt > (jig_range/iifac)); + int r_valid = 0; + + //if (r_timeUp && r_tooWeak && iifac == 1) { + if (r_timeUp && r_tooWeak) { + wdprintf("NO CHANGE #%02d\n", pki); + p->ignore = true; // situation not going to get any better + f_delete = true; + } + + // result priority: abort, minSync1, tooWeak, noChange, timeUp, image, decoded + + if (r_decoded) { + ndecodes_pass++; + p->ignore = true; + + // Unpack the decoded message, update the hashtable, apply + // sanity checks on grid and power, and return + // call_loc_pow string and also callsign (for de-duping). + int dBm; + r_valid = unpk_(w->decdata, w->call_loc_pow, w->callsign, w->grid, &dBm); + + // subtract even on last pass + #ifdef SUBTRACT_SIGNAL + if (w->subtraction && (ipass < npasses) && r_valid > 0) { + if (get_wspr_channel_symbols(w->call_loc_pow, w->channel_symbols)) { + subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, w->channel_symbols); + } else { + break; + } + + } + #endif + + // Remove dupes (same callsign and freq within 3 Hz) and images + bool r_dupe = false; + if (r_valid > 0) { + bool hash = (strcmp(w->callsign, "...") == 0); + for (i=0; i < uniques; i++) { + decode_t *dp = &w->deco[i]; + bool match = (strcmp(w->callsign, dp->call) == 0); + bool close = (fabs(f1 - dp->freq) < 3.0); + if ((match && close) || (match && !hash)) { + if (close) { + f_delete = true; + } else { + f_image = true; + } + wdprintf("%s #%02d with #%02d %s, %.3f secs\n", f_image? "IMAGE" : "DUPE ", + pki, i, dp->call, (float)(timer_ms()-decode_start)/1e3); + r_dupe = true; + break; + } + } + } + + if (r_valid <= 0 && !r_dupe) { + if (r_valid < 0) { + wdprintf("UNPK ERR! #%02d error code %d, %.3f secs\n", + pki, r_valid, (float)(timer_ms()-decode_start)/1e3); + } else { + wdprintf("NOT VALID #%02d %.3f secs\n", + pki, (float)(timer_ms()-decode_start)/1e3); + } + f_delete = true; + } + + if (r_valid > 0 && !r_dupe) { + f_decoded = true; + + if (w->wspr_type == WSPR_TYPE_15MIN) { + freq_print = w->dialfreq_MHz + (bfo+112.5+f1/8.0)/1e6; + dt_print = shift1*8*dt-1.0; + } else { + freq_print = w->dialfreq_MHz + (bfo+f1)/1e6; + dt_print = shift1*dt-1.0; + } + + struct tm tm; + gmtime_r(&w->utc[w->decode_ping_pong], &tm); + + wdprintf("TYPE%d %02d%02d %3.0f %4.1f %10.6f %2d %-s %4s %2d [%s] in %.3f secs --------------------------------------------------------------------\n", + r_valid, tm.tm_hour, tm.tm_min, snr, dt_print, freq_print, (int) drift1, + w->callsign, w->grid, dBm, w->call_loc_pow, (float)(timer_ms()-decode_start)/1e3); + + double watts, factor; + char *W_s; + const char *units; + + watts = pow(10.0, (dBm - 30)/10.0); + if (watts >= 1.0) { + factor = 1.0; + units = "W"; + } else + if (watts >= 1e-3) { + factor = 1e3; + units = "mW"; + } else + if (watts >= 1e-6) { + factor = 1e6; + units = "uW"; + } else + if (watts >= 1e-9) { + factor = 1e9; + units = "nW"; + } else { + factor = 1e12; + units = "pW"; + } + + watts *= factor; + if (watts < 10.0) + asprintf(&W_s, "%.1f %s", watts, units); + else + asprintf(&W_s, "%.0f %s", watts, units); + + // Call Grid km dBm + // TYPE1 c6cccc g4gg kkkkk ppp (s) + // TYPE2 ... ppp (s) ; no hash yet + // TYPE2 ppp/c6cccc ppp (s) ; 1-3 char prefix + // TYPE2 c6cccc/ss ppp (s) ; 1-2 char suffix + // TYPE3 ... g6gggg kkkkk ppp (s) ; no hash yet + // TYPE3 ppp/c6cccc g6gggg kkkkk ppp (s) ; worst case, just let the fields float + + if (r_valid == 1) // TYPE1 + fprintf(w->demod_pipe, + "%02d%02d %3.0f %4.1f %9.6f %2d " + "%-6s " + "%s " + "%5d %3d (%s)\n", + tm.tm_hour, tm.tm_min, snr, dt_print, freq_print, (int) drift1, + w->callsign, w->callsign, w->grid, w->grid, + (int) grid_to_distance_km(w->grid), dBm, W_s); + else + + if (r_valid == 2) { // TYPE2 + if (strcmp(w->callsign, "...") == 0) + fprintf(w->demod_pipe, + "%02d%02d %3.0f %4.1f %9.6f %2d ... %3d (%s)\n", + tm.tm_hour, tm.tm_min, snr, dt_print, freq_print, (int) drift1, dBm, W_s); + else + fprintf(w->demod_pipe, + "%02d%02d %3.0f %4.1f %9.6f %2d " + "%-17s" + " %3d (%s)\n", + //kkkkk--ppp + tm.tm_hour, tm.tm_min, snr, dt_print, freq_print, (int) drift1, + w->callsign, w->callsign, dBm, W_s); + } else + + if (r_valid == 3) { // TYPE3 + if (strcmp(w->callsign, "...") == 0) + fprintf(w->demod_pipe, + "%02d%02d %3.0f %4.1f %9.6f %2d " + "... " + "%6s " + "%5d %3d (%s)\n", + tm.tm_hour, tm.tm_min, snr, dt_print, freq_print, (int) drift1, + w->grid, w->grid, (int) grid_to_distance_km(w->grid), dBm, W_s); + else + fprintf(w->demod_pipe, + "%02d%02d %3.0f %4.1f %9.6f %2d " + "%s " + "%s " + "%d %d (%s)\n", + tm.tm_hour, tm.tm_min, snr, dt_print, freq_print, (int) drift1, + w->callsign, w->callsign, w->grid, w->grid, + (int) grid_to_distance_km(w->grid), dBm, W_s); + } + + free(W_s); + fflush(w->demod_pipe); + + decode_t *dp = &w->deco[uniques]; + strcpy(dp->call, w->callsign); + dp->freq = f1; + dp->hour = tm.tm_hour; + dp->min = tm.tm_min; + dp->snr = snr; + dp->dt_print = dt_print; + dp->freq_print = freq_print; + dp->drift1 = drift1; + strcpy(dp->c_l_p, w->call_loc_pow); + uniques++; + } else { + r_valid = 0; + } + } else { + if (r_timeUp) { + wdprintf("TIME UP #%02d %.3f secs\n", pki, (float)(timer_ms()-decode_start)/1e3); + #if defined(MORE_EFFORT) + f_decoding = false; + #else + f_delete = true; + #endif + } + } // decoded + + if (f_decoded) { + strcpy(pk_freq[p->freq_idx].snr_call, w->callsign); + pk_freq[p->freq_idx].flags |= WSPR_F_DECODED; + } else + if (f_delete) { + pk_freq[p->freq_idx].flags |= WSPR_F_DELETE; + } else + if (f_image) { + strcpy(pk_freq[p->freq_idx].snr_call, "image"); + pk_freq[p->freq_idx].flags |= WSPR_F_IMAGE; + } else + if (!f_decoding) { + pk_freq[p->freq_idx].flags &= ~WSPR_F_DECODING; + } + + wspr_send_peaks(w, pk_freq, npk); + } // peak list + + if (candidates == 0) + break; // nothing left to do + + } // passes + + // when finished delete any unresolved peaks + for (i=0; i < npk; i++) { + pk_t *p = &pk[i]; + if (!(pk_freq[p->freq_idx].flags & WSPR_F_DECODED)) + pk_freq[p->freq_idx].flags |= WSPR_F_DELETE; + } + wspr_send_peaks(w, pk_freq, npk); + + // upload spots at the end of the decoding when there is less load on wsprnet.org + // FIXME OpenWebRX doesn't upload yet + for (i = 0; i < uniques; i++) { + decode_t *dp = &w->deco[i]; + ext_send_msg_encoded(w->rx_chan, WSPR_DEBUG_MSG, "EXT", "WSPR_UPLOAD", + "%02d%02d %3.0f %4.1f %9.6f %2d %s", + dp->hour, dp->min, dp->snr, dp->dt_print, dp->freq_print, (int) dp->drift1, dp->c_l_p); + wdprintf("WSPR_UPLOAD U%d/%d " + "%02d%02d %3.0f %4.1f %9.6f %2d %s" "\n", i, uniques, + dp->hour, dp->min, dp->snr, dp->dt_print, dp->freq_print, (int) dp->drift1, dp->c_l_p); + } +} diff --git a/wspr/wspr.h b/wspr/wspr.h new file mode 100644 index 00000000..cc85ed6e --- /dev/null +++ b/wspr/wspr.h @@ -0,0 +1,263 @@ +/* + This file is part of program wsprd, a detector/demodulator/decoder + for the Weak Signal Propagation Reporter (WSPR) mode. + + Copyright 2001-2015, Joe Taylor, K1JT + + Much of the present code is based on work by Steven Franke, K9AN, + which in turn was based on earlier work by K1JT. + + Copyright 2014-2015, Steven Franke, K9AN + + License: GNU GPL v3 + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + */ + +#pragma once + +#include "types.h" +#include "misc.h" +#include "../libcsdr.h" +#include "fano.h" +#include "jelinek.h" + +#include +#include +#include +#include +#include +#include + +#ifndef TRY_YIELD + #define YIELD_EVERY_N_TIMES 64 + #define TRY_YIELD wspr_try_yield("TRY_YIELD", 0) + #define TRY_YIELD_DECODE wspr_try_yield("TRY_YIELD_DECODE", 1) +#endif + +// FIXME +#define ext_send_msg(...) +#define ext_send_msg_encoded(...) +#define ext_send_msg_data(...) 0 + +//#define WSPR_DEBUG_MSG true +#define WSPR_DEBUG_MSG false + +//#define WSPR_PRINTF +#ifdef WSPR_PRINTF + #define wprintf(fmt, ...) \ + fprintf(stderr, fmt, ## __VA_ARGS__) + + #define wdprintf(fmt, ...) \ + fprintf(stderr, "WSPR %3ds ", timer_sec() - passes_start); \ + fprintf(stderr, fmt, ## __VA_ARGS__) +#else + #define wprintf(fmt, ...) + #define wdprintf(fmt, ...) +#endif + +#define WSPR_FLOAT +#ifdef WSPR_FLOAT + typedef float WSPR_REAL_t; + typedef float WSPR_CPX_t; + typedef struct { + float re; + float im; + } WSPR_COMPLEX_t; + + #define WSPR_FFTW_COMPLEX fftwf_complex + #define WSPR_FFTW_MALLOC fftwf_malloc + #define WSPR_FFTW_FREE fftwf_free + #define WSPR_FFTW_PLAN fftwf_plan + #define WSPR_FFTW_PLAN_DFT_1D fftwf_plan_dft_1d + #define WSPR_FFTW_DESTROY_PLAN fftwf_destroy_plan + #define WSPR_FFTW_EXECUTE fftwf_execute +#else + typedef double WSPR_REAL_t; + typedef double WSPR_CPX_t; + typedef struct { + double re; + double im; + } WSPR_COMPLEX_t; + + #define WSPR_FFTW_COMPLEX fftw_complex + #define WSPR_FFTW_MALLOC fftw_malloc + #define WSPR_FFTW_FREE fftw_free + #define WSPR_FFTW_PLAN fftw_plan + #define WSPR_FFTW_PLAN_DFT_1D fftw_plan_dft_1d + #define WSPR_FFTW_DESTROY_PLAN fftw_destroy_plan + #define WSPR_FFTW_EXECUTE fftw_execute +#endif + +#define K_PI (3.14159265358979323846) + +#define SYMTIME (FSPS / FSRATE) // symbol time: 256 samps @ 375 srate, 683 ms, 1.46 Hz + +#define SRATE 375 // design sample rate +#define FSRATE 375.0 +#define CTIME 120 // capture time secs +#define TPOINTS (SRATE * CTIME) + +#define FMIN -110 // frequency range to search +#define FMAX 110 +//#define FMIN -150 // frequency range to search +//#define FMAX 150 +#define BW_MAX 300.0 // +/- 150 Hz + +// samples per symbol (at SRATE) +#define FSPS 256.0 // round(SYMTIME * SRATE) +#define SPS 256 // (int) FSPS +#define NFFT (SPS*2) +#define HSPS (SPS/2) + +// groups +#define GROUPS (TPOINTS/NFFT) +#define FPG 4 // FFTs per group + +#define NSYM_162 162 +#define FNSYM_162 162.0 +#define HSYM_81 (NSYM_162/2) +#define FHSYM_81 (FNSYM_162/2) + +#define NBITS HSYM_81 + +#define LEN_DECODE ((NBITS + 7) / 8) +#define LEN_CALL (12 + SPACE_FOR_NULL) // max of AANLL, ppp/AANLLL, AANLL/ss, plus 2 +#define LEN_C_L_P (22 + SPACE_FOR_NULL) // 22 = 12 call, sp, 6 grid, sp, 2 pwr +#define LEN_GRID (6 + SPACE_FOR_NULL) + +#define WSPR_STACKSIZE 2000 + +#define NPK 256 +#define MAX_NPK 12 +#define MAX_NPK_OLD 8 + +typedef struct { + bool ignore; + float freq0, snr0, drift0, sync0; + int shift0, bin0; + int freq_idx, flags; + char snr_call[LEN_CALL]; +} pk_t; + +#define WSPR_F_BIN 0x0fff +#define WSPR_F_DECODING 0x1000 +#define WSPR_F_DELETE 0x2000 +#define WSPR_F_DECODED 0x4000 +#define WSPR_F_IMAGE 0x8000 + +// assigned constants +extern int nffts; +extern int nbins_411; +extern int hbins_205; + +typedef struct { + float freq; + char call[LEN_CALL]; + int hour, min; + float snr, dt_print, drift1; + double freq_print; + char c_l_p[LEN_C_L_P]; +} decode_t; + +typedef struct { + bool init; + int rx_chan; + int ping_pong, fft_ping_pong, decode_ping_pong; + int capture; + int status, status_resume; + bool send_error, abort_decode; + + // csdr + bool need_decode; + float *input_buffer; + int the_bufsize; + FILE *demod_pipe; + + // options + int quickmode, medium_effort, more_candidates, stackdecoder, subtraction; + #define WSPR_TYPE_2MIN 2 + #define WSPR_TYPE_15MIN 15 + int wspr_type; + + // sampler + bool reset, tsync; + int last_min, last_sec; + int decim, didx, group; + double fi; + + // FFT task + WSPR_FFTW_COMPLEX *fftin, *fftout; + WSPR_FFTW_PLAN fftplan; + int FFTtask_group; + + // computed by sampler or FFT task, processed by decode task + #define N_PING_PONG 2 + time_t utc[N_PING_PONG]; + WSPR_CPX_t i_data[N_PING_PONG][TPOINTS], q_data[N_PING_PONG][TPOINTS]; + float pwr_samp[N_PING_PONG][NFFT][FPG*GROUPS]; + float pwr_sampavg[N_PING_PONG][NFFT]; + + // decode task + float min_snr, snr_scaling_factor; + struct snode *stack; + float dialfreq_MHz, cf_offset; + u1_t symbols[NSYM_162], decdata[LEN_DECODE], channel_symbols[NSYM_162]; + char callsign[LEN_CALL], call_loc_pow[LEN_C_L_P], grid[LEN_GRID]; + decode_t deco[NPK]; +} wspr_t; + +// configuration +extern int bfo; + +void wspr_try_yield(const char *from, int check_for_data); +void wspr_decode_init(); +void wspr_init(char *demod_pipe_name, int decimate, int the_bufsize); +void wspr_data(WSPR_REAL_t *isamps, int nsamps); +void wspr_file(char *filename, int the_bufsize); +void wspr_decode_old(wspr_t *w); +void wspr_decode(wspr_t *w); +void wspr_send_peaks(wspr_t *w, pk_t *pk, int npk); +void wspr_hash_init(); + +void sync_and_demodulate( + WSPR_CPX_t *id, WSPR_CPX_t *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 renormalize(wspr_t *w, float psavg[], float smspec[]); + +void unpack50(u1_t *dat, u4_t *call_28b, u4_t *grid_pwr_22b, u4_t *grid_15b, u4_t *pwr_7b); +int unpackcall(u4_t call_28b, char *call); +int unpackgrid(u4_t grid_15b, char *grid); +int unpackpfx(int32_t nprefix, char *call); +void deinterleave(unsigned char *sym); +int unpk_(u1_t *decdata, char *call_loc_pow, char *callsign, char *grid, int *dBm); +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); + +int snr_comp(const void *elem1, const void *elem2); +int freq_comp(const void *elem1, const void *elem2); + +typedef struct { + double lat, lon; +} latLon_t; + +void set_reporter_grid(char *grid); +double grid_to_distance_km(char *grid); +int latLon_to_grid6(latLon_t *loc, char *grid); diff --git a/wspr/wspr_main.c b/wspr/wspr_main.c new file mode 100644 index 00000000..a72a357b --- /dev/null +++ b/wspr/wspr_main.c @@ -0,0 +1,385 @@ +/* + This file is part of program wsprd, a detector/demodulator/decoder + for the Weak Signal Propagation Reporter (WSPR) mode. + + Copyright 2001-2015, Joe Taylor, K1JT + + Much of the present code is based on work by Steven Franke, K9AN, + which in turn was based on earlier work by K1JT. + + Copyright 2014-2015, Steven Franke, K9AN + + License: GNU GPL v3 + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + */ + +#include "wspr.h" + +#include +#include +#include +#include +#include + +// wspr_status +#define NONE 0 +#define IDLE 1 +#define SYNC 2 +#define RUNNING 3 +#define DECODING 4 + +static wspr_t wspr; + +// assigned constants +int nffts, nbins_411, hbins_205; + +// computed constants +static float window[NFFT]; + +// loaded from admin configuration +int bfo; + +void wspr_try_yield(const char *from, int check_for_data) +{ + wspr_t *w = &wspr; + + if (check_for_data) { + int n; + do { + // stdin was set non-blocking by csdr calling code + n = fread(w->input_buffer, sizeof(float), w->the_bufsize, stdin); + if (n > 0) { + wspr_data(w->input_buffer, n/2); // n/2 for number of complex samples + } else { + if (w->status != DECODING) { + // Not decoding: Sleep a little while so we don't hog the CPU. + // Otherwise the decoder needs to run at full speed and use regular + // Linux scheduling/time-slicing to share the CPU. + // This could be improved by not using read polling (FIXME). + usleep(100); + } + } + } while (n > 0); + } + + fflush(stdout); + sched_yield(); +} + +const char *status_str[] = { "none", "idle", "sync", "running", "decoding" }; + +static void wspr_status(wspr_t *w, int status, int resume) +{ + wprintf("WSPR wspr_status: current %d-%s, set to %d-%s\n", + w->status, status_str[w->status], status, status_str[status]); + + // sends status as a control character (0x01 - 0x04) in the HTML stream + fprintf(w->demod_pipe, "%c", status); fflush(w->demod_pipe); + ext_send_msg(w->rx_chan, WSPR_DEBUG_MSG, "EXT WSPR_STATUS=%d", status); + + w->status = status; + if (resume != NONE) { + wprintf("WSPR wspr_status: will resume to %d-%s\n", resume, status_str[resume]); + w->status_resume = resume; + } +} + +// compute FFTs incrementally during data capture rather than all at once at the beginning of the decode +void WSPR_FFT() +{ + int i,j,k; + + wspr_t *w = &wspr; + int grp = w->FFTtask_group; + int first = grp*FPG, last = first+FPG; + //wprintf("WSPR FFT pp=%d grp=%d (%d-%d)\n", w->fft_ping_pong, grp, first, last); + + // Do ffts over 2 symbols, stepped by half symbols + WSPR_CPX_t *id = w->i_data[w->fft_ping_pong], *qd = w->q_data[w->fft_ping_pong]; + + float maxiq = 1e-66, maxpwr = 1e-66; + int maxi=0; + float savg[NFFT]; + memset(savg, 0, sizeof(savg)); + + // FIXME: A large noise burst can washout the w->pwr_sampavg which prevents a proper + // peak list being created in the decoder. An individual signal can be decoded fine + // in the presence of the burst. But if there is no peak in the list the decoding + // process is never started! We've seen this problem fairly frequently. + + for (i=first; ifftin[j][0] = id[k] * window[j]; + w->fftin[j][1] = qd[k] * window[j]; + //if (i==0) wprintf("IN %d %fi %fq\n", j, w->fftin[j][0], w->fftin[j][1]); + } + + //u4_t start = timer_us(); + TRY_YIELD; + WSPR_FFTW_EXECUTE(w->fftplan); + TRY_YIELD; + //if (i==0) wprintf("512 FFT %.1f us\n", (float)(timer_us()-start)); + + // NFFT = SPS*2 + // unwrap fftout: + // |1---> SPS + // 2--->| SPS + + for (j=0; j (NFFT-1)) + k -= NFFT; + float ii = w->fftout[k][0]; + float qq = w->fftout[k][1]; + float pwr = ii*ii + qq*qq; + //if (i==0) wprintf("OUT %d %fi %fq\n", j, ii, qq); + if (ii > maxiq) { maxiq = ii; } + if (pwr > maxpwr) { maxpwr = pwr; maxi = k; } + w->pwr_samp[w->fft_ping_pong][j][i] = pwr; + w->pwr_sampavg[w->fft_ping_pong][j] += pwr; + savg[j] += pwr; + } + TRY_YIELD; + } + + // send spectrum data to client + float smspec[nbins_411]; + renormalize(w, savg, smspec); + + for (j=0; jsnr_scaling_factor; + } + + #define DATW nbins_411 + #define OUTW 1024 + #define LEVEL_COMP -30 + float out[OUTW]; + + for (i = 0; i < OUTW; i++) { + out[i] = smspec[(int) round((float) i / OUTW * (DATW-1))] + LEVEL_COMP; + } + + fwrite(&out, sizeof(float), OUTW, stdout); + + //wprintf("WSPR_FFTtask group %d:%d %d-%d(%d) %f %f(%d)\n", w->fft_ping_pong, grp, first, last, nffts, maxiq, maxpwr, maxi); + + time_t t; time(&t); struct tm tm; gmtime_r(&t, &tm); + wprintf(" %03d-%03d(%d) pp%d %02d:%02d\r", first, last, nffts, w->fft_ping_pong, tm.tm_min, tm.tm_sec); fflush(stderr); +} + +void wspr_send_peaks(wspr_t *w, pk_t *pk, int npk) +{ + char peaks_s[NPK*(6+1 + LEN_CALL) + 16]; + char *s = peaks_s; + *s = '\0'; + int j, n; + for (j=0; j < npk; j++) { + int bin_flags = pk[j].bin0 | pk[j].flags; + n = sprintf(s, "%d:%s:", bin_flags, pk[j].snr_call); s += n; + } + ext_send_msg_encoded(w->rx_chan, WSPR_DEBUG_MSG, "EXT", "WSPR_PEAKS", "%s", peaks_s); +} + +void WSPR_Decoder() +{ + u4_t start=0; + + wspr_t *w = &wspr; + wspr_status(w, SYNC, SYNC); + + while (1) { + if (w->need_decode) { + w->need_decode = false; + w->abort_decode = false; + + wspr_status(w, DECODING, NONE); + start = timer_ms(); + wspr_decode(w); + wprintf("WSPR Decoder TOTAL %.1f sec\n", (float)(timer_ms()-start)/1000.0); + + if (w->abort_decode) + wprintf("WSPR decoder aborted\n"); + + wspr_status(w, w->status_resume, NONE); + } + + // important that this is the decode version that checks for input data + TRY_YIELD_DECODE; + } +} + +static int int_decimate; + +void wspr_data(WSPR_REAL_t *isamps, int nsamps) +{ + wspr_t *w = &wspr; + int i; + WSPR_COMPLEX_t *samps = (WSPR_COMPLEX_t *) isamps; + + //wprintf("WD%d didx %d send_error %d reset %d\n", w->capture, w->didx, w->send_error, w->reset); + if (w->send_error) { + wprintf("WSPR STOP send_error %d\n", w->send_error); + w->send_error = FALSE; + w->capture = FALSE; + w->reset = TRUE; + } + + if (w->reset) { + w->ping_pong = w->decim = w->didx = w->group = 0; + w->fi = 0; + w->tsync = FALSE; + w->status_resume = IDLE; // decoder finishes after we stop capturing + w->reset = FALSE; + } + + if (!w->capture) { + return; + } + + // because of buffering this routine can be called with gaps > 1 sec (i.e. data gets bunched up) + time_t t; time(&t); struct tm tm; gmtime_r(&t, &tm); + + if (tm.tm_sec != w->last_sec) { + //fprintf(stderr, "WSPR %02d:%02d sync=%d\n", tm.tm_min, tm.tm_sec, w->tsync); fflush(stderr); + if (tm.tm_min&1 && (tm.tm_sec >= 40 && tm.tm_sec <= 43)) // leave enough time for upload + w->abort_decode = true; + } + + if (w->tsync == FALSE) { // sync to even minute boundary + #define START_DELAY_SECS 1 + #define START_WINDOW 2 + if (!(tm.tm_min&1) && tm.tm_sec >= START_DELAY_SECS && tm.tm_sec <= (START_DELAY_SECS + START_WINDOW)) { + + if (w->status == RUNNING && !w->need_decode) { + w->decode_ping_pong = w->ping_pong; + w->need_decode = true; + } + + w->ping_pong ^= 1; + wprintf("\nWSPR SYNC ping_pong %d, %s", w->ping_pong, ctime(&t)); + w->decim = w->didx = w->group = 0; + w->fi = 0; + if (w->status != DECODING) + wspr_status(w, RUNNING, RUNNING); + w->tsync = TRUE; + } + } + + if (tm.tm_sec != w->last_sec) { + w->last_min = tm.tm_min; + w->last_sec = tm.tm_sec; + } + + if (w->didx == 0) { + memset(&w->pwr_sampavg[w->ping_pong][0], 0, sizeof(w->pwr_sampavg[0])); + } + + if (w->group == 0) w->utc[w->ping_pong] = t; + + WSPR_CPX_t *idat = w->i_data[w->ping_pong], *qdat = w->q_data[w->ping_pong]; + + for (i=0; idecim++ < (int_decimate-1)) + continue; + w->decim = 0; + + if (w->didx >= TPOINTS) + return; + + if (w->group == 3) w->tsync = FALSE; // arm re-sync + + WSPR_CPX_t re = (WSPR_CPX_t) samps[i].re; + WSPR_CPX_t im = (WSPR_CPX_t) samps[i].im; + idat[w->didx] = re; + qdat[w->didx] = im; + + if ((w->didx % NFFT) == (NFFT-1)) { + w->fft_ping_pong = w->ping_pong; + w->FFTtask_group = w->group-1; + if (w->group) WSPR_FFT(); // skip first to pipeline + w->group++; + } + w->didx++; + } +} + +void wspr_init(char *demod_pipe_name, int decimate, int the_bufsize) +{ + int i; + wspr_t *w; + + assert(FSPS == round(SYMTIME * FSRATE)); + assert(SPS == (int) FSPS); + assert(HSPS == (SPS/2)); + + nffts = FPG * floor(GROUPS-1) -1; + nbins_411 = ceilf(NFFT * BW_MAX / FSRATE) +1; + hbins_205 = (nbins_411-1)/2; + + wspr_decode_init(); + + w = &wspr; + memset(w, 0, sizeof(wspr_t)); + + w->medium_effort = 1; + w->wspr_type = WSPR_TYPE_2MIN; + + w->fftin = (WSPR_FFTW_COMPLEX*) WSPR_FFTW_MALLOC(sizeof(WSPR_FFTW_COMPLEX)*NFFT); + w->fftout = (WSPR_FFTW_COMPLEX*) WSPR_FFTW_MALLOC(sizeof(WSPR_FFTW_COMPLEX)*NFFT); + w->fftplan = WSPR_FFTW_PLAN_DFT_1D(NFFT, w->fftin, w->fftout, FFTW_FORWARD, FFTW_ESTIMATE); + + w->status_resume = IDLE; + w->tsync = FALSE; + w->capture = 0; + time_t t; time(&t); struct tm tm; gmtime_r(&t, &tm); + w->last_min = tm.tm_min; + w->last_sec = tm.tm_sec; + w->abort_decode = false; + w->send_error = false; + + for (i=0; i < NFFT; i++) { + window[i] = sin(i * K_PI/(NFFT-1)); + } + + // sdr source sampling rate determines if we rationally decimate here or rely on csdr fractional decimation + int_decimate = decimate; + wprintf("WSPR %s int_decimate=%d sps=%d NFFT=%d nbins_411=%d\n", + (int_decimate != 1)? "LOCAL INTEGER DECIMATION" : "CSDR FRACTIONAL DECIMATION", + int_decimate, SPS, NFFT, nbins_411); + + w->input_buffer = (float *) malloc(the_bufsize * sizeof(float)); + w->the_bufsize = the_bufsize; + + w->demod_pipe = fopen(demod_pipe_name, "w"); + if (w->demod_pipe == NULL) { + wprintf("WSPR fopen FAIL %s\n", demod_pipe_name); + exit(-1); + } + + //fprintf(w->demod_pipe, "WSPR decoder running..\n"); fflush(w->demod_pipe); + + bfo = 750; + w->dialfreq_MHz = 7040.1/1e3 - bfo/1e6; + w->cf_offset = 100; + wprintf("WSPR dialfreq=%f bfo=%d cf_offset=%f\n", w->dialfreq_MHz, bfo, w->cf_offset); + w->rx_chan = 0; + w->reset = true; + w->capture = true; + + WSPR_Decoder(); +} diff --git a/wspr/wspr_util.c b/wspr/wspr_util.c new file mode 100644 index 00000000..49f246cc --- /dev/null +++ b/wspr/wspr_util.c @@ -0,0 +1,525 @@ +/* + This file is part of program wsprd, a detector/demodulator/decoder + for the Weak Signal Propagation Reporter (WSPR) mode. + + File name: wsprd_util.c + + Copyright 2001-2015, Joe Taylor, K1JT + + Most of the code is based on work by Steven Franke, K9AN, which + in turn was based on earlier work by K1JT. + + Copyright 2014-2015, Steven Franke, K9AN + + License: GNU GPL v3 + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + */ + +#include "wspr.h" +#include "nhash.h" + +#include +#include +#include +#include +#include +#include +#include + +// call_28b bits 27-0 = first 28 bits of *d (big-endian) +// ......d0 ......d1 ......d2 ......d3 +// 22222222 11111111 11000000 0000.... +// 76543210 98765432 10987654 3210.... + +// grid_pwr_22b bits 21-0 = 22 subsequent bits +// ......d3 ......d4 ......d5 ......d6 +// ....2211 11111111 00000000 00...... +// ....1098 76543210 98765432 10...... +// gggg gggggggg gggppppp pp g: 15-bit grid, p: 7-bit pwr + +void unpack50(u1_t *d, u4_t *call_28b, u4_t *grid_pwr_22b, u4_t *grid_15b, u4_t *pwr_7b) +{ + *call_28b = (d[0]<<20) | (d[1]<<12) | (d[2]<<4) | (d[3]>>4); + *grid_pwr_22b = ((d[3]&0xf)<<18) | (d[4]<<10) | (d[5]<<2) | (d[6]>>6); + *grid_15b = *grid_pwr_22b >> 7; + *pwr_7b = *grid_pwr_22b & 0x7f; +} + +static const char *c = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ "; // 37 characters + +// valid call: AANLLL, A = alpha-numeric, N = numeric, L = letter or space +int unpackcall(u4_t call_28b, char *call) +{ + int i; + char tmp[7]; + + memset(call, 0, LEN_CALL); + + if (call_28b < 262177560) { // = 37*36*10*27*27*27 = 0xfa08318 (28-bits) + tmp[5] = c[call_28b%27+10]; // letter, space + call_28b /= 27; + tmp[4] = c[call_28b%27+10]; // letter, space + call_28b /= 27; + tmp[3] = c[call_28b%27+10]; // letter, space + call_28b /= 27; + tmp[2] = c[call_28b%10]; // number + call_28b /= 10; + tmp[1] = c[call_28b%36]; // letter, number + call_28b /= 36; + tmp[0] = c[call_28b]; // letter, number, space + tmp[6] = '\0'; + + // remove leading whitespace + for (i=0; i<5; i++) { + if (tmp[i] != ' ') + break; + } + sprintf(call, "%-6s", &tmp[i]); + + // remove trailing whitespace + for (i=0; i<6; i++) { + if (call[i] == ' ') { + call[i] = '\0'; + } + } + } else { + return 0; + } + return 1; +} + +int unpackgrid(u4_t grid_15b, char *grid) +{ + int dlat, dlong; + + memset(grid, 0, LEN_GRID); + + if (grid_15b < 32400) { // = 0x7e90 (15-bits) + dlat = (grid_15b%180)-90; + dlong = (grid_15b/180)*2 - 180 + 2; + if (dlong < -180) + dlong = dlong+360; + if (dlong > 180) + dlong = dlong+360; + int nlong = 60.0*(180.0-dlong)/5.0; + int n1 = nlong/240; + int n2 = (nlong - 240*n1)/24; + //int n3 = nlong -40*n1 - 24*n2; + grid[0] = c[10+n1]; + grid[2] = c[n2]; + + int nlat = 60.0*(dlat+90)/2.5; + n1 = nlat/240; + n2 = (nlat-240*n1)/24; + //n3 = nlong - 240*n1 - 24*n2; + grid[1] = c[10+n1]; + grid[3] = c[n2]; + } else { + strcpy(grid,"XXXX"); + return 0; + } + return 1; +} + +int unpackpfx(int32_t nprefix, char *call) +{ + char nc, pfx[4]={'\0'}, tmpcall[7]; + int i; + int32_t n; + + strcpy(tmpcall,call); + if( nprefix < 60000 ) { // < 0xea60 + // add a prefix of 1 to 3 characters + n=nprefix; + for (i=2; i>=0; i--) { + nc = n%37; // letter, number, space + if( (nc >= 0) & (nc <= 9) ) { + pfx[i] = nc + '0'; + } + else if( (nc >= 10) & (nc <= 35) ) { + pfx[i] = nc-10 + 'A'; + } + else { + pfx[i] = ' '; + } + n = n/37; + } + + char * p = strrchr(pfx,' '); + strcpy(call, p ? p + 1 : pfx); + strncat(call,"/",1); + strncat(call,tmpcall,strlen(tmpcall)); + + } else { + // add a suffix of 1 or 2 characters + nc=nprefix-60000; + if( (nc >= 0) & (nc <= 9) ) { + pfx[0] = nc + '0'; + strcpy(call,tmpcall); + strncat(call,"/",1); + strncat(call,pfx,1); + } + else if( (nc >= 10) & (nc <= 35) ) { + pfx[0] = nc-10 + 'A'; + strcpy(call,tmpcall); + strncat(call,"/",1); + strncat(call,pfx,1); + } + else if( (nc >= 36) & (nc <= 125) ) { + pfx[0]=(nc-26)/10 + '0'; + pfx[1]=(nc-26)%10 + '0'; + strcpy(call,tmpcall); + strncat(call,"/",1); + strncat(call,pfx,2); + } + else { + return 0; + } + } + return 1; +} + +void deinterleave(unsigned char *sym) +{ + unsigned char tmp[NSYM_162]; + unsigned char p, i, j; + + p=0; + i=0; + while (p> 32; + if (j < NSYM_162 ) { + tmp[p]=sym[j]; + p=p+1; + } + i=i+1; + } + for (i=0; i %d\n", htsize, htsize*2); + htsize *= 2; + ht = (hashtab_t *) realloc(ht, sizeof(hashtab_t) * htsize); + memset(ht + htsize/2, 0, sizeof(hashtab_t) * htsize/2); + //wprintf("W-HASH %d 0x%04x exp new %s\n", htsize/2, hash, call); + ht[htsize/2].hash = hash; + strcpy(ht[htsize/2].call, call); + } +} + +static char *hash_lookup(int hash) +{ + int i; + + for (i=0; i < htsize; i++) { + if (ht[i].call[0] == '\0') + break; + if (ht[i].hash == hash) { + //wprintf("W-HASH %d 0x%04x lookup %s\n", i, hash, ht[i].call); + return ht[i].call; + } + } + + //wprintf("W-HASH 0x%04x lookup FAIL\n", hash); + return NULL; +} + +int unpk_(u1_t *decdata, char *call_loc_pow, char *callsign, char *grid, int *dBm) +{ + int rtn = 0; + u4_t call_28b, grid_pwr_22b, grid_15b, pwr_7b; + int n3, ndbm, nadd; + + memset(call_loc_pow, 0, LEN_C_L_P); + + unpack50(decdata, &call_28b, &grid_pwr_22b, &grid_15b, &pwr_7b); + if (!unpackcall(call_28b, callsign)) return -1; + if (!unpackgrid(grid_15b, grid)) return -2; + + // ntype -64..63, but this is NOT twos complement (just biasing) + int ntype = pwr_7b - 64; + + /* + Based on the value of ntype, decide whether this is a Type 1, 2, or 3 message. + + * Type 1: 6 digit call, grid, power - ntype is positive and is a member + of the set {0,3,7,10,13,17,20...60} + + * Type 2: extended callsign, power - ntype is positive but not + a member of the set of allowed powers + + * Type 3: hash, 6 digit grid, power - ntype is negative. + */ + + if ((ntype >= 0) && (ntype <= 62)) { + int nu = ntype%10; + if (nu == 0 || nu == 3 || nu == 7) { + *dBm = ndbm = ntype; + sprintf(call_loc_pow, "%s %s %2d", callsign, grid, ndbm); + hash_update(callsign); + rtn = 1; + } else { + nadd = nu; + if( nu > 3 ) nadd=nu-3; + if( nu > 7 ) nadd=nu-7; + // Nggggggg gggggggg + n3 = grid_15b + 32768*(nadd-1); // 32768 = 0x8000, (nadd-1) = 0..2 + if (!unpackpfx(n3, callsign)) return -3; + + grid[0] = '\0'; // no grid info for TYPE2 + + *dBm = ndbm = ntype-nadd; + sprintf(call_loc_pow, "%s no_grid %2d", callsign, ndbm); + + int nu = ndbm%10; + if (nu == 0 || nu == 3 || nu == 7 || nu == 10) { // make sure power is OK + hash_update(callsign); + } else return -4; + rtn = 2; + } + } else + + if (ntype < 0) { + *dBm = ndbm = -(ntype+1); + + memset(grid, 0, LEN_GRID); + // this because the grid6 L1L2N1N2L3L4 was packed as L2N1N2L3L4L1 to be + // compatible with the AANLLL required by pack_call() / unpackcall() + sprintf(grid, "%c%.5s", callsign[5], callsign); + + int nu = ndbm%10; + if ((nu != 0 && nu != 3 && nu != 7 && nu != 10) || + !isalpha(grid[0]) || !isalpha(grid[1]) || + !isdigit(grid[2]) || !isdigit(grid[3])) { + // not testing 4'th and 5'th chars because of this case: JO33 40 + // grid is only 4 chars even though this is a hashed callsign... + // isalpha(grid[4]) && isalpha(grid[5]) ) ) { + return -5; + } + + // get hashed callsign + // wspr_t:callsign shouldn't contain HTML-confusing angle brackets + char *callp = hash_lookup((grid_pwr_22b-ntype-64) / 128); + sprintf(callsign, "%s", callp? callp : "..."); + sprintf(call_loc_pow, "<%s> %s %2d", callsign, grid, ndbm); + + // I don't know what to do with these... They show up as "A000AA" grids. + if (ntype == -64) return -6; + rtn = 3; + } else { + return -7; + } + return rtn; +} + +int snr_comp(const void *elem1, const void *elem2) +{ + const pk_t *e1 = (const pk_t*) elem1, *e2 = (const pk_t*) elem2; + int r = (e1->snr0 < e2->snr0)? 1 : ((e1->snr0 > e2->snr0)? -1:0); // NB: comparison reversed to sort in descending order + return r; +} + +int freq_comp(const void *elem1, const void *elem2) +{ + const pk_t *e1 = (const pk_t*) elem1, *e2 = (const pk_t*) elem2; + int r = (e1->freq0 < e2->freq0)? -1 : ((e1->freq0 > e2->freq0)? 1:0); + return r; +} + +#define DEG_2_RAD(deg) ((deg) * K_PI / 180.0) + +#define latLon_deg_to_rad(loc) \ + loc.lat = DEG_2_RAD(loc.lat); \ + loc.lon = DEG_2_RAD(loc.lon); + +// field square subsquare (extended square) +// A-R 0-9 a-x 0-9 +// #18 #10 #24 #10 + +#define SQ_LON_DEG 2.0 +#define SQ_LAT_DEG 1.0 +#define SUBSQ_PER_SQ 24.0 +#define SUBSQ_LON_DEG (SQ_LON_DEG / SUBSQ_PER_SQ) +#define SUBSQ_LAT_DEG (SQ_LAT_DEG / SUBSQ_PER_SQ) + +#define SQ_PER_FLD 10.0 +#define FLD_DEG_LON (SQ_PER_FLD * SQ_LON_DEG) +#define FLD_DEG_LAT (SQ_PER_FLD * SQ_LAT_DEG) + +static void grid_to_latLon(char *grid, latLon_t *loc) +{ + double lat, lon; + char c; + int slen = strlen(grid); + + loc->lat = loc->lon = 999.0; + if (slen < 4) return; + + c = tolower(grid[0]); + if (c < 'a' || c > 'r') return; + lon = (c-'a')*20 - 180; + + c = tolower(grid[1]); + if (c < 'a' || c > 'r') return; + lat = (c-'a')*10 - 90; + + c = grid[2]; + if (c < '0' || c > '9') return; + lon += (c-'0') * SQ_LON_DEG; + + c = grid[3]; + if (c < '0' || c > '9') return; + lat += (c-'0') * SQ_LAT_DEG; + + if (slen != 6) { // assume center of square (i.e. "....ll") + lon += SQ_LON_DEG /2.0; + lat += SQ_LAT_DEG /2.0; + } else { + c = tolower(grid[4]); + if (c < 'a' || c > 'x') return; + lon += (c-'a') * SUBSQ_LON_DEG; + + c = tolower(grid[5]); + if (c < 'a' || c > 'x') return; + lat += (c-'a') * SUBSQ_LAT_DEG; + + lon += SUBSQ_LON_DEG /2.0; // assume center of sub-square (i.e. "......44") + lat += SUBSQ_LAT_DEG /2.0; + } + + loc->lat = lat; + loc->lon = lon; + //wprintf("GRID %s%s = (%f, %f)\n", grid, (slen != 6)? "[ll]":"", lat, lon); +} + +static const char *field = "ABCDEFGHIJKLMNOPQR"; +static const char *square = "0123456789"; +static const char *subsquare = "abcdefghijklmnopqrstuvwx"; + +int latLon_to_grid6(latLon_t *loc, char *grid6) +{ + int i; + double r, lat, lon; + + // longitude + lon = loc->lon + 180.0; + if (lon < 0 || lon >= 360.0) return -1; + i = (int) lon / FLD_DEG_LON; + grid6[0] = field[i]; + r = lon - (i * FLD_DEG_LON); + + i = (int) floor(r / SQ_LON_DEG); + grid6[2] = square[i]; + r = r - (i * SQ_LON_DEG); + + i = (int) floor(r * (SUBSQ_PER_SQ / SQ_LON_DEG)); + grid6[4] = subsquare[i]; + + // latitude + lat = loc->lat + 90.0; + if (lat < 0 || lat >= 180.0) return -1; + i = (int) lat / FLD_DEG_LAT; + grid6[1] = field[i]; + r = lat - (i * FLD_DEG_LAT); + + i = (int) floor(r / SQ_LAT_DEG); + grid6[3] = square[i]; + r = r - (i * SQ_LAT_DEG); + + i = (int) floor(r * (SUBSQ_PER_SQ / SQ_LAT_DEG)); + grid6[5] = subsquare[i]; + + return 0; +} + +static latLon_t r_loc; + +void set_reporter_grid(char *grid) +{ + grid_to_latLon(grid, &r_loc); + if (r_loc.lat != 999.0) + latLon_deg_to_rad(r_loc); +} + +double grid_to_distance_km(char *grid) +{ + if (r_loc.lat == 999.0 || *grid == '\0') + return 0; + + latLon_t loc; + grid_to_latLon(grid, &loc); + latLon_deg_to_rad(loc); + + double delta_lat = loc.lat - r_loc.lat; + delta_lat /= 2.0; + delta_lat = sin(delta_lat); + delta_lat *= delta_lat; + double delta_lon = loc.lon - r_loc.lon; + delta_lon /= 2.0; + delta_lon = sin(delta_lon); + delta_lon *= delta_lon; + + double t = delta_lat + (delta_lon * cos(loc.lat) * cos(r_loc.lat)); + #define EARTH_RADIUS_KM 6371.0 + double km = EARTH_RADIUS_KM * 2.0 * atan2(sqrt(t), sqrt(1.0-t)); + return km; +}