diff --git a/CMakeLists.txt b/CMakeLists.txt index 2aca941..8abd6e9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -545,26 +545,6 @@ set (wsjt_qt_UISRCS wf_palette_design_dialog.ui ) -set (wsprsim_CSRCS - lib/wsprd/wsprsim.c - lib/wsprd/wsprsim_utils.c - lib/wsprd/wsprd_utils.c - lib/wsprd/fano.c - lib/wsprd/tab.c - lib/wsprd/nhash.c - ) - -set (wsprd_CSRCS - lib/wsprd/wsprd.c - lib/wsprd/wsprsim_utils.c - lib/wsprd/wsprd_utils.c - lib/wsprd/fano.c - lib/wsprd/jelinek.c - lib/wsprd/tab.c - lib/wsprd/nhash.c - lib/init_random_seed.c - ) - set (wsjtx_UISRCS mainwindow.ui about.ui @@ -594,8 +574,6 @@ set (all_CXXSRCS set (all_C_and_CXXSRCS ${sqlite3_CSRCS} ${wsjt_CSRCS} - ${wsprsim_CSRCS} - ${wsprd_CSRCS} ${all_CXXSRCS} ) diff --git a/lib/wsprcode/go.sh b/lib/wsprcode/go.sh deleted file mode 100644 index dfce99e..0000000 --- a/lib/wsprcode/go.sh +++ /dev/null @@ -1 +0,0 @@ - gfortran -o wsprcode -fbounds-check wsprcode.f90 nhash.c diff --git a/lib/wsprcode/nhash.c b/lib/wsprcode/nhash.c deleted file mode 100644 index 675022b..0000000 --- a/lib/wsprcode/nhash.c +++ /dev/null @@ -1,384 +0,0 @@ -/* - *------------------------------------------------------------------------------- - * - * This file is part of the WSPR application, Weak Signal Propagation Reporter - * - * File Name: nhash.c - * Description: Functions to produce 32-bit hashes for hash table lookup - * - * Copyright (C) 2008-2014 Joseph Taylor, K1JT - * License: GPL-3 - * - * 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 */ -#ifdef Win32 -#include "win_stdint.h" /* defines uint32_t etc */ -#else -#include /* defines uint32_t etc */ -#endif -//#include /* attempt to define endianness */ -//#ifdef linux -//# include /* attempt to define endianness */ -//#endif - -#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); - return c; -} - -//uint32_t __stdcall NHASH(const void *key, size_t length, uint32_t initval) diff --git a/lib/wsprcode/wspr_old_subs.f90 b/lib/wsprcode/wspr_old_subs.f90 deleted file mode 100644 index b1bbb7e..0000000 --- a/lib/wsprcode/wspr_old_subs.f90 +++ /dev/null @@ -1,937 +0,0 @@ -!------------------------------------------------------------------------------- -! -! This file is part of the WSPR application, Weak Signal Propagation Reporter -! -! File Name: wspr_old_subs.f90 -! Description: Utility subroutines from WSPR 2.0 -! -! Copyright (C) 2001-2014 Joseph Taylor, K1JT -! License: GPL-3 -! -! 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. -! -!------------------------------------------------------------------------------- - -subroutine deg2grid(dlong0,dlat,grid) - - real dlong !West longitude (deg) - real dlat !Latitude (deg) - character grid*6 - - dlong=dlong0 - if(dlong.lt.-180.0) dlong=dlong+360.0 - if(dlong.gt.180.0) dlong=dlong-360.0 - -! Convert to units of 5 min of longitude, working east from 180 deg. - nlong=60.0*(180.0-dlong)/5.0 - n1=nlong/240 !20-degree field - n2=(nlong-240*n1)/24 !2 degree square - n3=nlong-240*n1-24*n2 !5 minute subsquare - grid(1:1)=char(ichar('A')+n1) - grid(3:3)=char(ichar('0')+n2) - grid(5:5)=char(ichar('a')+n3) - -! Convert to units of 2.5 min of latitude, working north from -90 deg. - nlat=60.0*(dlat+90)/2.5 - n1=nlat/240 !10-degree field - n2=(nlat-240*n1)/24 !1 degree square - n3=nlat-240*n1-24*n2 !2.5 minuts subsquare - grid(2:2)=char(ichar('A')+n1) - grid(4:4)=char(ichar('0')+n2) - grid(6:6)=char(ichar('a')+n3) - - return -end subroutine deg2grid - -subroutine encode232(dat,nbytes,symbol,maxsym) - -! Convolutional encoder for a K=32, r=1/2 code. - - integer*1 dat(nbytes) !User data, packed 8 bits per byte - integer*1 symbol(maxsym) !Channel symbols, one bit per byte - integer*1 i1 - -! Layland-Lushbaugh polynomials for a K=32, r=1/2 convolutional code, -! and 8-bit parity lookup table. - - data npoly1/-221228207/,npoly2/-463389625/ - integer*1 partab(0:255) - data 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/ - - nstate=0 - k=0 - do j=1,nbytes - do i=7,0,-1 - i1=dat(j) - i4=i1 - if (i4.lt.0) i4=i4+256 - nstate=ior(ishft(nstate,1),iand(ishft(i4,-i),1)) - n=iand(nstate,npoly1) - n=ieor(n,ishft(n,-16)) - k=k+1 - symbol(k)=partab(iand(ieor(n,ishft(n,-8)),255)) - n=iand(nstate,npoly2) - n=ieor(n,ishft(n,-16)) - k=k+1 - symbol(k)=partab(iand(ieor(n,ishft(n,-8)),255)) - enddo - enddo - - return -end subroutine encode232 - -subroutine fano232(symbol,nbits,mettab,ndelta,maxcycles,dat,ncycles,metric,ierr) - -! Sequential decoder for K=32, r=1/2 convolutional code using -! the Fano algorithm. Translated from C routine for same purpose -! written by Phil Karn, KA9Q. - - parameter (MAXBITS=103) - parameter (MAXDAT=(MAXBITS+7)/8) - integer*1 symbol(0:2*MAXBITS-1) - integer*1 dat(MAXDAT) !Decoded user data, 8 bits per byte - integer mettab(0:255,0:1) !Metric table - -! These were the "node" structure in Karn's C code: - integer nstate(0:MAXBITS-1) !Encoder state of next node - integer gamma(0:MAXBITS-1) !Cumulative metric to this node - integer metrics(0:3,0:MAXBITS-1) !Metrics indexed by all possible Tx syms - integer tm(0:1,0:MAXBITS-1) !Sorted metrics for current hypotheses - integer ii(0:MAXBITS-1) !Current branch being tested - - logical noback - -! Layland-Lushbaugh polynomials for a K=32, r=1/2 convolutional code, -! and 8-bit parity lookup table. - - data npoly1/-221228207/,npoly2/-463389625/ - integer*1 partab(0:255) - data 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/ - - ntail=nbits-31 - -! Compute all possible branch metrics for each symbol pair. -! This is the only place we actually look at the raw input symbols - i4a=0 - i4b=0 - do np=0,nbits-1 - j=2*np - i4a=symbol(j) - i4b=symbol(j+1) - if (i4a.lt.0) i4a=i4a+256 - if (i4b.lt.0) i4b=i4b+256 - metrics(0,np) = mettab(i4a,0) + mettab(i4b,0) - metrics(1,np) = mettab(i4a,0) + mettab(i4b,1) - metrics(2,np) = mettab(i4a,1) + mettab(i4b,0) - metrics(3,np) = mettab(i4a,1) + mettab(i4b,1) - enddo - - np=0 - nstate(np)=0 - -! Compute and sort branch metrics from the root node - n=iand(nstate(np),npoly1) - n=ieor(n,ishft(n,-16)) - lsym=partab(iand(ieor(n,ishft(n,-8)),255)) - n=iand(nstate(np),npoly2) - n=ieor(n,ishft(n,-16)) - lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255)) - m0=metrics(lsym,np) - m1=metrics(ieor(3,lsym),np) - if(m0.gt.m1) then - tm(0,np)=m0 !0-branch has better metric - tm(1,np)=m1 - else - tm(0,np)=m1 !1-branch is better - tm(1,np)=m0 - nstate(np)=nstate(np) + 1 !Set low bit - endif - -! Start with best branch - ii(np)=0 - gamma(np)=0 - nt=0 - -! Start the Fano decoder - do i=1,nbits*maxcycles -! Look forward - ngamma=gamma(np) + tm(ii(np),np) - if(ngamma.ge.nt) then - -! Node is acceptable. If first time visiting this node, tighten threshold: - if(gamma(np).lt.(nt+ndelta)) nt=nt + & - ndelta * ((ngamma-nt)/ndelta) - -! Move forward - gamma(np+1)=ngamma - nstate(np+1)=ishft(nstate(np),1) - np=np+1 - if(np.eq.nbits-1) go to 100 !We're done! - - n=iand(nstate(np),npoly1) - n=ieor(n,ishft(n,-16)) - lsym=partab(iand(ieor(n,ishft(n,-8)),255)) - n=iand(nstate(np),npoly2) - n=ieor(n,ishft(n,-16)) - lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255)) - - if(np.ge.ntail) then - tm(0,np)=metrics(lsym,np) !We're in the tail, all zeros - else - m0=metrics(lsym,np) - m1=metrics(ieor(3,lsym),np) - if(m0.gt.m1) then - tm(0,np)=m0 !0-branch has better metric - tm(1,np)=m1 - else - tm(0,np)=m1 !1-branch is better - tm(1,np)=m0 - nstate(np)=nstate(np) + 1 !Set low bit - endif - endif - - ii(np)=0 !Start with best branch - go to 99 - endif - -! Threshold violated, can't go forward -10 noback=.false. - if(np.eq.0) noback=.true. - if(np.gt.0) then - if(gamma(np-1).lt.nt) noback=.true. - endif - - if(noback) then -! Can't back up, either. Relax threshold and look forward again -! to a better branch. - nt=nt-ndelta - if(ii(np).ne.0) then - ii(np)=0 - nstate(np)=ieor(nstate(np),1) - endif - go to 99 - endif - -! Back up - np=np-1 - if(np.lt.ntail .and. ii(np).ne.1) then -! Search the next best branch - ii(np)=ii(np)+1 - nstate(np)=ieor(nstate(np),1) - go to 99 - endif - go to 10 -99 continue - enddo - i=nbits*maxcycles - -100 metric=gamma(np) !Final path metric - -! Copy decoded data to user's buffer - nbytes=(nbits+7)/8 - np=7 - do j=1,nbytes-1 - i4a=nstate(np) - dat(j)=i4a - np=np+8 - enddo - dat(nbytes)=0 - - ncycles=i+1 - ierr=0 - if(i.ge.maxcycles*nbits) ierr=-1 - - return -end subroutine fano232 - -subroutine grid2deg(grid0,dlong,dlat) - -! Converts Maidenhead grid locator to degrees of West longitude -! and North latitude. - - character*6 grid0,grid - character*1 g1,g2,g3,g4,g5,g6 - - grid=grid0 - i=ichar(grid(5:5)) - if(grid(5:5).eq.' ' .or. i.le.64 .or. i.ge.128) grid(5:6)='mm' - - if(grid(1:1).ge.'a' .and. grid(1:1).le.'z') grid(1:1)= & - char(ichar(grid(1:1))+ichar('A')-ichar('a')) - if(grid(2:2).ge.'a' .and. grid(2:2).le.'z') grid(2:2)= & - char(ichar(grid(2:2))+ichar('A')-ichar('a')) - if(grid(5:5).ge.'A' .and. grid(5:5).le.'Z') grid(5:5)= & - char(ichar(grid(5:5))-ichar('A')+ichar('a')) - if(grid(6:6).ge.'A' .and. grid(6:6).le.'Z') grid(6:6)= & - char(ichar(grid(6:6))-ichar('A')+ichar('a')) - - g1=grid(1:1) - g2=grid(2:2) - g3=grid(3:3) - g4=grid(4:4) - g5=grid(5:5) - g6=grid(6:6) - - nlong = 180 - 20*(ichar(g1)-ichar('A')) - n20d = 2*(ichar(g3)-ichar('0')) - xminlong = 5*(ichar(g5)-ichar('a')+0.5) - dlong = nlong - n20d - xminlong/60.0 - nlat = -90+10*(ichar(g2)-ichar('A')) + ichar(g4)-ichar('0') - xminlat = 2.5*(ichar(g6)-ichar('a')+0.5) - dlat = nlat + xminlat/60.0 - - return -end subroutine grid2deg - -subroutine hash(string,len,ihash) - - parameter (MASK15=32767) - character*(*) string - integer*1 ic(12) - - do i=1,len - ic(i)=ichar(string(i:i)) - enddo - i=nhash(ic,len,146) - ihash=iand(i,MASK15) - -! print*,'C',ihash,len,string - return -end subroutine hash - -subroutine inter_mept(id,ndir) - -! Interleave (ndir=1) or de-interleave (ndir=-1) the array id. - - integer*1 id(0:161),itmp(0:161) - integer j0(0:161) - logical first - data first/.true./ - save - - if(first) then -! Compute the interleave table using bit reversal. - k=-1 - do i=0,255 - n=0 - ii=i - do j=0,7 - n=n+n - if(iand(ii,1).ne.0) n=n+1 - ii=ii/2 - enddo - if(n.le.161) then - k=k+1 - j0(k)=n - endif - enddo - first=.false. - endif - - if(ndir.eq.1) then - do i=0,161 - itmp(j0(i))=id(i) - enddo - else - do i=0,161 - itmp(i)=id(j0(i)) - enddo - endif - - do i=0,161 - id(i)=itmp(i) - enddo - - return -end subroutine inter_mept - -function nchar(c) - -! Convert ASCII number, letter, or space to 0-36 for callsign packing. - - character c*1 - data n/0/ !Silence compiler warning - - if(c.ge.'0' .and. c.le.'9') then - n=ichar(c)-ichar('0') - else if(c.ge.'A' .and. c.le.'Z') then - n=ichar(c)-ichar('A') + 10 - else if(c.ge.'a' .and. c.le.'z') then - n=ichar(c)-ichar('a') + 10 - else if(c.ge.' ') then - n=36 - else - Print*,'Invalid character in callsign ',c,' ',ichar(c) - stop - endif - nchar=n - - return -end function nchar - -subroutine pack50(n1,n2,dat) - - integer*1 dat(11),i1 - - i1=iand(ishft(n1,-20),255) !8 bits - dat(1)=i1 - i1=iand(ishft(n1,-12),255) !8 bits - dat(2)=i1 - i1=iand(ishft(n1, -4),255) !8 bits - dat(3)=i1 - i1=16*iand(n1,15)+iand(ishft(n2,-18),15) !4+4 bits - dat(4)=i1 - i1=iand(ishft(n2,-10),255) !8 bits - dat(5)=i1 - i1=iand(ishft(n2, -2),255) !8 bits - dat(6)=i1 - i1=64*iand(n2,3) !2 bits - dat(7)=i1 - dat(8)=0 - dat(9)=0 - dat(10)=0 - dat(11)=0 - - return -end subroutine pack50 - -subroutine packcall(callsign,ncall,text) - -! Pack a valid callsign into a 28-bit integer. - - parameter (NBASE=37*36*10*27*27*27) - character callsign*6,c*1,tmp*6,digit*10 - logical text - data digit/'0123456789'/ - - text=.false. - -! Work-around for Swaziland prefix: - if(callsign(1:4).eq.'3DA0') callsign='3D0'//callsign(5:6) - - if(callsign(1:3).eq.'CQ ') then - ncall=NBASE + 1 - if(callsign(4:4).ge.'0' .and. callsign(4:4).le.'9' .and. & - callsign(5:5).ge.'0' .and. callsign(5:5).le.'9' .and. & - callsign(6:6).ge.'0' .and. callsign(6:6).le.'9') then - nfreq=100*(ichar(callsign(4:4))-48) + & - 10*(ichar(callsign(5:5))-48) + & - ichar(callsign(6:6))-48 - ncall=NBASE + 3 + nfreq - endif - return - else if(callsign(1:4).eq.'QRZ ') then - ncall=NBASE + 2 - return - endif - - tmp=' ' - if(callsign(3:3).ge.'0' .and. callsign(3:3).le.'9') then - tmp=callsign - else if(callsign(2:2).ge.'0' .and. callsign(2:2).le.'9') then - if(callsign(6:6).ne.' ') then - text=.true. - return - endif - tmp=' '//callsign - else - text=.true. - return - endif - - do i=1,6 - c=tmp(i:i) - if(c.ge.'a' .and. c.le.'z') & - tmp(i:i)=char(ichar(c)-ichar('a')+ichar('A')) - enddo - - n1=0 - if((tmp(1:1).ge.'A'.and.tmp(1:1).le.'Z').or.tmp(1:1).eq.' ') n1=1 - if(tmp(1:1).ge.'0' .and. tmp(1:1).le.'9') n1=1 - n2=0 - if(tmp(2:2).ge.'A' .and. tmp(2:2).le.'Z') n2=1 - if(tmp(2:2).ge.'0' .and. tmp(2:2).le.'9') n2=1 - n3=0 - if(tmp(3:3).ge.'0' .and. tmp(3:3).le.'9') n3=1 - n4=0 - if((tmp(4:4).ge.'A'.and.tmp(4:4).le.'Z').or.tmp(4:4).eq.' ') n4=1 - n5=0 - if((tmp(5:5).ge.'A'.and.tmp(5:5).le.'Z').or.tmp(5:5).eq.' ') n5=1 - n6=0 - if((tmp(6:6).ge.'A'.and.tmp(6:6).le.'Z').or.tmp(6:6).eq.' ') n6=1 - - if(n1+n2+n3+n4+n5+n6 .ne. 6) then - text=.true. - return - endif - - ncall=nchar(tmp(1:1)) - ncall=36*ncall+nchar(tmp(2:2)) - ncall=10*ncall+nchar(tmp(3:3)) - ncall=27*ncall+nchar(tmp(4:4))-10 - ncall=27*ncall+nchar(tmp(5:5))-10 - ncall=27*ncall+nchar(tmp(6:6))-10 - - return -end subroutine packcall - -subroutine packgrid(grid,ng,text) - - parameter (NGBASE=180*180) - character*4 grid - logical text - - text=.false. - if(grid.eq.' ') go to 90 !Blank grid is OK - -! Test for numerical signal report, etc. - if(grid(1:1).eq.'-') then - n=10*(ichar(grid(2:2))-48) + ichar(grid(3:3)) - 48 - ng=NGBASE+1+n - go to 100 - else if(grid(1:2).eq.'R-') then - n=10*(ichar(grid(3:3))-48) + ichar(grid(4:4)) - 48 - if(n.eq.0) go to 90 - ng=NGBASE+31+n - go to 100 - else if(grid(1:2).eq.'RO') then - ng=NGBASE+62 - go to 100 - else if(grid(1:3).eq.'RRR') then - ng=NGBASE+63 - go to 100 - else if(grid(1:2).eq.'73') then - ng=NGBASE+64 - go to 100 - endif - - if(grid(1:1).lt.'A' .or. grid(1:1).gt.'R') text=.true. - if(grid(2:2).lt.'A' .or. grid(2:2).gt.'R') text=.true. - if(grid(3:3).lt.'0' .or. grid(3:3).gt.'9') text=.true. - if(grid(4:4).lt.'0' .or. grid(4:4).gt.'9') text=.true. - if(text) go to 100 - - call grid2deg(grid//'mm',dlong,dlat) - long=dlong - lat=dlat+ 90.0 - ng=((long+180)/2)*180 + lat - go to 100 - -90 ng=NGBASE + 1 - -100 return -end subroutine packgrid - -subroutine packpfx(call1,n1,ng,nadd) - - character*12 call1,call0 - character*3 pfx - logical text - - i1=index(call1,'/') - if(call1(i1+2:i1+2).eq.' ') then -! Single-character add-on suffix (maybe also fourth suffix letter?) - call0=call1(:i1-1) - call packcall(call0,n1,text) - nadd=1 - nc=ichar(call1(i1+1:i1+1)) - if(nc.ge.48 .and. nc.le.57) then - n=nc-48 - else if(nc.ge.65 .and. nc.le.90) then - n=nc-65+10 - else - n=38 - endif - nadd=1 - ng=60000-32768+n - else if(call1(i1+3:i1+3).eq.' ') then -! Two-character numerical suffix, /10 to /99 - call0=call1(:i1-1) - call packcall(call0,n1,text) - nadd=1 - n=10*(ichar(call1(i1+1:i1+1))-48) + ichar(call1(i1+2:i1+2)) - 48 - nadd=1 - ng=60000 + 26 + n - else -! Prefix of 1 to 3 characters - pfx=call1(:i1-1) - if(pfx(3:3).eq.' ') pfx=' '//pfx - if(pfx(3:3).eq.' ') pfx=' '//pfx - call0=call1(i1+1:) - call packcall(call0,n1,text) - - ng=0 - do i=1,3 - nc=ichar(pfx(i:i)) - if(nc.ge.48 .and. nc.le.57) then - n=nc-48 - else if(nc.ge.65 .and. nc.le.90) then - n=nc-65+10 - else - n=36 - endif - ng=37*ng + n - enddo - nadd=0 - if(ng.ge.32768) then - ng=ng-32768 - nadd=1 - endif - endif - - return -end subroutine packpfx - -subroutine unpack50(dat,n1,n2) - - integer*1 dat(11) - - i=dat(1) - i4=iand(i,255) - n1=ishft(i4,20) - i=dat(2) - i4=iand(i,255) - n1=n1 + ishft(i4,12) - i=dat(3) - i4=iand(i,255) - n1=n1 + ishft(i4,4) - i=dat(4) - i4=iand(i,255) - n1=n1 + iand(ishft(i4,-4),15) - n2=ishft(iand(i4,15),18) - i=dat(5) - i4=iand(i,255) - n2=n2 + ishft(i4,10) - i=dat(6) - i4=iand(i,255) - n2=n2 + ishft(i4,2) - i=dat(7) - i4=iand(i,255) - n2=n2 + iand(ishft(i4,-6),3) - - return -end subroutine unpack50 - -subroutine unpackcall(ncall,word) - - character word*12,c*37 - - data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/ - - n=ncall - word='......' - if(n.ge.262177560) go to 999 !Plain text message ... - i=mod(n,27)+11 - word(6:6)=c(i:i) - n=n/27 - i=mod(n,27)+11 - word(5:5)=c(i:i) - n=n/27 - i=mod(n,27)+11 - word(4:4)=c(i:i) - n=n/27 - i=mod(n,10)+1 - word(3:3)=c(i:i) - n=n/10 - i=mod(n,36)+1 - word(2:2)=c(i:i) - n=n/36 - i=n+1 - word(1:1)=c(i:i) - do i=1,4 - if(word(i:i).ne.' ') go to 10 - enddo - go to 999 -10 word=word(i:) - -999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:) - return -end subroutine unpackcall - -subroutine unpackgrid(ng,grid) - - parameter (NGBASE=180*180) - character grid*4,grid6*6,digit*10 - data digit/'0123456789'/ - - grid=' ' - if(ng.ge.32400) go to 10 - dlat=mod(ng,180)-90 - dlong=(ng/180)*2 - 180 + 2 - call deg2grid(dlong,dlat,grid6) - grid=grid6(1:4) !XXX explicitly truncate this -db - go to 100 - -10 n=ng-NGBASE-1 - if(n.ge.1 .and.n.le.30) then - grid(1:1)='-' - grid(2:2)=char(48+n/10) - grid(3:3)=char(48+mod(n,10)) - else if(n.ge.31 .and.n.le.60) then - n=n-30 - grid(1:2)='R-' - grid(3:3)=char(48+n/10) - grid(4:4)=char(48+mod(n,10)) - else if(n.eq.61) then - grid='RO' - else if(n.eq.62) then - grid='RRR' - else if(n.eq.63) then - grid='73' - endif - -100 return -end subroutine unpackgrid - -subroutine unpackpfx(ng,call1) - - character*12 call1 - character*3 pfx - - if(ng.lt.60000) then -! Add-on prefix of 1 to 3 characters - n=ng - do i=3,1,-1 - nc=mod(n,37) - if(nc.ge.0 .and. nc.le.9) then - pfx(i:i)=char(nc+48) - else if(nc.ge.10 .and. nc.le.35) then - pfx(i:i)=char(nc+55) - else - pfx(i:i)=' ' - endif - n=n/37 - enddo - call1=pfx//'/'//call1 - if(call1(1:1).eq.' ') call1=call1(2:) - if(call1(1:1).eq.' ') call1=call1(2:) - else -! Add-on suffix, one or teo characters - i1=index(call1,' ') - nc=ng-60000 - if(nc.ge.0 .and. nc.le.9) then - call1=call1(:i1-1)//'/'//char(nc+48) - else if(nc.ge.10 .and. nc.le.35) then - call1=call1(:i1-1)//'/'//char(nc+55) - else if(nc.ge.36 .and. nc.le.125) then - nc1=(nc-26)/10 - nc2=mod(nc-26,10) - call1=call1(:i1-1)//'/'//char(nc1+48)//char(nc2+48) - endif - endif - - return -end subroutine unpackpfx - -subroutine wqdecode(data0,message,ntype) - - parameter (N15=32768) - integer*1 data0(11) - character*22 message - character*12 callsign - character*3 cdbm - character grid4*4,grid6*6 - logical first - character*12 dcall(0:N15-1) - data first/.true./ - save first,dcall - -! May want to have a timeout (say, one hour?) on calls fetched -! from the hash table. - - if(first) then - dcall=' ' - first=.false. - endif - - message=' ' - call unpack50(data0,n1,n2) - call unpackcall(n1,callsign) - i1=index(callsign,' ') - call unpackgrid(n2/128,grid4) - ntype=iand(n2,127) -64 - -! Standard WSPR message (types 0 3 7 10 13 17 ... 60) - if(ntype.ge.0 .and. ntype.le.62) then - nu=mod(ntype,10) - if(nu.eq.0 .or. nu.eq.3 .or. nu.eq.7) then - write(cdbm,'(i3)'),ntype - if(cdbm(1:1).eq.' ') cdbm=cdbm(2:) - if(cdbm(1:1).eq.' ') cdbm=cdbm(2:) - message=callsign(1:i1)//grid4//' '//cdbm - call hash(callsign,i1-1,ih) - dcall(ih)=callsign(:i1) - else - nadd=nu - if(nu.gt.3) nadd=nu-3 - if(nu.gt.7) nadd=nu-7 - ng=n2/128 + 32768*(nadd-1) - call unpackpfx(ng,callsign) - ndbm=ntype-nadd - write(cdbm,'(i3)'),ndbm - if(cdbm(1:1).eq.' ') cdbm=cdbm(2:) - if(cdbm(1:1).eq.' ') cdbm=cdbm(2:) - i2=index(callsign,' ') - message=callsign(:i2)//cdbm - call hash(callsign,i2-1,ih) - dcall(ih)=callsign(:i2) - endif - else if(ntype.lt.0) then - ndbm=-(ntype+1) - grid6=callsign(6:6)//callsign(1:5) - ih=(n2-ntype-64)/128 - callsign=dcall(ih) - write(cdbm,'(i3)'),ndbm - if(cdbm(1:1).eq.' ') cdbm=cdbm(2:) - if(cdbm(1:1).eq.' ') cdbm=cdbm(2:) - i2=index(callsign,' ') - if(dcall(ih)(1:1).ne.' ') then - message='<'//callsign(:i2-1)//'> '//grid6//' '//cdbm - else - message='<...> '//grid6//' '//cdbm - endif - endif - - return -end subroutine wqdecode - -subroutine wqencode(msg,ntype,data0) - -! Parse and encode a WSPR message. - - parameter (MASK15=32767) - character*22 msg - character*12 call1,call2 - character grid4*4,grid6*6 - logical lbad1,lbad2 - integer*1 data0(11) - integer nu(0:9) - data nu/0,-1,1,0,-1,2,1,0,-1,1/ - -! Standard WSPR message (types 0 3 7 10 13 17 ... 60) - i1=index(msg,' ') - i2=index(msg,'/') - i3=index(msg,'<') - call1=msg(:i1-1) - if(i1.lt.3 .or. i1.gt.7 .or. i2.gt.0 .or. i3.gt.0) go to 10 - grid4=msg(i1+1:i1+4) - call packcall(call1,n1,lbad1) - call packgrid(grid4,ng,lbad2) - if(lbad1 .or. lbad2) go to 10 - ndbm=0 - read(msg(i1+5:),*) ndbm - if(ndbm.lt.0) ndbm=0 - if(ndbm.gt.60) ndbm=60 - ndbm=ndbm+nu(mod(ndbm,10)) - n2=128*ng + (ndbm+64) - call pack50(n1,n2,data0) - ntype=ndbm - go to 900 - -10 if(i2.ge.2 .and. i3.lt.1) then - call packpfx(call1,n1,ng,nadd) - ndbm=0 - read(msg(i1+1:),*) ndbm - if(ndbm.lt.0) ndbm=0 - if(ndbm.gt.60) ndbm=60 - ndbm=ndbm+nu(mod(ndbm,10)) - ntype=ndbm + 1 + nadd - n2=128*ng + ntype + 64 - call pack50(n1,n2,data0) - else if(i3.eq.1) then - i4=index(msg,'>') - call1=msg(2:i4-1) - call hash(call1,i4-2,ih) - grid6=msg(i1+1:i1+6) - call2=grid6(2:6)//grid6(1:1)//' ' - call packcall(call2,n1,lbad1) - ndbm=0 - read(msg(i1+8:),*) ndbm - if(ndbm.lt.0) ndbm=0 - if(ndbm.gt.60) ndbm=60 - ndbm=ndbm+nu(mod(ndbm,10)) - ntype=-(ndbm+1) - n2=128*ih + ntype + 64 - call pack50(n1,n2,data0) - endif - go to 900 - -900 continue - return -end subroutine wqencode diff --git a/lib/wsprcode/wsprcode.f90 b/lib/wsprcode/wsprcode.f90 deleted file mode 100644 index 88b5590..0000000 --- a/lib/wsprcode/wsprcode.f90 +++ /dev/null @@ -1,158 +0,0 @@ -!------------------------------------------------------------------------------- -! -! This file is part of the WSPR application, Weak Signal Propagation Reporter -! -! File Name: wsprcode.f90 -! Description: This program provides examples of the source encoding, -! convulsional error-control coding, bit and symbol ordering, -! and synchronizing information contained in WSPR messages. -! -! Copyright (C) 2001-2014 Joseph Taylor, K1JT -! License: GPL-3 -! -! 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. -! -!------------------------------------------------------------------------------- - -program wsprcode - - parameter (NSYM=162) - parameter (MAXSYM=176) - character*22 msg,msg2 - integer*1 data0(13) - integer*1 data1(13) - integer*1 dat(206) - integer*1 softsym(206) - -! Define the sync vector: - integer*1 sync(NSYM) - data sync/ & - 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/ - -! Metric table for decoding from soft symbols - integer mettab(0:255,0:1) - data mettab/ & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, & - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & - 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, & - 2, 2, 2, 2, 1, 1, 1, 1, 0, 0, & - -1, -1, -1, -2, -2, -3, -4, -4, -5, -6, & - -7, -7, -8, -9, -10, -11, -12, -12, -13, -14, & - -15, -16, -17, -17, -18, -19, -20, -21, -22, -22, & - -23, -24, -25, -26, -26, -27, -28, -29, -30, -30, & - -31, -32, -33, -33, -34, -35, -36, -36, -37, -38, & - -38, -39, -40, -41, -41, -42, -43, -43, -44, -45, & - -45, -46, -47, -47, -48, -49, -49, -50, -51, -51, & - -52, -53, -53, -54, -54, -55, -56, -56, -57, -57, & - -58, -59, -59, -60, -60, -61, -62, -62, -62, -63, & - -64, -64, -65, -65, -66, -67, -67, -67, -68, -69, & - -69, -70, -70, -71, -72, -72, -72, -72, -73, -74, & - -75, -75, -75, -77, -76, -76, -78, -78, -80, -81, & - -80, -79, -83, -82, -81, -82, -82, -83, -84, -84, & - -84, -87, -86, -87, -88, -89, -89, -89, -88, -87, & - -86, -87, -84, -84, -84, -83, -82, -82, -81, -82, & - -83, -79, -80, -81, -80, -78, -78, -76, -76, -77, & - -75, -75, -75, -74, -73, -72, -72, -72, -72, -71, & - -70, -70, -69, -69, -68, -67, -67, -67, -66, -65, & - -65, -64, -64, -63, -62, -62, -62, -61, -60, -60, & - -59, -59, -58, -57, -57, -56, -56, -55, -54, -54, & - -53, -53, -52, -51, -51, -50, -49, -49, -48, -47, & - -47, -46, -45, -45, -44, -43, -43, -42, -41, -41, & - -40, -39, -38, -38, -37, -36, -36, -35, -34, -33, & - -33, -32, -31, -30, -30, -29, -28, -27, -26, -26, & - -25, -24, -23, -22, -22, -21, -20, -19, -18, -17, & - -17, -16, -15, -14, -13, -12, -12, -11, -10, -9, & - -8, -7, -7, -6, -5, -4, -4, -3, -2, -2, & - -1, -1, -1, 0, 0, 1, 1, 1, 1, 2, & - 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, & - 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, & - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & - 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5/ - -! Get command-line argument(s) - nargs=iargc() - if(nargs.ne.1) then - print*,'Usage: WSPRcode "message"' - go to 999 - endif - call getarg(1,msg) !Get message from command line - write(*,1000) msg -1000 format('Message: ',a22) - - nbits=50+31 !User bits=50, constraint length=32 - nbytes=(nbits+7)/8 - ndelta=50 - limit=20000 - - data0=0 - call wqencode(msg,ntype0,data0) !Source encoding - write(*,1002) data0(1:7),data0(1:6),data0(7)/64 -1002 format(/'Source-encoded message, 50 bits:'/'Hex: ',7z3.2/ & - 'Binary: ',6b9.8,b3.2) - - call encode232(data0,nbytes,dat,MAXSYM) !Convolutional encoding - call inter_mept(dat,1) !Interleaving - - write(*,1004) -1004 format(/'Data symbols:') - write(*,1006) (dat(i),i=1,NSYM) -1006 format(5x,30i2) - - write(*,1008) -1008 format(/'Sync symbols:') - write(*,1006) (sync(i),i=1,NSYM) - - write(*,1010) -1010 format(/'Channel symbols:') - write(*,1006) (2*dat(i)+sync(i),i=1,NSYM) - - call inter_mept(dat,-1) !Remove interleaving - softsym=-dat !Simulate soft symbols - -! Call the sequential (Fano algorithm) decoder - call fano232(softsym,nbits,mettab,ndelta,limit,data1,ncycles,metric,nerr) - call wqdecode(data1,msg2,ntype1) - - write(*,1020) msg2,ntype1 -1020 format(/'Decoded message: ',a22,' ntype:',i3) - -999 end program wsprcode - -include 'wspr_old_subs.f90' - diff --git a/lib/wsprd/Makefile b/lib/wsprd/Makefile deleted file mode 100644 index 9204073..0000000 --- a/lib/wsprd/Makefile +++ /dev/null @@ -1,34 +0,0 @@ -CC = gcc -#CC = clang-3.5 -FC = gfortran - -CFLAGS= -I/usr/include -Wall -Wno-missing-braces -O3 -ffast-math -LDFLAGS = -L/usr/lib -FFLAGS = -O2 -Wall -Wno-conversion -LIBS = -lfftw3f -lm - -# Default rules -%.o: %.c $(DEPS) - ${CC} ${CFLAGS} -c $< -%.o: %.f - ${FC} ${FFLAGS} -c $< -%.o: %.F - ${FC} ${FFLAGS} -c $< -%.o: %.f90 - ${FC} ${FFLAGS} -c $< -%.o: %.F90 - ${FC} ${FFLAGS} -c $< - -all: wsprd wsprsim wsprd_exp - -DEPS = wsprsim_utils.h wsprd_utils.h fano.h jelinek.h nhash.h -OBJS1 = wsprd.o wsprsim_utils.o wsprd_utils.o tab.o fano.o jelinek.o nhash.o -wsprd: $(OBJS1) - $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS) - -OBJS2 = wsprsim.o wsprsim_utils.o wsprd_utils.o tab.o fano.o nhash.o -wsprsim: $(OBJS2) - $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS) - -clean: - rm *.o wsprd wsprsim diff --git a/lib/wsprd/Makefile.MinGW b/lib/wsprd/Makefile.MinGW deleted file mode 100644 index d9fd6b8..0000000 --- a/lib/wsprd/Makefile.MinGW +++ /dev/null @@ -1,30 +0,0 @@ -CC = gcc -#CC = clang -FC = gfortran - -FFLAGS = -O2 -Wall -Wno-conversion -CFLAGS= -Wall -Wno-missing-braces -O2 -#LDFLAGS = -L/JTSDK/fftw3f -LIBS = c:/JTSDK/fftw3f/libfftw3-3.dll -lm - -# Default rules -%.o: %.c $(DEPS) - ${CC} ${CFLAGS} -c $< -%.o: %.f - ${FC} ${FFLAGS} -c $< -%.o: %.F - ${FC} ${FFLAGS} -c $< -%.o: %.f90 - ${FC} ${FFLAGS} -c $< -%.o: %.F90 - ${FC} ${FFLAGS} -c $< - -all: wsprd - -DEPS = fano.h -OBJS1 = wsprd.o wsprd_utils.o fano.o tab.o nhash.o -wsprd: $(OBJS1) - $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS) - -clean: - rm *.o wsprd diff --git a/lib/wsprd/Makefile.win32 b/lib/wsprd/Makefile.win32 deleted file mode 100644 index 8d600d7..0000000 --- a/lib/wsprd/Makefile.win32 +++ /dev/null @@ -1,40 +0,0 @@ -CC = gcc -FC = gfortran - -FFLAGS = -O2 -Wall -Wno-conversion -CFLAGS= -I/JTSDK/fftw3f -Wall -Wno-missing-braces -O2 -#LDFLAGS = -L/JTSDK/fftw3f -Wl,--stack,4000000 -LDFLAGS = -L/JTSDK/fftw3f -LIBS = c:/JTSDK/fftw3f/libfftw3-3.dll -lm - -# Default rules -%.o: %.c $(DEPS) - ${CC} ${CFLAGS} -c $< -%.o: %.f - ${FC} ${FFLAGS} -c $< -%.o: %.F - ${FC} ${FFLAGS} -c $< -%.o: %.f90 - ${FC} ${FFLAGS} -c $< -%.o: %.F90 - ${FC} ${FFLAGS} -c $< - -all: wsprd.exe wsprsim.exe wsprd_exp.exe - -DEPS = wsprsim_utils.h wsprd_utils.h fano.h jelinek.h nhash.h - -OBJS1 = wsprd.o wsprsim_utils.o wsprd_utils.o tab.o fano.o jelinek.o nhash.o -wsprd.exe: $(OBJS1) - $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS) - -OBJS2 = wsprsim.o wsprsim_utils.o wsprd_utils.o tab.o fano.o nhash.o -wsprsim.exe: $(OBJS2) - $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS) - -OBJS3 = wsprd_exp.o wsprsim_utils.o wsprd_utils.o tab.o fano.o jelinek.o \ - nhash.o -wsprd_exp.exe: $(OBJS3) - $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS) - -clean: - rm *.o wsprd.exe wsprsim.exe wsprd_exp.exe diff --git a/lib/wsprd/README b/lib/wsprd/README deleted file mode 100755 index 4579d2d..0000000 --- a/lib/wsprd/README +++ /dev/null @@ -1,54 +0,0 @@ -wsprd is a decoder for K1JT's Weak Signal Propagation Reporter (WSPR) mode. - -The program is written in C and is a command-line program that reads from a -.c2 file or .wav file and writes output to the console. It is used by WSJT-X -for wspr-mode decoding. - -USAGE: - wsprd [options...] infile - -OPTIONS: - -a path to writeable data files, default="." - -c write .c2 file at the end of the first pass - -e x (x is transceiver dial frequency error in Hz) - -f x (x is transceiver dial frequency in MHz) - -H do not use (or update) the hash table - -m decode wspr-15 .wav file - -q quick mode - doesn't dig deep for weak signals - -s single pass mode, no subtraction (same as original wsprd) - -v verbose mode (shows dupes) - -w wideband mode - decode signals within +/- 150 Hz of center - -z x (x is fano metric table bias, default is 0.42) - -infile can be either .wav or .c2 - -e.g. -./wsprd -wf 14.0956 140709_2258.wav - -Note that for .c2 files, the frequency within the file overrides the command -line value. - -FEATURES: -By default, wsprd reports signals that are within +/- 110 Hz of the -subband center frequency. The wideband option (-w) extends this to +/- 150 Hz. - -wsprd maintains a hashtable and will decode all three types of wspr -messages. An option (-H) is available to turn off use of the hashtable. - -The symbols are decoded using Phil Karn's sequential decoder routine, -fano.c. - -NOTES: -This program attempts to maximize the number of successful decodes per transmit -interval by trying to decode virtually every peak in the averaged spectrum. -The program also implements two-pass decoding, whereby signals that are successfully -decoded are subtracted one-by-one during the first decoding pass. Then, the -decoder is run again. In many cases the subtraction process will uncover signals -that can then be successfully decoded on the second pass. - -There will be occasional duplicate decodes when two closely spaced -peaks come from the same signal. The program removes dupes based on callsign -and frequency. Two decodes that have the same callsign and estimated frequencies -that are within 1 Hz will be treated as decodes of the same signal. This -dupechecking is turned off with the -v flag. - diff --git a/lib/wsprd/WSPRcode.f90 b/lib/wsprd/WSPRcode.f90 deleted file mode 100644 index 1ec9171..0000000 --- a/lib/wsprd/WSPRcode.f90 +++ /dev/null @@ -1,132 +0,0 @@ -program wsprcode - -! This program provides examples of the source encoding, convolutional -! error-control coding, bit and symbol ordering, and synchronizing -! information contained in WSPR messages. - - parameter (NSYM=162) - parameter (MAXSYM=176) - character*22 msg,msg2 - integer*1 data0(7) - integer*1 data1(7) - integer*1 dat(NSYM) - integer*1 softsym(NSYM) - -! Define the sync vector: - integer*1 sync(NSYM) - data sync/ & - 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/ - -! Metric table for decoding from soft symbols - integer mettab(0:255,0:1) - data mettab/ & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, & - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & - 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, & - 2, 2, 2, 2, 1, 1, 1, 1, 0, 0, & - -1, -1, -1, -2, -2, -3, -4, -4, -5, -6, & - -7, -7, -8, -9, -10, -11, -12, -12, -13, -14, & - -15, -16, -17, -17, -18, -19, -20, -21, -22, -22, & - -23, -24, -25, -26, -26, -27, -28, -29, -30, -30, & - -31, -32, -33, -33, -34, -35, -36, -36, -37, -38, & - -38, -39, -40, -41, -41, -42, -43, -43, -44, -45, & - -45, -46, -47, -47, -48, -49, -49, -50, -51, -51, & - -52, -53, -53, -54, -54, -55, -56, -56, -57, -57, & - -58, -59, -59, -60, -60, -61, -62, -62, -62, -63, & - -64, -64, -65, -65, -66, -67, -67, -67, -68, -69, & - -69, -70, -70, -71, -72, -72, -72, -72, -73, -74, & - -75, -75, -75, -77, -76, -76, -78, -78, -80, -81, & - -80, -79, -83, -82, -81, -82, -82, -83, -84, -84, & - -84, -87, -86, -87, -88, -89, -89, -89, -88, -87, & - -86, -87, -84, -84, -84, -83, -82, -82, -81, -82, & - -83, -79, -80, -81, -80, -78, -78, -76, -76, -77, & - -75, -75, -75, -74, -73, -72, -72, -72, -72, -71, & - -70, -70, -69, -69, -68, -67, -67, -67, -66, -65, & - -65, -64, -64, -63, -62, -62, -62, -61, -60, -60, & - -59, -59, -58, -57, -57, -56, -56, -55, -54, -54, & - -53, -53, -52, -51, -51, -50, -49, -49, -48, -47, & - -47, -46, -45, -45, -44, -43, -43, -42, -41, -41, & - -40, -39, -38, -38, -37, -36, -36, -35, -34, -33, & - -33, -32, -31, -30, -30, -29, -28, -27, -26, -26, & - -25, -24, -23, -22, -22, -21, -20, -19, -18, -17, & - -17, -16, -15, -14, -13, -12, -12, -11, -10, -9, & - -8, -7, -7, -6, -5, -4, -4, -3, -2, -2, & - -1, -1, -1, 0, 0, 1, 1, 1, 1, 2, & - 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, & - 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, & - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & - 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & - 5, 5/ - -! Get command-line argument(s) - nargs=iargc() - if(nargs.ne.1) then - print*,'Usage: WSPRcode "message"' - go to 999 - endif - call getarg(1,msg) !Get message from command line - write(*,1000) msg -1000 format('Message: ',a22) - - nbits=50+31 !User bits=50, constraint length=32 - nbytes=(nbits+7)/8 - ndelta=50 - limit=20000 - - data0=0 - call wqencode(msg,ntype0,data0) !Source encoding - write(*,1002) data0 -1002 format(/'Source-encoded message (50 bits, hex):',7z3.2) - - call encode232(data0,nbytes,dat,MAXSYM) !Convolutional encoding - call inter_mept(dat,1) !Interleaving - - write(*,1004) -1004 format(/'Data symbols:') - write(*,1006) (dat(i),i=1,NSYM) -1006 format(5x,30i2) - - write(*,1008) -1008 format(/'Sync symbols:') - write(*,1006) (sync(i),i=1,NSYM) - - write(*,1010) -1010 format(/'Channel symbols:') - write(*,1006) (2*dat(i)+sync(i),i=1,NSYM) - - call inter_mept(dat,-1) !Remove interleaving - softsym=-dat !Simulate soft symbols - -! Call the sequential (Fano algorithm) decoder - call fano232(softsym,nbits,mettab,ndelta,limit,data1,ncycles,metric,nerr) - call wqdecode(data1,msg2,ntype1) - - write(*,1020) ntype1 -1020 format(/'Message type: ',i7) - write(*,1030) msg2 -1030 format('Decoded message: ',a22) - -999 end program wsprcode diff --git a/lib/wsprd/fano.c b/lib/wsprd/fano.c deleted file mode 100644 index 8fe1a64..0000000 --- a/lib/wsprd/fano.c +++ /dev/null @@ -1,239 +0,0 @@ -/* - This file is part of wsprd. - - File name: fano.c - - Description: Soft decision Fano sequential decoder for K=32 r=1/2 - convolutional code. - - Copyright 1994, Phil Karn, KA9Q - Minor modifications by Joe Taylor, K1JT -*/ - -#define LL 1 // Select Layland-Lushbaugh code -#include -#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. - */ -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; -} - -/* Decode packet with the Fano algorithm. - * Return 0 on success, -1 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++) { - 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 -1; // Decoder timed out - return 0; // Successful completion -} diff --git a/lib/wsprd/fano.h b/lib/wsprd/fano.h deleted file mode 100644 index 3290a09..0000000 --- a/lib/wsprd/fano.h +++ /dev/null @@ -1,39 +0,0 @@ -/* - 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/lib/wsprd/fftw3.h b/lib/wsprd/fftw3.h deleted file mode 100644 index 58a2c73..0000000 --- a/lib/wsprd/fftw3.h +++ /dev/null @@ -1,410 +0,0 @@ -/* - * Copyright (c) 2003, 2007-11 Matteo Frigo - * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology - * - * The following statement of license applies *only* to this header file, - * and *not* to the other files distributed with FFTW or derived therefrom: - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS - * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE - * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS - * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -/***************************** NOTE TO USERS ********************************* - * - * THIS IS A HEADER FILE, NOT A MANUAL - * - * If you want to know how to use FFTW, please read the manual, - * online at http://www.fftw.org/doc/ and also included with FFTW. - * For a quick start, see the manual's tutorial section. - * - * (Reading header files to learn how to use a library is a habit - * stemming from code lacking a proper manual. Arguably, it's a - * *bad* habit in most cases, because header files can contain - * interfaces that are not part of the public, stable API.) - * - ****************************************************************************/ - -#ifndef FFTW3_H -#define FFTW3_H - -#include - -#ifdef __cplusplus -extern "C" -{ -#endif /* __cplusplus */ - -/* If is included, use the C99 complex type. Otherwise - define a type bit-compatible with C99 complex */ -#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) -# define FFTW_DEFINE_COMPLEX(R, C) typedef R _Complex C -#else -# define FFTW_DEFINE_COMPLEX(R, C) typedef R C[2] -#endif - -#define FFTW_CONCAT(prefix, name) prefix ## name -#define FFTW_MANGLE_DOUBLE(name) FFTW_CONCAT(fftw_, name) -#define FFTW_MANGLE_FLOAT(name) FFTW_CONCAT(fftwf_, name) -#define FFTW_MANGLE_LONG_DOUBLE(name) FFTW_CONCAT(fftwl_, name) -#define FFTW_MANGLE_QUAD(name) FFTW_CONCAT(fftwq_, name) - -/* IMPORTANT: for Windows compilers, you should add a line - #define FFTW_DLL - here and in kernel/ifftw.h if you are compiling/using FFTW as a - DLL, in order to do the proper importing/exporting, or - alternatively compile with -DFFTW_DLL or the equivalent - command-line flag. This is not necessary under MinGW/Cygwin, where - libtool does the imports/exports automatically. */ -#if defined(FFTW_DLL) && (defined(_WIN32) || defined(__WIN32__)) - /* annoying Windows syntax for shared-library declarations */ -# if defined(COMPILING_FFTW) /* defined in api.h when compiling FFTW */ -# define FFTW_EXTERN extern __declspec(dllexport) -# else /* user is calling FFTW; import symbol */ -# define FFTW_EXTERN extern __declspec(dllimport) -# endif -#else -# define FFTW_EXTERN extern -#endif - -enum fftw_r2r_kind_do_not_use_me { - FFTW_R2HC=0, FFTW_HC2R=1, FFTW_DHT=2, - FFTW_REDFT00=3, FFTW_REDFT01=4, FFTW_REDFT10=5, FFTW_REDFT11=6, - FFTW_RODFT00=7, FFTW_RODFT01=8, FFTW_RODFT10=9, FFTW_RODFT11=10 -}; - -struct fftw_iodim_do_not_use_me { - int n; /* dimension size */ - int is; /* input stride */ - int os; /* output stride */ -}; - -#include /* for ptrdiff_t */ -struct fftw_iodim64_do_not_use_me { - ptrdiff_t n; /* dimension size */ - ptrdiff_t is; /* input stride */ - ptrdiff_t os; /* output stride */ -}; - -typedef void (*fftw_write_char_func_do_not_use_me)(char c, void *); -typedef int (*fftw_read_char_func_do_not_use_me)(void *); - -/* - huge second-order macro that defines prototypes for all API - functions. We expand this macro for each supported precision - - X: name-mangling macro - R: real data type - C: complex data type -*/ - -#define FFTW_DEFINE_API(X, R, C) \ - \ -FFTW_DEFINE_COMPLEX(R, C); \ - \ -typedef struct X(plan_s) *X(plan); \ - \ -typedef struct fftw_iodim_do_not_use_me X(iodim); \ -typedef struct fftw_iodim64_do_not_use_me X(iodim64); \ - \ -typedef enum fftw_r2r_kind_do_not_use_me X(r2r_kind); \ - \ -typedef fftw_write_char_func_do_not_use_me X(write_char_func); \ -typedef fftw_read_char_func_do_not_use_me X(read_char_func); \ - \ -FFTW_EXTERN void X(execute)(const X(plan) p); \ - \ -FFTW_EXTERN X(plan) X(plan_dft)(int rank, const int *n, \ - C *in, C *out, int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_1d)(int n, C *in, C *out, int sign, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_2d)(int n0, int n1, \ - C *in, C *out, int sign, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_3d)(int n0, int n1, int n2, \ - C *in, C *out, int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_many_dft)(int rank, const int *n, \ - int howmany, \ - C *in, const int *inembed, \ - int istride, int idist, \ - C *out, const int *onembed, \ - int ostride, int odist, \ - int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_dft)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - C *in, C *out, \ - int sign, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru_split_dft)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *ri, R *ii, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_dft)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - C *in, C *out, \ - int sign, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru64_split_dft)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *ri, R *ii, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN void X(execute_dft)(const X(plan) p, C *in, C *out); \ -FFTW_EXTERN void X(execute_split_dft)(const X(plan) p, R *ri, R *ii, \ - R *ro, R *io); \ - \ -FFTW_EXTERN X(plan) X(plan_many_dft_r2c)(int rank, const int *n, \ - int howmany, \ - R *in, const int *inembed, \ - int istride, int idist, \ - C *out, const int *onembed, \ - int ostride, int odist, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_r2c)(int rank, const int *n, \ - R *in, C *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_r2c_1d)(int n,R *in,C *out,unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_r2c_2d)(int n0, int n1, \ - R *in, C *out, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_r2c_3d)(int n0, int n1, \ - int n2, \ - R *in, C *out, unsigned flags); \ - \ - \ -FFTW_EXTERN X(plan) X(plan_many_dft_c2r)(int rank, const int *n, \ - int howmany, \ - C *in, const int *inembed, \ - int istride, int idist, \ - R *out, const int *onembed, \ - int ostride, int odist, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_c2r)(int rank, const int *n, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_c2r_1d)(int n,C *in,R *out,unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_c2r_2d)(int n0, int n1, \ - C *in, R *out, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_c2r_3d)(int n0, int n1, \ - int n2, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_dft_r2c)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, C *out, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru_dft_c2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - C *in, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_split_dft_r2c)( \ - int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, R *ro, R *io, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru_split_dft_c2r)( \ - int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *ri, R *ii, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_dft_r2c)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, C *out, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru64_dft_c2r)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - C *in, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_split_dft_r2c)( \ - int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, R *ro, R *io, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru64_split_dft_c2r)( \ - int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *ri, R *ii, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN void X(execute_dft_r2c)(const X(plan) p, R *in, C *out); \ -FFTW_EXTERN void X(execute_dft_c2r)(const X(plan) p, C *in, R *out); \ - \ -FFTW_EXTERN void X(execute_split_dft_r2c)(const X(plan) p, \ - R *in, R *ro, R *io); \ -FFTW_EXTERN void X(execute_split_dft_c2r)(const X(plan) p, \ - R *ri, R *ii, R *out); \ - \ -FFTW_EXTERN X(plan) X(plan_many_r2r)(int rank, const int *n, \ - int howmany, \ - R *in, const int *inembed, \ - int istride, int idist, \ - R *out, const int *onembed, \ - int ostride, int odist, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_r2r)(int rank, const int *n, R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_r2r_1d)(int n, R *in, R *out, \ - X(r2r_kind) kind, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_r2r_2d)(int n0, int n1, R *in, R *out, \ - X(r2r_kind) kind0, X(r2r_kind) kind1, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_r2r_3d)(int n0, int n1, int n2, \ - R *in, R *out, X(r2r_kind) kind0, \ - X(r2r_kind) kind1, X(r2r_kind) kind2, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_r2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_r2r)(int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN void X(execute_r2r)(const X(plan) p, R *in, R *out); \ - \ -FFTW_EXTERN void X(destroy_plan)(X(plan) p); \ -FFTW_EXTERN void X(forget_wisdom)(void); \ -FFTW_EXTERN void X(cleanup)(void); \ - \ -FFTW_EXTERN void X(set_timelimit)(double t); \ - \ -FFTW_EXTERN void X(plan_with_nthreads)(int nthreads); \ -FFTW_EXTERN int X(init_threads)(void); \ -FFTW_EXTERN void X(cleanup_threads)(void); \ - \ -FFTW_EXTERN int X(export_wisdom_to_filename)(const char *filename); \ -FFTW_EXTERN void X(export_wisdom_to_file)(FILE *output_file); \ -FFTW_EXTERN char *X(export_wisdom_to_string)(void); \ -FFTW_EXTERN void X(export_wisdom)(X(write_char_func) write_char, \ - void *data); \ -FFTW_EXTERN int X(import_system_wisdom)(void); \ -FFTW_EXTERN int X(import_wisdom_from_filename)(const char *filename); \ -FFTW_EXTERN int X(import_wisdom_from_file)(FILE *input_file); \ -FFTW_EXTERN int X(import_wisdom_from_string)(const char *input_string); \ -FFTW_EXTERN int X(import_wisdom)(X(read_char_func) read_char, void *data); \ - \ -FFTW_EXTERN void X(fprint_plan)(const X(plan) p, FILE *output_file); \ -FFTW_EXTERN void X(print_plan)(const X(plan) p); \ - \ -FFTW_EXTERN void *X(malloc)(size_t n); \ -FFTW_EXTERN R *X(alloc_real)(size_t n); \ -FFTW_EXTERN C *X(alloc_complex)(size_t n); \ -FFTW_EXTERN void X(free)(void *p); \ - \ -FFTW_EXTERN void X(flops)(const X(plan) p, \ - double *add, double *mul, double *fmas); \ -FFTW_EXTERN double X(estimate_cost)(const X(plan) p); \ -FFTW_EXTERN double X(cost)(const X(plan) p); \ - \ -FFTW_EXTERN const char X(version)[]; \ -FFTW_EXTERN const char X(cc)[]; \ -FFTW_EXTERN const char X(codelet_optim)[]; - - -/* end of FFTW_DEFINE_API macro */ - -FFTW_DEFINE_API(FFTW_MANGLE_DOUBLE, double, fftw_complex) -FFTW_DEFINE_API(FFTW_MANGLE_FLOAT, float, fftwf_complex) -FFTW_DEFINE_API(FFTW_MANGLE_LONG_DOUBLE, long double, fftwl_complex) - -/* __float128 (quad precision) is a gcc extension on i386, x86_64, and ia64 - for gcc >= 4.6 (compiled in FFTW with --enable-quad-precision) */ -#if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) \ - && !(defined(__ICC) || defined(__INTEL_COMPILER)) \ - && (defined(__i386__) || defined(__x86_64__) || defined(__ia64__)) -# if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) -/* note: __float128 is a typedef, which is not supported with the _Complex - keyword in gcc, so instead we use this ugly __attribute__ version. - However, we can't simply pass the __attribute__ version to - FFTW_DEFINE_API because the __attribute__ confuses gcc in pointer - types. Hence redefining FFTW_DEFINE_COMPLEX. Ugh. */ -# undef FFTW_DEFINE_COMPLEX -# define FFTW_DEFINE_COMPLEX(R, C) typedef _Complex float __attribute__((mode(TC))) C -# endif -FFTW_DEFINE_API(FFTW_MANGLE_QUAD, __float128, fftwq_complex) -#endif - -#define FFTW_FORWARD (-1) -#define FFTW_BACKWARD (+1) - -#define FFTW_NO_TIMELIMIT (-1.0) - -/* documented flags */ -#define FFTW_MEASURE (0U) -#define FFTW_DESTROY_INPUT (1U << 0) -#define FFTW_UNALIGNED (1U << 1) -#define FFTW_CONSERVE_MEMORY (1U << 2) -#define FFTW_EXHAUSTIVE (1U << 3) /* NO_EXHAUSTIVE is default */ -#define FFTW_PRESERVE_INPUT (1U << 4) /* cancels FFTW_DESTROY_INPUT */ -#define FFTW_PATIENT (1U << 5) /* IMPATIENT is default */ -#define FFTW_ESTIMATE (1U << 6) -#define FFTW_WISDOM_ONLY (1U << 21) - -/* undocumented beyond-guru flags */ -#define FFTW_ESTIMATE_PATIENT (1U << 7) -#define FFTW_BELIEVE_PCOST (1U << 8) -#define FFTW_NO_DFT_R2HC (1U << 9) -#define FFTW_NO_NONTHREADED (1U << 10) -#define FFTW_NO_BUFFERING (1U << 11) -#define FFTW_NO_INDIRECT_OP (1U << 12) -#define FFTW_ALLOW_LARGE_GENERIC (1U << 13) /* NO_LARGE_GENERIC is default */ -#define FFTW_NO_RANK_SPLITS (1U << 14) -#define FFTW_NO_VRANK_SPLITS (1U << 15) -#define FFTW_NO_VRECURSE (1U << 16) -#define FFTW_NO_SIMD (1U << 17) -#define FFTW_NO_SLOW (1U << 18) -#define FFTW_NO_FIXED_RADIX_LARGE_N (1U << 19) -#define FFTW_ALLOW_PRUNING (1U << 20) - -#ifdef __cplusplus -} /* extern "C" */ -#endif /* __cplusplus */ - -#endif /* FFTW3_H */ diff --git a/lib/wsprd/genmet.f90 b/lib/wsprd/genmet.f90 deleted file mode 100644 index e5dcdab..0000000 --- a/lib/wsprd/genmet.f90 +++ /dev/null @@ -1,50 +0,0 @@ -program genmet - - character*12 arg - integer hist(-128:128) - lim(x)=min(127,max(-128,nint(scale*x))) - - nargs=iargc() - if(nargs.ne.4) then - print*,'Usage: genmet bw scale snr iters' - print*,'Example: genmet 1.46 20 -24 1000000' - go to 999 - endif - call getarg(1,arg) - read(arg,*) bw - call getarg(2,arg) - read(arg,*) scale - call getarg(3,arg) - read(arg,*) snr - call getarg(4,arg) - read(arg,*) iters - - hist=0 - s=sqrt(2500.0/bw) * 10.0**(0.05*snr) - fac=1.0/sqrt(2.0) - do iter=1,iters - x1=fac*gran() - y1=fac*gran() - x0=fac*gran() - y0=fac*gran() - r=(x1+s)**2 + y1*y1 - x0*x0 - y0*y0 - hist(lim(r))=hist(lim(r))+1 - enddo - - xln2=log(2.0) - do i=-128,127 - p1=hist(i)/dfloat(iters) - j=-i - if(j.gt.127) j=127 - p0=hist(j)/dfloat(iters) - xlhd0=log(max(0.001,2.0*p0/(p0+p1)))/xln2 - xlhd1=log(max(0.001,2.0*p1/(p0+p1)))/xln2 - write(13,1010) i/scale,hist(i)/dfloat(iters) -1010 format(f8.3,f12.9) - write(14,1012) i+128,xlhd0,xlhd1 -1012 format(i4,2f8.3) - enddo - -999 end program genmet - - diff --git a/lib/wsprd/gran.c b/lib/wsprd/gran.c deleted file mode 100644 index 24b9865..0000000 --- a/lib/wsprd/gran.c +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include - -/* Generate gaussian random float with mean=0 and std_dev=1 */ -float gran_() -{ - float fac,rsq,v1,v2; - static float gset; - static int iset; - - if(iset){ - /* Already got one */ - iset = 0; - return gset; - } - /* Generate two evenly distributed numbers between -1 and +1 - * that are inside the unit circle - */ - do { - v1 = 2.0 * (float)rand() / RAND_MAX - 1; - v2 = 2.0 * (float)rand() / RAND_MAX - 1; - rsq = v1*v1 + v2*v2; - } while(rsq >= 1.0 || rsq == 0.0); - fac = sqrt(-2.0*log(rsq)/rsq); - gset = v1*fac; - iset++; - return v2*fac; -} diff --git a/lib/wsprd/jelinek.c b/lib/wsprd/jelinek.c deleted file mode 100644 index 8fe0d2f..0000000 --- a/lib/wsprd/jelinek.c +++ /dev/null @@ -1,163 +0,0 @@ -/* - 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 0 on success, -1 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 -1; - } - return 0; //success -} diff --git a/lib/wsprd/jelinek.h b/lib/wsprd/jelinek.h deleted file mode 100644 index 64ecdb9..0000000 --- a/lib/wsprd/jelinek.h +++ /dev/null @@ -1,26 +0,0 @@ -#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; -}; - -struct snode *stack; - -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/lib/wsprd/metric_tables.c b/lib/wsprd/metric_tables.c deleted file mode 100644 index 973651d..0000000 --- a/lib/wsprd/metric_tables.c +++ /dev/null @@ -1,138 +0,0 @@ -/******************************************************************************* -* 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/lib/wsprd/mettab.c b/lib/wsprd/mettab.c deleted file mode 100644 index f188c0c..0000000 --- a/lib/wsprd/mettab.c +++ /dev/null @@ -1,76 +0,0 @@ -/* - This file is part of wsprd. - - File name: mettab.c - Description: Metric table for sequential Fano decoder. - - Copyright 2008-2015, 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, see . -*/ - -int mettab[2][256]={ - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, - 2, 2, 2, 2, 1, 1, 1, 1, 0, 0, - -1, -1, -1, -2, -2, -3, -4, -4, -5, -6, - -7, -7, -8, -9, -10, -11, -12, -12, -13, -14, - -15, -16, -17, -17, -18, -19, -20, -21, -22, -22, - -23, -24, -25, -26, -26, -27, -28, -29, -30, -30, - -31, -32, -33, -33, -34, -35, -36, -36, -37, -38, - -38, -39, -40, -41, -41, -42, -43, -43, -44, -45, - -45, -46, -47, -47, -48, -49, -49, -50, -51, -51, - -52, -53, -53, -54, -54, -55, -56, -56, -57, -57, - -58, -59, -59, -60, -60, -61, -62, -62, -62, -63, - -64, -64, -65, -65, -66, -67, -67, -67, -68, -69, - -69, -70, -70, -71, -72, -72, -72, -72, -73, -74, - -75, -75, -75, -77, -76, -76, -78, -78, -80, -81, - -80, -79, -83, -82, -81, -82, -82, -83, -84, -84, - -84, -87, -86, -87, -88,-105, -94,-105, -88, -87, - -86, -87, -84, -84, -84, -83, -82, -82, -81, -82, - -83, -79, -80, -81, -80, -78, -78, -76, -76, -77, - -75, -75, -75, -74, -73, -72, -72, -72, -72, -71, - -70, -70, -69, -69, -68, -67, -67, -67, -66, -65, - -65, -64, -64, -63, -62, -62, -62, -61, -60, -60, - -59, -59, -58, -57, -57, -56, -56, -55, -54, -54, - -53, -53, -52, -51, -51, -50, -49, -49, -48, -47, - -47, -46, -45, -45, -44, -43, -43, -42, -41, -41, - -40, -39, -38, -38, -37, -36, -36, -35, -34, -33, - -33, -32, -31, -30, -30, -29, -28, -27, -26, -26, - -25, -24, -23, -22, -22, -21, -20, -19, -18, -17, - -17, -16, -15, -14, -13, -12, -12, -11, -10, -9, - -8, -7, -7, -6, -5, -4, -4, -3, -2, -2, - -1, -1, -1, 0, 0, 1, 1, 1, 1, 2, - 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, - 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5 }; diff --git a/lib/wsprd/t2.f90 b/lib/wsprd/t2.f90 deleted file mode 100644 index 0f214db..0000000 --- a/lib/wsprd/t2.f90 +++ /dev/null @@ -1,11 +0,0 @@ -program t2 - - df=375.0/65536.0 - do i=1,65536 - w=1.0/(1.0 + ((i-32768)/26214.0)**20) - f=(i-32768)*df - write(13,1010) f,w -1010 format(2f15.6) - enddo - -end program t2 diff --git a/lib/wsprd/tab.c b/lib/wsprd/tab.c deleted file mode 100644 index e330c75..0000000 --- a/lib/wsprd/tab.c +++ /dev/null @@ -1,41 +0,0 @@ -/* - 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/lib/wsprd/test_wspr.f90 b/lib/wsprd/test_wspr.f90 deleted file mode 100644 index 86e2918..0000000 --- a/lib/wsprd/test_wspr.f90 +++ /dev/null @@ -1,61 +0,0 @@ -program test_wspr - -! This program provides examples of the source encoding, convolutional -! error-control coding, bit and symbol ordering, and synchronizing -! information contained in WSPR messages. - - character*22 msg,msg2 - character*23 msg3 - character*1 err2,err3 - integer*1 data0(11) - logical lfile - -! Get command-line argument(s) - nargs=iargc() - if(nargs.ne.1) then - print*,'Usage: test_wspr "message"' - go to 999 - endif - call getarg(1,msg) !Get message from command line - call unpk(data0,1,msg3) !Read the C hashtable - lfile=msg(1:2).eq."-t" - if(lfile) open(10,file="messages.txt",status="old") - - do imsg=1,999 - if(lfile) read(10,1001,end=900) msg -1001 format(a22) - - data0=0 - call wqencode(msg,ntype0,data0) !Source encoding -! write(*,1002) data0(1:7) -!1002 format('Source-encoded message (50 bits, hex):',7z3.2) -! data0(8:11)=0 - - call wqdecode(data0,msg2,ntype1) - -! write(*,1020) ntype1 -!1020 format('Message type: ',i7) -! write(*,1030) msg2 -!1030 format('Decoded message: ',a22) - - call unpk(data0,0,msg3) - do i=1,23 - if(ichar(msg3(i:i)).eq.0) then - msg3(i:)=" " - exit - endif - enddo - - err2=' ' - err3=' ' - if(msg2.ne.msg) err2='*' - if(msg3.ne.msg) err3='*' - - write(*,1040) msg,err2,msg2,err3,msg3 -1040 format(a22,1x,a1,1x,a22,1x,a1,1x,a22) - if(.not.lfile) exit - enddo -900 call unpk(data0,2,msg3) - - -999 end program test_wspr diff --git a/lib/wsprd/unpk.c.obsolete b/lib/wsprd/unpk.c.obsolete deleted file mode 100644 index fae2bcf..0000000 --- a/lib/wsprd/unpk.c.obsolete +++ /dev/null @@ -1,124 +0,0 @@ -/* the routine unpk() is not in wsprd_utils.c */ -#include -#include -#include -#include -#include -#include -#include - -#include "wsprd_utils.h" - -unsigned int nhash_( const void *key, size_t length, uint32_t initval); - -void unpk_(signed char message[], int *nhashtab, char call_loc_pow[]) -{ - int i,n1,n2,n3,ndbm,ihash,nadd,noprint,nh; - char callsign[13],grid[5],grid6[7],cdbm[3]; - static char hashtab[32768][13]; - FILE *fhash; - - if(*nhashtab==1) { - char line[80], hcall[12]; - if( (fhash=fopen("hashtable.txt","r+")) ) { - while (fgets(line, sizeof(line), fhash) != NULL) { - sscanf(line,"%d %s",&nh,hcall); - strcpy(*hashtab+nh*13,hcall); - } - } else { - fhash=fopen("hashtable.txt","w+"); - } - fclose(fhash); - return; - } - - if(*nhashtab==2) { - fhash=fopen("hashtable.txt","w"); - for (i=0; i<32768; i++) { - if( strncmp(hashtab[i],"\0",1) != 0 ) { - fprintf(fhash,"%5d %s\n",i,*hashtab+i*13); - } - } - fclose(fhash); - return; - } - - unpack50(message,&n1,&n2); - unpackcall(n1,callsign); - unpackgrid(n2, grid); - int ntype = (n2&127) - 64; - callsign[12]=0; - grid[4]=0; - -/* - 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 ) { - ndbm=ntype; - memset(call_loc_pow,0,sizeof(char)*23); - sprintf(cdbm,"%2d",ndbm); - strncat(call_loc_pow,callsign,strlen(callsign)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,grid,4); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,cdbm,2); - strncat(call_loc_pow,"\0",1); - ihash=nhash_(callsign,strlen(callsign),(uint32_t)146); - strcpy(*hashtab+ihash*13,callsign); - } else { - nadd=nu; - if( nu > 3 ) nadd=nu-3; - if( nu > 7 ) nadd=nu-7; - n3=n2/128+32768*(nadd-1); - unpackpfx(n3,callsign); - ndbm=ntype-nadd; - memset(call_loc_pow,0,sizeof(char)*23); - sprintf(cdbm,"%2d",ndbm); - strncat(call_loc_pow,callsign,strlen(callsign)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,cdbm,2); - strncat(call_loc_pow,"\0",1); - ihash=nhash_(callsign,strlen(callsign),(uint32_t)146); - strcpy(*hashtab+ihash*13,callsign); - noprint=0; - } - } else if ( ntype < 0 ) { - ndbm=-(ntype+1); - memset(grid6,0,sizeof(char)*7); - strncat(grid6,callsign+5,1); - strncat(grid6,callsign,5); - ihash=(n2-ntype-64)/128; - if( strncmp(hashtab[ihash],"\0",1) != 0 ) { - sprintf(callsign,"<%s>",hashtab[ihash]); - } else { - sprintf(callsign,"%5s","<...>"); - } - - memset(call_loc_pow,0,sizeof(char)*23); - sprintf(cdbm,"%2d",ndbm); - strncat(call_loc_pow,callsign,strlen(callsign)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,grid6,strlen(grid6)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,cdbm,2); - strncat(call_loc_pow,"\0",1); - - noprint=0; - -// I don't know what to do with these... They show up as "A000AA" grids. - if( ntype == -64 ) noprint=1; - } - // printf("\nUnpacked in C: %s\n",call_loc_pow); -} diff --git a/lib/wsprd/wsprd.c b/lib/wsprd/wsprd.c deleted file mode 100644 index 48ec1fa..0000000 --- a/lib/wsprd/wsprd.c +++ /dev/null @@ -1,1453 +0,0 @@ -/* - This file is part of program wsprd, a detector/demodulator/decoder - for the Weak Signal Propagation Reporter (WSPR) mode. - - File name: wsprd.c - - Copyright 2001-2018, 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-2018, 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 -#include -#include -#include -#include -#include -#include -#include - -#include "fano.h" -#include "jelinek.h" -#include "nhash.h" -#include "wsprd_utils.h" -#include "wsprsim_utils.h" - -#define max(x,y) ((x) > (y) ? (x) : (y)) -// Possible PATIENCE options: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT, -// FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE -#define PATIENCE FFTW_ESTIMATE -fftwf_plan PLAN1,PLAN2,PLAN3; - -unsigned char pr3[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}; - -unsigned long nr; - -int printdata=0; - -//*************************************************************************** -unsigned long readc2file(char *ptr_to_infile, float *idat, float *qdat, - double *freq, int *wspr_type) -{ - float *buffer; - double dfreq; - int i,ntrmin; - char *c2file[15]; - FILE* fp; - - buffer=calloc(2*65536,sizeof(float)); - - fp = fopen(ptr_to_infile,"rb"); - if (fp == NULL) { - fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile); - return 1; - } - unsigned long nread=fread(c2file,sizeof(char),14,fp); - nread=fread(&ntrmin,sizeof(int),1,fp); - nread=fread(&dfreq,sizeof(double),1,fp); - *freq=dfreq; - nread=fread(buffer,sizeof(float),2*45000,fp); - fclose(fp); - - *wspr_type=ntrmin; - - for(i=0; i<45000; i++) { - idat[i]=buffer[2*i]; - qdat[i]=-buffer[2*i+1]; - } - - if( nread == 2*45000 ) { - return nread/2; - } else { - return 1; - } - free(buffer); -} - -//*************************************************************************** -unsigned long readwavfile(char *ptr_to_infile, int ntrmin, float *idat, float *qdat ) -{ - size_t i, j, npoints; - int nfft1, nfft2, nh2, i0; - double df; - - nfft2=46080; //this is the number of downsampled points that will be returned - nh2=nfft2/2; - - if( ntrmin == 2 ) { - nfft1=nfft2*32; //need to downsample by a factor of 32 - df=12000.0/nfft1; - i0=1500.0/df+0.5; - npoints=114*12000; - } else if ( ntrmin == 15 ) { - nfft1=nfft2*8*32; - df=12000.0/nfft1; - i0=(1500.0+112.5)/df+0.5; - npoints=8*114*12000; - } else { - fprintf(stderr,"This should not happen\n"); - return 1; - } - - float *realin; - fftwf_complex *fftin, *fftout; - - FILE *fp; - short int *buf2; - buf2 = calloc(npoints,sizeof(short int)); - - fp = fopen(ptr_to_infile,"rb"); - if (fp == NULL) { - fprintf(stderr, "Cannot open data file '%s'\n", ptr_to_infile); - return 1; - } - nr=fread(buf2,2,22,fp); //Read and ignore header - nr=fread(buf2,2,npoints,fp); //Read raw data - fclose(fp); - - realin=(float*) fftwf_malloc(sizeof(float)*nfft1); - fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*(nfft1/2+1)); - PLAN1 = fftwf_plan_dft_r2c_1d(nfft1, realin, fftout, PATIENCE); - - for (i=0; i(size_t)nh2 ) j=j-nfft2; - fftin[i][0]=fftout[j][0]; - fftin[i][1]=fftout[j][1]; - } - - fftwf_free(fftout); - fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*nfft2); - PLAN2 = fftwf_plan_dft_1d(nfft2, fftin, fftout, FFTW_BACKWARD, PATIENCE); - fftwf_execute(PLAN2); - - for (i=0; i<(size_t)nfft2; i++) { - idat[i]=fftout[i][0]/1000.0; - qdat[i]=fftout[i][1]/1000.0; - } - - fftwf_free(fftin); - fftwf_free(fftout); - return nfft2; -} - -//*************************************************************************** -void sync_and_demodulate(float *id, float *qd, long np, - unsigned char *symbols, float *f1, int ifmin, int ifmax, float fstep, - int *shift1, int lagmin, int lagmax, int lagstep, - float *drift1, int symfac, float *sync, int mode) -{ - /*********************************************************************** - * mode = 0: no frequency or drift search. find best time lag. * - * 1: no time lag or drift search. find best frequency. * - * 2: no frequency or time lag search. calculate soft-decision * - * symbols using passed frequency and shift. * - ************************************************************************/ - - static float fplast=-10000.0; - static float dt=1.0/375.0, df=375.0/256.0; - static float pi=3.14159265358979323846; - float twopidt, df15=df*1.5, df05=df*0.5; - - int i, j, k, lag; - float i0[162],q0[162],i1[162],q1[162],i2[162],q2[162],i3[162],q3[162]; - float p0,p1,p2,p3,cmet,totp,syncmax,fac; - float c0[256],s0[256],c1[256],s1[256],c2[256],s2[256],c3[256],s3[256]; - float 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[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; i<162; i++) { - fp = f0 + (*drift1/2.0)*((float)i-81.0)/81.0; - if( i==0 || (fp != fplast) ) { // only calculate sin/cos if necessary - dphi0=twopidt*(fp-df15); - cdphi0=cos(dphi0); - sdphi0=sin(dphi0); - - dphi1=twopidt*(fp-df05); - cdphi1=cos(dphi1); - sdphi1=sin(dphi1); - - dphi2=twopidt*(fp+df05); - cdphi2=cos(dphi2); - sdphi2=sin(dphi2); - - dphi3=twopidt*(fp+df15); - cdphi3=cos(dphi3); - sdphi3=sin(dphi3); - - c0[0]=1; s0[0]=0; - c1[0]=1; s1[0]=0; - c2[0]=1; s2[0]=0; - c3[0]=1; s3[0]=0; - - for (j=1; j<256; j++) { - c0[j]=c0[j-1]*cdphi0 - s0[j-1]*sdphi0; - s0[j]=c0[j-1]*sdphi0 + s0[j-1]*cdphi0; - c1[j]=c1[j-1]*cdphi1 - s1[j-1]*sdphi1; - s1[j]=c1[j-1]*sdphi1 + s1[j-1]*cdphi1; - c2[j]=c2[j-1]*cdphi2 - s2[j-1]*sdphi2; - s2[j]=c2[j-1]*sdphi2 + s2[j-1]*cdphi2; - c3[j]=c3[j-1]*cdphi3 - s3[j-1]*sdphi3; - s3[j]=c3[j-1]*sdphi3 + s3[j-1]*cdphi3; - } - fplast = fp; - } - - i0[i]=0.0; q0[i]=0.0; - i1[i]=0.0; q1[i]=0.0; - i2[i]=0.0; q2[i]=0.0; - i3[i]=0.0; q3[i]=0.0; - - for (j=0; j<256; j++) { - k=lag+i*256+j; - if( (k>0) && (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<162; i++) { //Normalize the soft symbols - fsum=fsum+fsymb[i]/162.0; - f2sum=f2sum+fsymb[i]*fsymb[i]/162.0; - } - fac=sqrt(f2sum-fsum*fsum); - for (i=0; i<162; i++) { - fsymb[i]=symfac*fsymb[i]/fac; - if( fsymb[i] > 127) fsymb[i]=127.0; - if( fsymb[i] < -128 ) fsymb[i]=-128.0; - symbols[i]=fsymb[i] + 128; - } - return; - } - return; -} - -void noncoherent_sequence_detection(float *id, float *qd, long np, - unsigned char *symbols, float *f1, int *shift1, - float *drift1, int symfac, int *nblocksize) -{ - /************************************************************************ - * Noncoherent sequence detection for wspr. * - * Allowed block lengths are nblock=1,2,3,6, or 9 symbols. * - * Longer block lengths require longer channel coherence time. * - * The whole block is estimated at once. * - * nblock=1 corresponds to noncoherent detection of individual symbols * - * like the original wsprd symbol demodulator. * - ************************************************************************/ - static float fplast=-10000.0; - static float dt=1.0/375.0, df=375.0/256.0; - static float pi=3.14159265358979323846; - float twopidt, df15=df*1.5, df05=df*0.5; - - int i, j, k, lag, itone, ib, b, nblock, nseq, imask; - float xi[512],xq[512]; - float is[4][162],qs[4][162],cf[4][162],sf[4][162],cm,sm,cmp,smp; - float p[512],fac,xm1,xm0; - float c0[257],s0[257],c1[257],s1[257],c2[257],s2[257],c3[257],s3[257]; - float dphi0, cdphi0, sdphi0, dphi1, cdphi1, sdphi1, dphi2, cdphi2, sdphi2, - dphi3, cdphi3, sdphi3; - float f0, fp, fsum=0.0, f2sum=0.0, fsymb[162]; - - twopidt=2*pi*dt; - f0=*f1; - lag=*shift1; - nblock=*nblocksize; - nseq=1<0) && (k>(nblock-1-ib); - itone=pr3[i+ib]+2*b; - xi[j]=xi[j]+is[itone][i+ib]*cm + qs[itone][i+ib]*sm; - xq[j]=xq[j]+qs[itone][i+ib]*cm - is[itone][i+ib]*sm; - cmp=cf[itone][i+ib]*cm - sf[itone][i+ib]*sm; - smp=sf[itone][i+ib]*cm + cf[itone][i+ib]*sm; - cm=cmp; sm=smp; - } - p[j]=xi[j]*xi[j]+xq[j]*xq[j]; - p[j]=sqrt(p[j]); - } - for (ib=0; ib xm1) xm1=p[j]; - } - if((j & imask)==0) { - if(p[j]>xm0) xm0=p[j]; - } - } - fsymb[i+ib]=xm1-xm0; - } - } - for (i=0; i<162; i++) { //Normalize the soft symbols - fsum=fsum+fsymb[i]/162.0; - f2sum=f2sum+fsymb[i]*fsymb[i]/162.0; - } - fac=sqrt(f2sum-fsum*fsum); - for (i=0; i<162; i++) { - fsymb[i]=symfac*fsymb[i]/fac; - if( fsymb[i] > 127) fsymb[i]=127.0; - if( fsymb[i] < -128 ) fsymb[i]=-128.0; - symbols[i]=fsymb[i] + 128; - } - return; -} - -/*************************************************************************** - symbol-by-symbol signal subtraction - ****************************************************************************/ -void subtract_signal(float *id, float *qd, long np, - float f0, int shift0, float drift0, unsigned char* channel_symbols) -{ - float dt=1.0/375.0, df=375.0/256.0; - int i, j, k; - float pi=4.*atan(1.0),twopidt, fp; - - float i0,q0; - float c0[256],s0[256]; - float dphi, cdphi, sdphi; - - twopidt=2*pi*dt; - - for (i=0; i<162; i++) { - fp = f0 + ((float)drift0/2.0)*((float)i-81.0)/81.0; - - dphi=twopidt*(fp+((float)channel_symbols[i]-1.5)*df); - cdphi=cos(dphi); - sdphi=sin(dphi); - - c0[0]=1; s0[0]=0; - - for (j=1; j<256; j++) { - c0[j]=c0[j-1]*cdphi - s0[j-1]*sdphi; - s0[j]=c0[j-1]*sdphi + s0[j-1]*cdphi; - } - - i0=0.0; q0=0.0; - - for (j=0; j<256; j++) { - k=shift0+i*256+j; - if( (k>0) & (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) && (k path to writeable data files, default=\".\"\n"); - printf(" -B disable block demodulation - use single-symbol noncoherent demod\n"); - printf(" -c write .c2 file at the end of the first pass\n"); - printf(" -C maximum number of decoder cycles per bit, default 10000\n"); - printf(" -d deeper search. Slower, a few more decodes\n"); - printf(" -e x (x is transceiver dial frequency error in Hz)\n"); - printf(" -f x (x is transceiver dial frequency in MHz)\n"); - printf(" -H do not use (or update) the hash table\n"); - printf(" -J use the stack decoder instead of Fano decoder\n"); - printf(" -m decode wspr-15 .wav file\n"); - printf(" -q quick mode - doesn't dig deep for weak signals\n"); - printf(" -s single pass mode, no subtraction (same as original wsprd)\n"); - printf(" -v verbose mode (shows dupes)\n"); - printf(" -w wideband mode - decode signals within +/- 150 Hz of center\n"); - printf(" -z x (x is fano metric table bias, default is 0.45)\n"); -} - -//*************************************************************************** -int main(int argc, char *argv[]) -{ - char cr[] = "(C) 2018, Steven Franke - K9AN"; - (void)cr; - extern char *optarg; - extern int optind; - int i,j,k; - unsigned char *symbols, *decdata, *channel_symbols; - signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0}; - char *callsign, *call_loc_pow; - char *ptr_to_infile,*ptr_to_infile_suffix; - char *data_dir=NULL; - char wisdom_fname[200],all_fname[200],spots_fname[200]; - char timer_fname[200],hash_fname[200]; - char uttime[5],date[7]; - int c,delta,maxpts=65536,verbose=0,quickmode=0,more_candidates=0, stackdecoder=0; - int writenoise=0,usehashtable=1,wspr_type=2, ipass, nblocksize; - int writec2=0,maxdrift; - int shift1, lagmin, lagmax, lagstep, ifmin, ifmax, worth_a_try, not_decoded; - unsigned int nbits=81, stacksize=200000; - unsigned int npoints, metric, cycles, maxnp; - float df=375.0/256.0/2; - float freq0[200],snr0[200],drift0[200],sync0[200]; - int shift0[200]; - float dt=1.0/375.0, dt_print; - double dialfreq_cmdline=0.0, dialfreq, freq_print; - double dialfreq_error=0.0; - float fmin=-110, fmax=110; - float f1, fstep, sync1, drift1; - float psavg[512]; - float *idat, *qdat; - clock_t t0,t00; - float tfano=0.0,treadwav=0.0,tcandidates=0.0,tsync0=0.0; - float tsync1=0.0,tsync2=0.0,ttotal=0.0; - - struct result { char date[7]; char time[5]; float sync; float snr; - float dt; double freq; char message[23]; float drift; - unsigned int cycles; int jitter; int blocksize; unsigned int metric; }; - struct result decodes[50]; - - char *hashtab; - hashtab=calloc(32768*13,sizeof(char)); - int nh; - symbols=calloc(nbits*2,sizeof(char)); - decdata=calloc(11,sizeof(char)); - channel_symbols=calloc(nbits*2,sizeof(char)); - callsign=calloc(13,sizeof(char)); - call_loc_pow=calloc(23,sizeof(char)); - float allfreqs[100]; - char allcalls[100][13]; - for (i=0; i<100; i++) allfreqs[i]=0.0; - memset(allcalls,0,sizeof(char)*100*13); - - int uniques=0, noprint=0, ndecodes_pass=0; - - // Parameters used for performance-tuning: - unsigned int maxcycles=10000; //Decoder timeout limit - float minsync1=0.10; //First sync limit - float minsync2=0.12; //Second sync limit - int iifac=8; //Step size in final DT peakup - int symfac=50; //Soft-symbol normalizing factor - int block_demod=1; //Default is to use block demod on pass 2 - int subtraction=1; - int npasses=2; - - float minrms=52.0 * (symfac/64.0); //Final test for plausible decoding - delta=60; //Fano threshold step - float bias=0.45; //Fano metric bias (used for both Fano and stack algorithms) - - t00=clock(); - fftwf_complex *fftin, *fftout; -#include "./metric_tables.c" - - int mettab[2][256]; - - idat=calloc(maxpts,sizeof(float)); - qdat=calloc(maxpts,sizeof(float)); - - while ( (c = getopt(argc, argv, "a:BcC:de:f:HJmqstwvz:")) !=-1 ) { - switch (c) { - case 'a': - data_dir = optarg; - break; - case 'B': - block_demod=0; - break; - case 'c': - writec2=1; - break; - case 'C': - maxcycles=(unsigned int) strtoul(optarg,NULL,10); - break; - case 'd': - more_candidates=1; - break; - case 'e': - dialfreq_error = strtod(optarg,NULL); // units of Hz - // dialfreq_error = dial reading - actual, correct frequency - break; - case 'f': - dialfreq_cmdline = strtod(optarg,NULL); // units of MHz - break; - case 'H': - usehashtable = 0; - break; - case 'J': //Stack (Jelinek) decoder, Fano decoder is the default - stackdecoder = 1; - break; - case 'm': //15-minute wspr mode - wspr_type = 15; - break; - case 'q': //no shift jittering - quickmode = 1; - break; - case 's': //single pass mode - subtraction = 0; - npasses = 1; - break; - case 'v': - verbose = 1; - break; - case 'w': - fmin=-150.0; - fmax=150.0; - break; - case 'z': - bias=strtod(optarg,NULL); //fano metric bias (default is 0.45) - break; - case '?': - usage(); - return 1; - } - } - - if( stackdecoder ) { - stack=calloc(stacksize,sizeof(struct snode)); - } - - if( optind+1 > argc) { - usage(); - return 1; - } else { - ptr_to_infile=argv[optind]; - } - - // setup metric table - for(i=0; i<256; i++) { - mettab[0][i]=round( 10*(metric_tables[2][i]-bias) ); - mettab[1][i]=round( 10*(metric_tables[2][255-i]-bias) ); - } - - FILE *fp_fftwf_wisdom_file, *fall_wspr, *fwsprd, *fhash, *ftimer; - strcpy(wisdom_fname,"."); - strcpy(all_fname,"."); - strcpy(spots_fname,"."); - strcpy(timer_fname,"."); - strcpy(hash_fname,"."); - if(data_dir != NULL) { - strcpy(wisdom_fname,data_dir); - strcpy(all_fname,data_dir); - strcpy(spots_fname,data_dir); - strcpy(timer_fname,data_dir); - strcpy(hash_fname,data_dir); - } - strncat(wisdom_fname,"/wspr_wisdom.dat",20); - strncat(all_fname,"/ALL_WSPR.TXT",20); - strncat(spots_fname,"/wspr_spots.txt",20); - strncat(timer_fname,"/wspr_timer.out",20); - strncat(hash_fname,"/hashtable.txt",20); - if ((fp_fftwf_wisdom_file = fopen(wisdom_fname, "r"))) { //Open FFTW wisdom - fftwf_import_wisdom_from_file(fp_fftwf_wisdom_file); - fclose(fp_fftwf_wisdom_file); - } - - fall_wspr=fopen(all_fname,"a"); - fwsprd=fopen(spots_fname,"w"); - // FILE *fdiag; - // fdiag=fopen("wsprd_diag","a"); - - if((ftimer=fopen(timer_fname,"r"))) { - //Accumulate timing data - nr=fscanf(ftimer,"%f %f %f %f %f %f %f", - &treadwav,&tcandidates,&tsync0,&tsync1,&tsync2,&tfano,&ttotal); - fclose(ftimer); - } - ftimer=fopen(timer_fname,"w"); - - if( strstr(ptr_to_infile,".wav") ) { - ptr_to_infile_suffix=strstr(ptr_to_infile,".wav"); - - t0 = clock(); - npoints=readwavfile(ptr_to_infile, wspr_type, idat, qdat); - treadwav += (float)(clock()-t0)/CLOCKS_PER_SEC; - - if( npoints == 1 ) { - return 1; - } - dialfreq=dialfreq_cmdline - (dialfreq_error*1.0e-06); - } else if ( strstr(ptr_to_infile,".c2") !=0 ) { - ptr_to_infile_suffix=strstr(ptr_to_infile,".c2"); - npoints=readc2file(ptr_to_infile, idat, qdat, &dialfreq, &wspr_type); - if( npoints == 1 ) { - return 1; - } - dialfreq -= (dialfreq_error*1.0e-06); - } else { - printf("Error: Failed to open %s\n",ptr_to_infile); - printf("WSPR file must have suffix .wav or .c2\n"); - return 1; - } - - // Parse date and time from given filename - strncpy(date,ptr_to_infile_suffix-11,6); - strncpy(uttime,ptr_to_infile_suffix-4,4); - date[6]='\0'; - uttime[4]='\0'; - - // Do windowed ffts over 2 symbols, stepped by half symbols - int nffts=4*floor(npoints/512)-1; - fftin=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512); - fftout=(fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*512); - PLAN3 = fftwf_plan_dft_1d(512, fftin, fftout, FFTW_FORWARD, PATIENCE); - - float ps[512][nffts]; - float w[512]; - for(i=0; i<512; i++) { - w[i]=sin(0.006147931*i); - } - - if( usehashtable ) { - char line[80], hcall[12]; - if( (fhash=fopen(hash_fname,"r+")) ) { - while (fgets(line, sizeof(line), fhash) != NULL) { - sscanf(line,"%d %s",&nh,hcall); - strcpy(hashtab+nh*13,hcall); - } - } else { - fhash=fopen(hash_fname,"w+"); - } - fclose(fhash); - } - - //*************** main loop starts here ***************** - for (ipass=0; ipass511 ) - k=k-512; - ps[j][i]=fftout[k][0]*fftout[k][0]+fftout[k][1]*fftout[k][1]; - } - } - - // Compute average spectrum - for (i=0; i<512; i++) psavg[i]=0.0; - for (i=0; imin_snr) && (npk<200); - if ( candidate ) { - freq0[npk]=(j-205)*df; - snr0[npk]=10*log10(smspec[j])-snr_scaling_factor; - npk++; - } - } - } else { - for(j=1; j<410; j++) { - candidate = (smspec[j]>smspec[j-1]) && - (smspec[j]>smspec[j+1]) && - (npk<200); - if ( candidate ) { - freq0[npk]=(j-205)*df; - snr0[npk]=10*log10(smspec[j])-snr_scaling_factor; - npk++; - } - } - } - - // Compute corrected fmin, fmax, accounting for dial frequency error - fmin += dialfreq_error; // dialfreq_error is in units of Hz - fmax += dialfreq_error; - - // Don't waste time on signals outside of the range [fmin,fmax]. - i=0; - for( j=0; j= fmin && freq0[j] <= fmax ) { - freq0[i]=freq0[j]; - snr0[i]=snr0[j]; - i++; - } - } - npk=i; - - // bubble sort on snr, bringing freq along for the ride - int pass; - float tmp; - for (pass = 1; pass <= npk - 1; pass++) { - for (k = 0; k < npk - pass ; k++) { - if (snr0[k] < snr0[k+1]) { - tmp = snr0[k]; - snr0[k] = snr0[k+1]; - snr0[k+1] = tmp; - tmp = freq0[k]; - freq0[k] = freq0[k+1]; - freq0[k+1] = tmp; - } - } - } - - t0=clock(); - - /* 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,pow,p0,p1,p2,p3; - for(j=0; j smax ) { //Save coarse parameters - smax=sync1; - shift0[j]=128*(k0+1); - drift0[j]=idrift; - freq0[j]=(ifr-256)*df; - sync0[j]=sync1; - } - } - } - } - } - tcandidates += (float)(clock()-t0)/CLOCKS_PER_SEC; - - /* - 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. - */ - for (j=0; jminsync1 continue - fstep=0.0; ifmin=0; ifmax=0; - lagmin=shift1-128; - lagmax=shift1+128; - lagstep=64; - t0 = clock(); - sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1, - lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0); - tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC; - - fstep=0.25; ifmin=-2; ifmax=2; - t0 = clock(); - sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1, - lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1); - - if(ipass == 0) { - // 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, symbols, &f1, ifmin, ifmax, fstep, &shift1, - lagmin, lagmax, lagstep, &driftp, symfac, &syncp, 1); - - driftm=drift1-0.5; - sync_and_demodulate(idat, qdat, npoints, 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; - } - } - tsync1 += (float)(clock()-t0)/CLOCKS_PER_SEC; - - // fine-grid lag and freq search - if( sync1 > minsync1 ) { - - lagmin=shift1-32; lagmax=shift1+32; lagstep=16; - t0 = clock(); - sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1, - lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 0); - tsync0 += (float)(clock()-t0)/CLOCKS_PER_SEC; - - // fine search over frequency - fstep=0.05; ifmin=-2; ifmax=2; - t0 = clock(); - sync_and_demodulate(idat, qdat, npoints, symbols, &f1, ifmin, ifmax, fstep, &shift1, - lagmin, lagmax, lagstep, &drift1, symfac, &sync1, 1); - tsync1 += (float)(clock()-t0)/CLOCKS_PER_SEC; - - worth_a_try = 1; - } else { - worth_a_try = 0; - } - - int idt, ii, jittered_shift; - float y,sq,rms; - not_decoded=1; - int ib=1, blocksize; - while( ib <= nblocksize && not_decoded ) { - blocksize=ib; - idt=0; ii=0; - while ( worth_a_try && not_decoded && idt<=(128/iifac)) { - ii=(idt+1)/2; - if( idt%2 == 1 ) ii=-ii; - ii=iifac*ii; - jittered_shift=shift1+ii; - - // Use mode 2 to get soft-decision symbols - t0 = clock(); - noncoherent_sequence_detection(idat, qdat, npoints, symbols, &f1, - &jittered_shift, &drift1, symfac, &blocksize); - tsync2 += (float)(clock()-t0)/CLOCKS_PER_SEC; - - sq=0.0; - for(i=0; i<162; i++) { - y=(float)symbols[i] - 128.0; - sq += y*y; - } - rms=sqrt(sq/162.0); - - if((sync1 > minsync2) && (rms > minrms)) { - deinterleave(symbols); - t0 = clock(); - - if ( stackdecoder ) { - not_decoded = jelinek(&metric, &cycles, decdata, symbols, nbits, - stacksize, stack, mettab,maxcycles); - } else { - not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits, - mettab,delta,maxcycles); - } - - tfano += (float)(clock()-t0)/CLOCKS_PER_SEC; - - } - idt++; - if( quickmode ) break; - } - ib++; - } - - if( worth_a_try && !not_decoded ) { - ndecodes_pass++; - - for(i=0; i<11; i++) { - - if( decdata[i]>127 ) { - message[i]=decdata[i]-256; - } else { - message[i]=decdata[i]; - } - - } - - // 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). - noprint=unpk_(message,hashtab,call_loc_pow,callsign); - - // subtract even on last pass - if( subtraction && (ipass < npasses ) && !noprint ) { - if( get_wspr_channel_symbols(call_loc_pow, hashtab, channel_symbols) ) { - subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, channel_symbols); - } else { - break; - } - } - - // Remove dupes (same callsign and freq within 3 Hz) - int dupe=0; - for (i=0; i decodes[k+1].freq) { - temp = decodes[k]; - decodes[k]=decodes[k+1];; - decodes[k+1] = temp; - } - } - } - - for (i=0; i\n"); - - fftwf_free(fftin); - fftwf_free(fftout); - - if ((fp_fftwf_wisdom_file = fopen(wisdom_fname, "w"))) { - fftwf_export_wisdom_to_file(fp_fftwf_wisdom_file); - fclose(fp_fftwf_wisdom_file); - } - - ttotal += (float)(clock()-t00)/CLOCKS_PER_SEC; - - fprintf(ftimer,"%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n\n", - treadwav,tcandidates,tsync0,tsync1,tsync2,tfano,ttotal); - - fprintf(ftimer,"Code segment Seconds Frac\n"); - fprintf(ftimer,"-----------------------------------\n"); - fprintf(ftimer,"readwavfile %7.2f %7.2f\n",treadwav,treadwav/ttotal); - fprintf(ftimer,"Coarse DT f0 f1 %7.2f %7.2f\n",tcandidates, - tcandidates/ttotal); - fprintf(ftimer,"sync_and_demod(0) %7.2f %7.2f\n",tsync0,tsync0/ttotal); - fprintf(ftimer,"sync_and_demod(1) %7.2f %7.2f\n",tsync1,tsync1/ttotal); - fprintf(ftimer,"sync_and_demod(2) %7.2f %7.2f\n",tsync2,tsync2/ttotal); - fprintf(ftimer,"Stack/Fano decoder %7.2f %7.2f\n",tfano,tfano/ttotal); - fprintf(ftimer,"-----------------------------------\n"); - fprintf(ftimer,"Total %7.2f %7.2f\n",ttotal,1.0); - - fclose(fall_wspr); - fclose(fwsprd); - // fclose(fdiag); - fclose(ftimer); - fftwf_destroy_plan(PLAN1); - fftwf_destroy_plan(PLAN2); - fftwf_destroy_plan(PLAN3); - - if( usehashtable ) { - fhash=fopen(hash_fname,"w"); - for (i=0; i<32768; i++) { - if( strncmp(hashtab+i*13,"\0",1) != 0 ) { - fprintf(fhash,"%5d %s\n",i,hashtab+i*13); - } - } - fclose(fhash); - } - - free(hashtab); - free(symbols); - free(decdata); - free(channel_symbols); - free(callsign); - free(call_loc_pow); - free(idat); - free(qdat); - if( stackdecoder ) { - free(stack); - } - - if(writenoise == 999) return -1; //Silence compiler warning - return 0; -} diff --git a/lib/wsprd/wsprd_stats.txt b/lib/wsprd/wsprd_stats.txt deleted file mode 100644 index 0d925aa..0000000 --- a/lib/wsprd/wsprd_stats.txt +++ /dev/null @@ -1,24 +0,0 @@ - Linux Windows -Program Time Decodes Time Decodes -------------------------------------------------- -wsprd (Mar 2013) 2413 1451 2718 1451 - -k9an-wsprd 1800 2122 -k9an_wsprd -q 354 1939 - -wsprd 399 2190 356 2190 -wsprd -q 214 2034 192 2034 - -wsprd* 1240 2215 -wsprd# 1599 2220 - -------------------------------------------------- -* maxcycles=30000 -# maxcycles=20000, iifac=1 -------------------------------------------------- -Test data: 638 *.wav files (recorded by WSJT-X) -------------------------------------------------- -Linux machine: Core 2 Duo, E6750 CPU -Windows machine: 4-Core i5-2500 CPU -wsprd git commit: eecc274 -------------------------------------------------- diff --git a/lib/wsprd/wsprd_utils.c b/lib/wsprd/wsprd_utils.c deleted file mode 100644 index 65485d4..0000000 --- a/lib/wsprd/wsprd_utils.c +++ /dev/null @@ -1,339 +0,0 @@ -/* - This file is part of program wsprd, a detector/demodulator/decoder - for the Weak Signal Propagation Reporter (WSPR) mode. - - File name: wsprd_utils.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 "wsprd_utils.h" - -#ifndef int32_t -#define int32_t int -#endif - -void unpack50( signed char *dat, int32_t *n1, int32_t *n2 ) -{ - int32_t i,i4; - - i=dat[0]; - i4=i&255; - *n1=i4<<20; - - i=dat[1]; - i4=i&255; - *n1=*n1+(i4<<12); - - i=dat[2]; - i4=i&255; - *n1=*n1+(i4<<4); - - i=dat[3]; - i4=i&255; - *n1=*n1+((i4>>4)&15); - *n2=(i4&15)<<18; - - i=dat[4]; - i4=i&255; - *n2=*n2+(i4<<10); - - i=dat[5]; - i4=i&255; - *n2=*n2+(i4<<2); - - i=dat[6]; - i4=i&255; - *n2=*n2+((i4>>6)&3); -} - -int unpackcall( int32_t ncall, char *call ) -{ - char c[]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E', - 'F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T', - 'U','V','W','X','Y','Z',' '}; - int32_t n; - int i; - char tmp[7]; - - n=ncall; - strcpy(call,"......"); - if (n < 262177560 ) { - i=n%27+10; - tmp[5]=c[i]; - n=n/27; - i=n%27+10; - tmp[4]=c[i]; - n=n/27; - i=n%27+10; - tmp[3]=c[i]; - n=n/27; - i=n%10; - tmp[2]=c[i]; - n=n/10; - i=n%36; - tmp[1]=c[i]; - n=n/36; - i=n; - tmp[0]=c[i]; - tmp[6]='\0'; - // remove leading whitespace - for(i=0; i<5; i++) { - if( tmp[i] != c[36] ) - break; - } - sprintf(call,"%-6s",&tmp[i]); - // remove trailing whitespace - for(i=0; i<6; i++) { - if( call[i] == c[36] ) { - call[i]='\0'; - } - } - } else { - return 0; - } - return 1; -} - -int unpackgrid( int32_t ngrid, char *grid) -{ - char c[]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E', - 'F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T', - 'U','V','W','X','Y','Z',' '}; - int dlat, dlong; - - ngrid=ngrid>>7; - if( ngrid < 32400 ) { - dlat=(ngrid%180)-90; - dlong=(ngrid/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; - 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; - 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 ) { - // add a prefix of 1 to 3 characters - n=nprefix; - for (i=2; i>=0; i--) { - nc=n%37; - if( (nc >= 0) & (nc <= 9) ) { - pfx[i]=nc+48; - } - else if( (nc >= 10) & (nc <= 35) ) { - pfx[i]=nc+55; - } - 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+48; - strcpy(call,tmpcall); - strncat(call,"/",1); - strncat(call,pfx,1); - } - else if( (nc >= 10) & (nc <= 35) ) { - pfx[0]=nc+55; - strcpy(call,tmpcall); - strncat(call,"/",1); - strncat(call,pfx,1); - } - else if( (nc >= 36) & (nc <= 125) ) { - pfx[0]=(nc-26)/10+48; - pfx[1]=(nc-26)%10+48; - strcpy(call,tmpcall); - strncat(call,"/",1); - strncat(call,pfx,2); - } - else { - return 0; - } - } - return 1; -} - -void deinterleave(unsigned char *sym) -{ - unsigned char tmp[162]; - unsigned char p, i, j; - - p=0; - i=0; - while (p<162) { - j=((i * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32; - if (j < 162 ) { - tmp[p]=sym[j]; - p=p+1; - } - i=i+1; - } - for (i=0; i<162; i++) { - sym[i]=tmp[i]; - } -} - -// used by qsort -int doublecomp(const void* elem1, const void* elem2) -{ - if(*(const double*)elem1 < *(const double*)elem2) - return -1; - return *(const double*)elem1 > *(const double*)elem2; -} - -int floatcomp(const void* elem1, const void* elem2) -{ - if(*(const float*)elem1 < *(const float*)elem2) - return -1; - return *(const float*)elem1 > *(const float*)elem2; -} - -int unpk_(signed char *message, char *hashtab, char *call_loc_pow, char *callsign) -{ - int n1,n2,n3,ndbm,ihash,nadd,noprint=0; - char grid[5],grid6[7],cdbm[3]; - - unpack50(message,&n1,&n2); - if( !unpackcall(n1,callsign) ) return 1; - if( !unpackgrid(n2, grid) ) return 1; - int ntype = (n2&127) - 64; - callsign[12]=0; - grid[4]=0; - - /* - 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 ) { - ndbm=ntype; - memset(call_loc_pow,0,sizeof(char)*23); - sprintf(cdbm,"%2d",ndbm); - strncat(call_loc_pow,callsign,strlen(callsign)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,grid,4); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,cdbm,2); - strncat(call_loc_pow,"\0",1); - ihash=nhash(callsign,strlen(callsign),(uint32_t)146); - strcpy(hashtab+ihash*13,callsign); - } else { - nadd=nu; - if( nu > 3 ) nadd=nu-3; - if( nu > 7 ) nadd=nu-7; - n3=n2/128+32768*(nadd-1); - if( !unpackpfx(n3,callsign) ) return 1; - ndbm=ntype-nadd; - memset(call_loc_pow,0,sizeof(char)*23); - sprintf(cdbm,"%2d",ndbm); - strncat(call_loc_pow,callsign,strlen(callsign)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,cdbm,2); - strncat(call_loc_pow,"\0",1); - int nu=ndbm%10; - if( nu == 0 || nu == 3 || nu == 7 || nu == 10 ) { //make sure power is OK - ihash=nhash(callsign,strlen(callsign),(uint32_t)146); - strcpy(hashtab+ihash*13,callsign); - } else noprint=1; - } - } else if ( ntype < 0 ) { - ndbm=-(ntype+1); - memset(grid6,0,sizeof(char)*7); -// size_t len=strlen(callsign); - size_t len=6; - strncat(grid6,callsign+len-1,1); - strncat(grid6,callsign,len-1); - int nu=ndbm%10; - if ((nu != 0 && nu != 3 && nu != 7 && nu != 10) || - !isalpha(grid6[0]) || !isalpha(grid6[1]) || - !isdigit(grid6[2]) || !isdigit(grid6[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(grid6[4]) && isalpha(grid6[5]) ) ) { - noprint=1; - } - - ihash=(n2-ntype-64)/128; - if( strncmp(hashtab+ihash*13,"\0",1) != 0 ) { - sprintf(callsign,"<%s>",hashtab+ihash*13); - } else { - sprintf(callsign,"%5s","<...>"); - } - - memset(call_loc_pow,0,sizeof(char)*23); - sprintf(cdbm,"%2d",ndbm); - strncat(call_loc_pow,callsign,strlen(callsign)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,grid6,strlen(grid6)); - strncat(call_loc_pow," ",1); - strncat(call_loc_pow,cdbm,2); - strncat(call_loc_pow,"\0",1); - - - // I don't know what to do with these... They show up as "A000AA" grids. - if( ntype == -64 ) noprint=1; - } - return noprint; -} diff --git a/lib/wsprd/wsprd_utils.h b/lib/wsprd/wsprd_utils.h deleted file mode 100644 index d670836..0000000 --- a/lib/wsprd/wsprd_utils.h +++ /dev/null @@ -1,29 +0,0 @@ -#ifndef WSPRD_UTILS_H -#define WSPRD_UTILS_H - -#include -#include -#include -#include -#include -#include -#include -#include "nhash.h" - -void unpack50( signed char *dat, int32_t *n1, int32_t *n2 ); - -int unpackcall( int32_t ncall, char *call ); - -int unpackgrid( int32_t ngrid, char *grid); - -int unpackpfx( int32_t nprefix, char *call); - -void deinterleave(unsigned char *sym); - -// used by qsort -int doublecomp(const void* elem1, const void* elem2); -int floatcomp(const void* elem1, const void* elem2); - -int unpk_( signed char *message, char* hashtab, char *call_loc_pow, char *callsign); - -#endif diff --git a/lib/wsprd/wsprsim.c b/lib/wsprd/wsprsim.c deleted file mode 100644 index b3af998..0000000 --- a/lib/wsprd/wsprsim.c +++ /dev/null @@ -1,209 +0,0 @@ -/* - File name: wsprsim.c (first committed to wsjtx June 13, 2015) - */ -#include -#include -#include -#include - -#include "wsprsim_utils.h" -#include "wsprd_utils.h" -#include "fano.h" - -int printdata=0; - -void usage() { - printf("Usage: wsprsim [options] message\n"); - printf(" message format: \"K1ABC FN42 33\"\n"); - printf(" \"PJ4/K1ABC 33\"\n"); - printf(" \" FK52UD 33\"\n"); - printf("Options:\n"); - printf(" -c (print channel symbols)\n"); - printf(" -d (print packed data with zero tail - 11 bytes)\n"); - printf(" -f x (-100 Hz < f < 100 Hz)\n"); - printf(" -o filename (write a c2 file with this name)\n"); - printf(" -s x (x is snr of signal that is written to .c2 file)\n"); - printf("\n"); - printf(" e.g. ./wsprsim -cds -28 -o 150613_1920.c2 \"K1ABC FN42 33\"\n"); - printf(" then ./wsprd 150613_1920.c2\n"); -} - -int add_signal_vector(float f0, float t0, float amp, unsigned char* symbols - , double* isig, double* qsig) -{ - int i, j, ii, idelay; - double phi=0.0, twopidt, df, dt, dphi; - twopidt=8.0*atan(1.0)/375.0; - df=375.0/256.0; - dt=1/375.0; - idelay=t0/dt; - - for (i=0; i<162; i++) { - dphi=twopidt*(f0 + ( (double)symbols[i]-1.5)*df ); - for ( j=0; j<256; j++ ) { - ii=idelay+256*i+j; - isig[ii]=isig[ii]+amp*cos(phi); - qsig[ii]=qsig[ii]+amp*sin(phi); - phi=phi+dphi; - } - } - return 1; -} - -char* tobinary(int x) -{ - static char b[33]; - b[0] = '\0'; - - long unsigned int z; - for (z = 0x80000000; z > 0; z >>= 1) - { - strcat(b, ((x & z) == z) ? "1" : "0"); - } - - return b; -} - -double gaussrand() -{ - static double V1, V2, S; - static int phase = 0; - double X; - - if(phase == 0) { - do { - double U1 = (double)rand() / RAND_MAX; - double U2 = (double)rand() / RAND_MAX; - - V1 = 2 * U1 - 1; - V2 = 2 * U2 - 1; - S = V1 * V1 + V2 * V2; - } while(S >= 1 || S == 0); - - X = V1 * sqrt(-2 * log(S) / S); - } else - X = V2 * sqrt(-2 * log(S) / S); - - phase = 1 - phase; - - return X; -} - -unsigned long writec2file(char *c2filename, int trmin, double freq - , double *idat, double *qdat) -{ - int i; - float buffer[2*45000]; - memset(buffer,0,sizeof(float)*2*45000); - FILE *fp; - - fp = fopen(c2filename,"wb"); - if( fp == NULL ) { - fprintf(stderr, "Could not open c2 file '%s'\n", c2filename); - return 0; - } - unsigned long nwrite = fwrite(c2filename,sizeof(char),14,fp); - nwrite = fwrite(&trmin, sizeof(int), 1, fp); - nwrite = fwrite(&freq, sizeof(double), 1, fp); - - for(i=0; i<45000; i++) { - buffer[2*i]=idat[i]; - buffer[2*i+1]=-qdat[i]; - } - - nwrite = fwrite(buffer, sizeof(float), 2*45000, fp); - if( nwrite == 2*45000 ) { - return nwrite; - } else { - return 0; - } -} - - -//******************************************************************** -int main(int argc, char *argv[]) -{ - extern char *optarg; - extern int optind; - int i, c, printchannel=0, writec2=0; - float snr=50.0; - float f0=0.0, t0=1.0; - char *message, *c2filename, *hashtab; - c2filename=malloc(sizeof(char)*15); - hashtab=malloc(sizeof(char)*32768*13); - memset(hashtab,0,sizeof(char)*32768*13); - - // message length is 22 characters - message=malloc(sizeof(char)*23); - - strcpy(c2filename,"000000_0001.c2"); - - srand(getpid()); - - while ( (c = getopt(argc, argv, "cdf:o:s:")) !=-1 ) { - switch (c) { - case 'c': - printchannel=1; - break; - case 'd': - printdata=1; - break; - case 'f': - f0 = atof(optarg); - case 'o': - c2filename = optarg; - writec2=1; - break; - case 's': -// snr = (float)atoi(optarg); - snr = atof(optarg); - break; - } - } - - if( optind+1 > argc ) { - usage(); - return 0; - } else { - message=argv[optind]; - } - - unsigned char channel_symbols[162]; - get_wspr_channel_symbols(message, hashtab, channel_symbols); - - if( printchannel ) { - printf("Channel symbols:\n"); - for (i=0; i<162; i++) { - printf("%d ",channel_symbols[i]); - } - printf("\n"); - } - - // add noise, then signal - double isig[45000], qsig[45000]; - memset(isig,0,sizeof(double)*45000); - memset(qsig,0,sizeof(double)*45000); - - if( snr < 40 ) { - // snr in 375Hz is 8.2 dB higher than in 2500 Hz. - snr=snr+8.2; - snr=pow(10,snr/20.0)*pow(2,0.5); - - for (i = 0; i<45000; i++) { - isig[i]=isig[i]+gaussrand(); - qsig[i]=qsig[i]+gaussrand(); - } - } else { - snr=1.0; - } - - add_signal_vector(f0, t0, snr, channel_symbols, isig, qsig); - if( writec2) { - // write a .c2 file - double carrierfreq=10.1387; - int wsprtype=2; - printf("Writing %s\n",c2filename); - writec2file(c2filename, wsprtype, carrierfreq, isig, qsig); - } - return 1; -} diff --git a/lib/wsprd/wsprsim_utils.c b/lib/wsprd/wsprsim_utils.c deleted file mode 100644 index f65efba..0000000 --- a/lib/wsprd/wsprsim_utils.c +++ /dev/null @@ -1,312 +0,0 @@ -/* - Functions used by wsprsim - */ -#include "wsprsim_utils.h" -#include "wsprd_utils.h" -#include "nhash.h" -#include "fano.h" - -static char get_locator_character_code(char ch); -static char get_callsign_character_code(char ch); -static long unsigned int pack_grid4_power(char const *grid4, int power); -static long unsigned int pack_call(char const *callsign); -static void pack_prefix(char *callsign, int32_t *n, int32_t *m, int32_t *nadd ); -static void interleave(unsigned char *sym); - -char get_locator_character_code(char ch) { - if( ch >=48 && ch <=57 ) { //0-9 - return ch-48; - } - if( ch == 32 ) { //space - return 36; - } - if( ch >= 65 && ch <= 82 ) { //A-Z - return ch-65; - } - return -1; -} - -char get_callsign_character_code(char ch) { - if( ch >=48 && ch <=57 ) { //0-9 - return ch-48; - } - if( ch == 32 ) { //space - return 36; - } - if( ch >= 65 && ch <= 90 ) { //A-Z - return ch-55; - } - return -1; -} - -long unsigned int pack_grid4_power(char const *grid4, int power) { - long unsigned int m; - - m=(179-10*grid4[0]-grid4[2])*180+10*grid4[1]+grid4[3]; - m=m*128+power+64; - return m; -} - -long unsigned int pack_call(char const *callsign) { - unsigned int i; - long unsigned int n; - char call6[6]; - memset(call6,' ',sizeof(call6)); - // callsign is 6 characters in length. Exactly. - size_t call_len = strlen(callsign); - if( call_len > 6 ) { - return 0; - } - if( isdigit(callsign[2]) ) { - for (i=0; i= 48 && nc <= 57 ) { - *m=nc-48; - } else if ( nc >= 65 && nc <= 90 ) { - *m=nc-65+10; - } else { - *m=38; - } - *m=60000-32768+*m; - } else if( callsign[i1+3]==0 ) { - //two char suffix - for (i=0; i= 48 && nc <= 57 ) { - nc=nc-48; - } else if ( nc >= 65 && nc <= 90 ) { - nc=nc-65+10; - } else { - nc=36; - } - *m=37*(*m)+nc; - } - *nadd=0; - if( *m > 32768 ) { - *m=*m-32768; - *nadd=1; - } - } - free (call6); -} - -void interleave(unsigned char *sym) -{ - unsigned char tmp[162]; - unsigned char p, i, j; - - p=0; - i=0; - while (p<162) { - j=((i * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32; - if (j < 162 ) { - tmp[j]=sym[p]; - p=p+1; - } - i=i+1; - } - for (i=0; i<162; i++) { - sym[i]=tmp[i]; - } -} - -int get_wspr_channel_symbols(char* rawmessage, char* hashtab, unsigned char* symbols) { - int m=0, ntype=0; - long unsigned int n=0; - int i, j, ihash; - unsigned char pr3[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}; - int nu[10]={0,-1,1,0,-1,2,1,0,-1,1}; - char *callsign, *grid, *powstr; - char grid4[5], message[23]; - - memset(message,0,sizeof(char)*23); - i=0; - while ( rawmessage[i] != 0 && i<23 ) { - message[i]=rawmessage[i]; - i++; - } - - size_t i1=strcspn(message," "); - size_t i2=strcspn(message,"/"); - size_t i3=strcspn(message,"<"); - size_t i4=strcspn(message,">"); - size_t mlen=strlen(message); - - // Use the presence and/or absence of "<" and "/" to decide what - // type of message. No sanity checks! Beware! - - if( i1 > 3 && i1 < 7 && i2 == mlen && i3 == mlen ) { - // Type 1 message: K9AN EN50 33 - // xxnxxxx xxnn nn - callsign = strtok(message," "); - grid = strtok(NULL," "); - powstr = strtok(NULL," "); - int power = atoi(powstr); - n = pack_call(callsign); - - for (i=0; i<4; i++) { - grid4[i]=get_locator_character_code(*(grid+i)); - } - m = pack_grid4_power(grid4,power); - - } else if ( i3 == 0 && i4 < mlen ) { - // Type 3: EN50WC 33 - // FK52UD 37 - // send hash instead of callsign to make room for 6 char grid. - // if 4-digit locator is specified, 2 spaces are added to the end. - callsign=strtok(message,"<> "); - grid=strtok(NULL," "); - powstr=strtok(NULL," "); - int power = atoi(powstr); - if( power < 0 ) power=0; - if( power > 60 ) power=60; - power=power+nu[power%10]; - ntype=-(power+1); - ihash=nhash(callsign,strlen(callsign),(uint32_t)146); - m=128*ihash + ntype + 64; - - char grid6[7]; - memset(grid6,0,sizeof(char)*7); - j=strlen(grid); - for(i=0; istrlen(callsign) ) return 0; //guards against pathological case - powstr = strtok (NULL," "); - int power = atoi (powstr); - if( power < 0 ) power=0; - if( power > 60 ) power=60; - power=power+nu[power%10]; - int n1, ng, nadd; - pack_prefix(callsign, &n1, &ng, &nadd); - ntype=power + 1 + nadd; - m=128*ng+ntype+64; - n=n1; - } else { - return 0; - } - - // pack 50 bits + 31 (0) tail bits into 11 bytes - unsigned char it, data[11]; - memset(data,0,sizeof(char)*11); - it=0xFF & (n>>20); - data[0]=it; - it=0xFF & (n>>12); - data[1]=it; - it=0xFF & (n>>4); - data[2]=it; - it= ((n&(0x0F))<<4) + ((m>>18)&(0x0F)); - data[3]=it; - it=0xFF & (m>>10); - data[4]=it; - it=0xFF & (m>>2); - data[5]=it; - it=(m & 0x03)<<6 ; - data[6]=it; - data[7]=0; - data[8]=0; - data[9]=0; - data[10]=0; - - if( printdata ) { - printf("Data is :"); - for (i=0; i<11; i++) { - printf("%02X ",data[i]); - } - printf("\n"); - } - - // make sure that the 11-byte data vector is unpackable - // unpack it with the routine that the decoder will use and display - // the result. let the operator decide whether it worked. - - char *check_call_loc_pow, *check_callsign; - check_call_loc_pow=malloc(sizeof(char)*23); - check_callsign=malloc(sizeof(char)*13); - signed char check_data[11]; - memcpy(check_data,data,sizeof(char)*11); - - unpk_(check_data,hashtab,check_call_loc_pow,check_callsign); -// printf("Will decode as: %s\n",check_call_loc_pow); - - unsigned int nbytes=11; // The message with tail is packed into almost 11 bytes. - unsigned char channelbits[nbytes*8*2]; /* 162 rounded up */ - memset(channelbits,0,sizeof(char)*nbytes*8*2); - - encode(channelbits,data,nbytes); - - interleave(channelbits); - - for (i=0; i<162; i++) { - symbols[i]=2*channelbits[i]+pr3[i]; - } - free(check_call_loc_pow); - free(check_callsign); - return 1; -} diff --git a/lib/wsprd/wsprsim_utils.h b/lib/wsprd/wsprsim_utils.h deleted file mode 100644 index acffdee..0000000 --- a/lib/wsprd/wsprsim_utils.h +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef WSPRSIM_UTILS_H -#define WSPRSIM_UTILS_H - -#include -#include -#include -#include -#include -#include -#include - -extern int printdata; - -int get_wspr_channel_symbols(char* message, char* hashtab, unsigned char* symbols); - -#endif