chris/js8call
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removed some legacy wspr code

Browse Source
This commit is contained in:
Jordan Sherer 2019-01-28 22:59:50 -05:00
parent c4e02121f1
commit 896ca374bb
30 changed files with 0 additions and 5580 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
CMakeLists.txt
View File

@ -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}
)

1
lib/wsprcode/go.sh
View File

@ -1 +0,0 @@
gfortran -o wsprcode -fbounds-check wsprcode.f90 nhash.c

384
lib/wsprcode/nhash.c
View File

@ -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 <bob_jenkins@burtleburtle.net>
* 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 <stdio.h> /* defines printf for tests */
#include <time.h> /* defines time_t for timings in the test */
#ifdef Win32
#include "win_stdint.h" /* defines uint32_t etc */
#else
#include <stdint.h> /* defines uint32_t etc */
#endif
//#include <sys/param.h> /* attempt to define endianness */
//#ifdef linux
//# include <endian.h> /* 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<n; ++i) h = hashlittle( k[i], len[i], h);
By Bob Jenkins, 2006. bob_jenkins@burtleburtle.net. You may use this
code any way you wish, private, educational, or commercial. It's free.
Use for hash table lookup, or anything where one collision in 2^^32 is
acceptable. Do NOT use for cryptographic purposes.
-------------------------------------------------------------------------------
*/
//uint32_t hashlittle( const void *key, size_t length, uint32_t initval)
#ifdef STDCALL
uint32_t __stdcall NHASH( const void *key, size_t *length0, uint32_t *initval0)
#else
uint32_t nhash_( const void *key, int *length0, uint32_t *initval0)
#endif
{
uint32_t a,b,c; /* internal state */
size_t length;
uint32_t initval;
union { const void *ptr; size_t i; } u; /* needed for Mac Powerbook G4 */
length=*length0;
initval=*initval0;
/* Set up the internal state */
a = b = c = 0xdeadbeef + ((uint32_t)length) + initval;
u.ptr = key;
if (HASH_LITTLE_ENDIAN && ((u.i & 0x3) == 0)) {
const uint32_t *k = (const uint32_t *)key; /* read 32-bit chunks */
const uint8_t *k8;
k8=0; //Silence compiler warning
/*------ all but last block: aligned reads and affect 32 bits of (a,b,c) */
while (length > 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)

937
lib/wsprcode/wspr_old_subs.f90
View File

@ -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

158
lib/wsprcode/wsprcode.f90
View File

@ -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'

34
lib/wsprd/Makefile
View File

@ -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

30
lib/wsprd/Makefile.MinGW
View File

@ -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

40
lib/wsprd/Makefile.win32
View File

@ -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

54
lib/wsprd/README
View File

@ -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> 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.

132
lib/wsprd/WSPRcode.f90
View File

@ -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

239
lib/wsprd/fano.c
View File

@ -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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#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
}

39
lib/wsprd/fano.h
View File

@ -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

410
lib/wsprd/fftw3.h
View File

@ -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 <stdio.h>
#ifdef __cplusplus
extern "C"
{
#endif /* __cplusplus */
/* If <complex.h> 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 <stddef.h> /* 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 */

50
lib/wsprd/genmet.f90
View File

@ -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

28
lib/wsprd/gran.c
View File

@ -1,28 +0,0 @@
#include <stdlib.h>
#include <math.h>
/* 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;
}

163
lib/wsprd/jelinek.c
View File

@ -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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h> /* 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<nbits; i++){
metrics[i][0] = mettab[0][symbols[0]] + mettab[0][symbols[1]];
metrics[i][1] = mettab[0][symbols[0]] + mettab[1][symbols[1]];
metrics[i][2] = mettab[1][symbols[0]] + mettab[0][symbols[1]];
metrics[i][3] = mettab[1][symbols[0]] + mettab[1][symbols[1]];
symbols += 2;
}
// zero the stack
memset(stack,0,stacksize*sizeof(struct snode));
// initialize the loop variables
unsigned int lsym, ntail=31;
uint64_t encstate=0;
unsigned int nbuckets=1000;
unsigned int low_bucket=nbuckets-1; //will be set on first run-through
unsigned int high_bucket=0;
unsigned int *buckets, bucket;
buckets=malloc(nbuckets*sizeof(unsigned int));
memset(buckets,0,nbuckets*sizeof(unsigned int));
unsigned int ptr=1;
unsigned int stackptr=1; //pointer values of 0 are reserved (they mean that a bucket is empty)
unsigned int depth=0, nbits_minus_ntail=nbits-ntail;
unsigned int stacksize_minus_1=stacksize-1;
long int totmet0, totmet1, gamma=0;
unsigned int ncycles=maxcycles*nbits;
/********************* Start the stack decoder *****************/
for (i=1; i <= ncycles; i++) {
#ifdef DEBUG
printf("***stackptr=%ld, depth=%d, gamma=%d, encstate=%lx, bucket %d, low_bucket %d, high_bucket %d\n",
stackptr, depth, gamma, encstate, bucket, low_bucket, high_bucket);
#endif
// no need to store more than 7 bytes (56 bits) for encoder state because
// only 50 bits are not 0's.
if( depth < 56 ) {
encstate=encstate<<1;
ENCODE(lsym,encstate); // get channel symbols associated with the 0 branch
} else {
ENCODE(lsym,encstate<<(depth-55));
}
// lsym are the 0-branch channel symbols and 3^lsym are the 1-branch
// channel symbols (due to a special property of our generator polynomials)
totmet0 = gamma+metrics[depth][lsym]; // total metric for 0-branch daughter node
totmet1 = gamma+metrics[depth][3^lsym]; // total metric for 1-branch daughter node
depth++; //the depth of the daughter nodes
bucket=(totmet0>>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
}

26
lib/wsprd/jelinek.h
View File

@ -1,26 +0,0 @@
#ifndef JELINEK_H
#define JELINEK_H
#include <stdint.h>
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

138
lib/wsprd/metric_tables.c
View File

@ -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}};

76
lib/wsprd/mettab.c
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
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 };

11
lib/wsprd/t2.f90
View File

@ -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

41
lib/wsprd/tab.c
View File

@ -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,
};

61
lib/wsprd/test_wspr.f90
View File

@ -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

124
lib/wsprd/unpk.c.obsolete
View File

@ -1,124 +0,0 @@
/* the routine unpk() is not in wsprd_utils.c */
#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdint.h>
#include <time.h>
#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);
}

1453
lib/wsprd/wsprd.c
View File

File diff suppressed because it is too large Load Diff

24
lib/wsprd/wsprd_stats.txt
View File

@ -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
-------------------------------------------------

339
lib/wsprd/wsprd_utils.c
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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: <PA0SKT/2> 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;
}

29
lib/wsprd/wsprd_utils.h
View File

@ -1,29 +0,0 @@
#ifndef WSPRD_UTILS_H
#define WSPRD_UTILS_H
#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <stdint.h>
#include <time.h>
#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

209
lib/wsprd/wsprsim.c
View File

@ -1,209 +0,0 @@
/*
File name: wsprsim.c (first committed to wsjtx June 13, 2015)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#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(" \"<PJ4/K1ABC> 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;
}

312
lib/wsprd/wsprsim_utils.c
View File

@ -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<call_len; i++) {
call6[i]=callsign[i];
}
} else if( isdigit(callsign[1]) ) {
for (i=1; i<call_len+1; i++) {
call6[i]=callsign[i-1];
}
}
for (i=0; i<6; i++) {
call6[i]=get_callsign_character_code(call6[i]);
}
n = call6[0];
n = n*36+call6[1];
n = n*10+call6[2];
n = n*27+call6[3]-10;
n = n*27+call6[4]-10;
n = n*27+call6[5]-10;
return n;
}
void pack_prefix(char *callsign, int32_t *n, int32_t *m, int32_t *nadd ) {
size_t i;
char * call6 = calloc(7,sizeof (char));
size_t i1=strcspn(callsign,"/");
if( callsign[i1+2] == 0 ) {
//single char suffix
for (i=0; i<i1; i++) {
call6[i]=callsign[i];
}
call6[i] = '\0';
*n=pack_call(call6);
*nadd=1;
int nc = callsign[i1+1];
if( nc >= 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<i1; i++) {
call6[i]=callsign[i];
}
*n=pack_call(call6);
*nadd=1;
*m=10*(callsign[i1+1]-48)+(callsign[i1+2]-48);
*m=60000 + 26 + *m;
} else {
char const * pfx = strtok (callsign,"/");
char const * call = strtok(NULL," ");
*n = pack_call (call);
size_t plen=strlen (pfx);
if( plen ==1 ) {
*m=36;
*m=37*(*m)+36;
} else if( plen == 2 ) {
*m=36;
} else {
*m=0;
}
for (i=0; i<plen; i++) {
int nc = callsign[i];
if( nc >= 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: <K1ABC> EN50WC 33
// <PJ4/K1ABC> 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; i<j-1; i++) {
grid6[i]=grid[i+1];
}
grid6[5]=grid[0];
n = pack_call(grid6);
} else if ( i2 < mlen ) { // just looks for a right slash
// Type 2: PJ4/K1ABC 37
callsign = strtok (message," ");
if( i2==0 || i2>strlen(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;
}

16
lib/wsprd/wsprsim_utils.h
View File

@ -1,16 +0,0 @@
#ifndef WSPRSIM_UTILS_H
#define WSPRSIM_UTILS_H
#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <stdint.h>
#include <time.h>
extern int printdata;
int get_wspr_channel_symbols(char* message, char* hashtab, unsigned char* symbols);
#endif