diff --git a/CMakeLists.txt b/CMakeLists.txt index ed31da3d4..13f55763f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -353,6 +353,10 @@ set (wsjt_FSRCS lib/wavhdr.f90 lib/qra/q65/q65_encoding_modules.f90 lib/ft8/ft8_a7.f90 + lib/superfox/sfox_mod.f90 + lib/superfox/julian.f90 + lib/superfox/popen_module.f90 + lib/superfox/qpc/qpc_mod.f90 # remaining non-module sources lib/addit.f90 @@ -597,8 +601,18 @@ set (wsjt_FSRCS lib/fst4/get_crc24.f90 lib/fst4/fst4_baseline.f90 lib/77bit/hash22calc.f90 + + lib/superfox/qpc_decode2.f90 + lib/superfox/qpc_likelihoods2.f90 + lib/superfox/qpc_snr.f90 + lib/superfox/qpc_sync.f90 + lib/superfox/sfox_ana.f90 + lib/superfox/sfox_demod.f90 + lib/superfox/sfox_remove_ft8.f90 + lib/superfox/sfox_unpack.f90 lib/superfox/sfox_wave.f90 lib/superfox/sfox_wave_gfsk.f90 + lib/superfox/twkfreq2.f90 ) # temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit @@ -641,6 +655,15 @@ set (wsjt_CSRCS lib/wrapkarn.c ${ldpc_CSRCS} ${qra_CSRCS} + + lib/superfox/qpc/dbgprintf.c + lib/superfox/qpc/nhash2.c + lib/superfox/qpc/np_qpc.c + lib/superfox/qpc/np_rnd.c + lib/superfox/qpc/qpc_fwht.c + lib/superfox/qpc/qpc_n127k50q128.c + lib/superfox/qpc/qpc_subs.c + ) set (wsjt_qt_UISRCS @@ -1226,6 +1249,9 @@ target_link_libraries (echosim wsjt_fort wsjt_cxx) add_executable (ft8sim lib/ft8/ft8sim.f90) target_link_libraries (ft8sim wsjt_fort wsjt_cxx) +add_executable (sfrx lib/superfox/sfrx.f90) +target_link_libraries (sfrx wsjt_fort wsjt_cxx) + add_executable (msk144sim lib/msk144sim.f90) target_link_libraries (msk144sim wsjt_fort wsjt_cxx) @@ -1613,7 +1639,7 @@ install (TARGETS udp_daemon message_aggregator wsjtx_app_version BUNDLE DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT runtime ) -install (TARGETS jt9 wsprd fmtave fcal fmeasure +install (TARGETS jt9 wsprd fmtave fcal fmeasure sfrx RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT runtime BUNDLE DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT runtime ) @@ -1686,7 +1712,7 @@ if (WIN32) endif () install (FILES - ${sfox_dir}/sfrx.exe +# ${sfox_dir}/sfrx.exe ${sfox_dir}/sftx.exe ${sfox_dir}/foxchk.exe DESTINATION ${CMAKE_INSTALL_BINDIR} @@ -1696,7 +1722,7 @@ endif (WIN32) if (UNIX) install (FILES - lib/superfox/linux/sfrx +# lib/superfox/linux/sfrx lib/superfox/linux/sftx lib/superfox/linux/foxchk DESTINATION ${CMAKE_INSTALL_BINDIR} @@ -1709,7 +1735,7 @@ endif (UNIX) if (APPLE) install (FILES - lib/superfox/mac/sfrx +# lib/superfox/mac/sfrx lib/superfox/mac/sftx lib/superfox/mac/foxchk DESTINATION ${CMAKE_INSTALL_BINDIR} @@ -1729,7 +1755,7 @@ execute_process(COMMAND if(${CMAKE_DEB_HOST_ARCH} MATCHES "arm64") install (FILES - lib/superfox/arm/sfrx +# lib/superfox/arm/sfrx lib/superfox/arm/sftx lib/superfox/arm/foxchk DESTINATION ${CMAKE_INSTALL_BINDIR} @@ -1742,7 +1768,7 @@ endif() if(${CMAKE_DEB_HOST_ARCH} MATCHES "armhf") install (FILES - lib/superfox/arm32/sfrx +# lib/superfox/arm32/sfrx lib/superfox/arm32/sftx lib/superfox/arm32/foxchk DESTINATION ${CMAKE_INSTALL_BINDIR} diff --git a/lib/ft8/ft8_downsample.f90 b/lib/ft8/ft8_downsample.f90 index f35522583..a57b5324a 100644 --- a/lib/ft8/ft8_downsample.f90 +++ b/lib/ft8/ft8_downsample.f90 @@ -8,10 +8,10 @@ subroutine ft8_downsample(dd,newdat,f0,c1) logical newdat,first complex c1(0:NFFT2-1) complex cx(0:NFFT1/2) - real dd(NMAX),x(NFFT1),taper(0:100) - equivalence (x,cx) + real dd(NMAX),x(NFFT1+2),taper(0:100) data first/.true./ - save cx,first,taper + save x,cx,first,taper + equivalence (x,cx) if(first) then pi=4.0*atan(1.0) @@ -23,11 +23,10 @@ subroutine ft8_downsample(dd,newdat,f0,c1) if(newdat) then ! Data in dd have changed, recompute the long FFT x(1:NMAX)=dd - x(NMAX+1:NFFT1)=0. !Zero-pad the x array + x(NMAX+1:NFFT1+2)=0. !Zero-pad the x array call four2a(cx,NFFT1,1,-1,0) !r2c FFT to freq domain newdat=.false. endif - df=12000.0/NFFT1 baud=12000.0/NSPS i0=nint(f0/df) diff --git a/lib/superfox/ft8_params.f90 b/lib/superfox/ft8_params.f90 new file mode 100644 index 000000000..139ff964e --- /dev/null +++ b/lib/superfox/ft8_params.f90 @@ -0,0 +1,12 @@ +! LDPC (174,91) code +parameter (KK=91) !Information bits (77 + CRC14) +parameter (ND=58) !Data symbols +parameter (NS=21) !Sync symbols (3 @ Costas 7x7) +parameter (NN=NS+ND) !Total channel symbols (79) +parameter (NSPS=1920) !Samples per symbol at 12000 S/s +parameter (NZ=NSPS*NN) !Samples in full 15 s waveform (151,680) +parameter (NMAX=15*12000) !Samples in iwave (180,000) +parameter (NFFT1=2*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra +parameter (NSTEP=NSPS/4) !Rough time-sync step size +parameter (NHSYM=NMAX/NSTEP-3) !Number of symbol spectra (1/4-sym steps) +parameter (NDOWN=60) !Downsample factor diff --git a/lib/superfox/gtag.f90 b/lib/superfox/gtag.f90 new file mode 100644 index 000000000..9cf93ab37 --- /dev/null +++ b/lib/superfox/gtag.f90 @@ -0,0 +1,3 @@ + integer ntag + data ntag/z"1234567"/ + !----------------------------------------------------------------------- diff --git a/lib/superfox/julian.f90 b/lib/superfox/julian.f90 new file mode 100644 index 000000000..002167fca --- /dev/null +++ b/lib/superfox/julian.f90 @@ -0,0 +1,39 @@ +module julian + +contains + + integer*8 function itime8() + + implicit integer (a-z) +! 1 2 3 4 5 6 7 + integer it(8) !y m d nmin h m s + integer*8 secday + data secday/86400/ + + call date_and_time(values=it) + iyr=it(1) + imo=it(2) + iday=it(3) + days=JD(iyr,imo,iday) - 2440588 !Days since epoch Jan 1, 1970 + ih=it(5) + im=it(6)-it(4) !it(4) corrects for time zone + is=it(7) + nsec=3600*ih + 60*im + is + itime8=days*secday + nsec + + return + end function itime8 + + integer function JD(I,J,K) + +! Return Julian Date for I=year, J=month, K=day +! Reference: Fliegel and Van Flandern, Comm. ACM 11, 10, 1968 + + JD = K - 32075 + 1461*(I + 4800 + (J - 14)/12)/4 & + + 367*(J - 2 - (J - 14)/12*12)/12 - 3 & + *((I + 4900 + (J - 14)/12)/100)/4 + + return + end function JD + +end module julian diff --git a/lib/superfox/popen_module.f90 b/lib/superfox/popen_module.f90 new file mode 100644 index 000000000..9aea9ff55 --- /dev/null +++ b/lib/superfox/popen_module.f90 @@ -0,0 +1,103 @@ +!******************************************************************************* +!> +! A simple module for `popen`. + + module popen_module + + use,intrinsic :: iso_c_binding + + implicit none + + private + + ! interfaces to C functions + interface + + function popen(command, mode) bind(C,name='popen') + !! initiate pipe streams to or from a process + import :: c_char, c_ptr + implicit none + character(kind=c_char),dimension(*) :: command + character(kind=c_char),dimension(*) :: mode + type(c_ptr) :: popen + end function popen + + function fgets(s, siz, stream) bind(C,name='fgets') + !! get a string from a stream + import :: c_char, c_ptr, c_int + implicit none + type (c_ptr) :: fgets + character(kind=c_char),dimension(*) :: s + integer(kind=c_int),value :: siz + type(c_ptr),value :: stream + end function fgets + + function pclose(stream) bind(C,name='pclose') + !! close a pipe stream to or from a process + import :: c_ptr, c_int + implicit none + integer(c_int) :: pclose + type(c_ptr),value :: stream + end function pclose + + end interface + + public :: c2f_string, get_command_as_string + +contains + + !******************************************************************************* + !> + ! Convert a C string to a Fortran string. + + function c2f_string(c) result(f) + + character(len=*),intent(in) :: c !! C string + character(len=:),allocatable :: f !! Fortran string + + integer :: i + + i = index(c,c_null_char) + + if (i<=0) then + f = c + else if (i==1) then + f = '' + else if (i>1) then + f = c(1:i-1) + end if + + end function c2f_string + + !******************************************************************************* + !> + ! Return the result of the command as a string. + + function get_command_as_string(command) result(str) + + character(len=*),intent(in) :: command !! the command to run + character(len=:),allocatable :: str !! the result of that command + + integer,parameter :: buffer_length = 1000 !! read stream in chunks of this size + + type(c_ptr) :: h !! for `popen` + integer(c_int) :: istat !! `pclose` status + character(kind=c_char,len=buffer_length) :: line !! buffer to read from `fgets` + + str = '' + h = c_null_ptr + h = popen(command//c_null_char,'r'//c_null_char) + + if (c_associated(h)) then + do while (c_associated(fgets(line,buffer_length,h))) + str = str//c2f_string(line) + end do + istat = pclose(h) + end if + + end function get_command_as_string + +!******************************************************************************* + end module popen_module +!******************************************************************************* + diff --git a/lib/superfox/qpc/dbgprintf.c b/lib/superfox/qpc/dbgprintf.c new file mode 100644 index 000000000..4a86e5556 --- /dev/null +++ b/lib/superfox/qpc/dbgprintf.c @@ -0,0 +1,67 @@ +// ------------------------------------------------------------------------------ +// dbgprintf.c +// Functions for printing debug information +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . +// ------------------------------------------------------------------------------ + +#include "dbgprintf.h" + +// print functions for debug purposes +void dbgprintf_vector_uchar(const char* name, const unsigned char* v, int vsize) +{ + int k; + + printf("%s=", name); + for (k = 0; k < vsize; k++) { + if ((k & 0x0F) == 0) + printf("\n"); + printf("%02X ", v[k]); + } + printf("\n"); +} + +void dbgprintf_vector_float(const char* name, const float* v, int vsize) +{ + int k; + + printf("%s=", name); + for (k = 0; k < vsize; k++) { + if ((k & 0x07) == 0) + printf("\n"); + printf("%10.3f ", v[k]); + } + printf("\n"); +} + +void dbgprintf_rows_float(const char* name, const float* v, int vsize, int nrows) +{ + int k; + int j; + + printf("%s=\n", name); + for (j = 0; j < nrows; j++) { + printf("r%d:", j); + for (k = 0; k < vsize; k++) { + if ((k & 0x07) == 0) + printf("\n"); + printf("%7.3f ", v[k]); + } + printf("\n"); + v += vsize; + } +} diff --git a/lib/superfox/qpc/dbgprintf.h b/lib/superfox/qpc/dbgprintf.h new file mode 100644 index 000000000..dc906f097 --- /dev/null +++ b/lib/superfox/qpc/dbgprintf.h @@ -0,0 +1,64 @@ +// ------------------------------------------------------------------------------ +// dbgprintf.h +// Functions for printing debug information +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . +// ------------------------------------------------------------------------------ + +#ifndef _dbgprintf_h_ +#define _dbgprintf_h_ + +// Define DBGPRINTF_ANYWAY +// in order to print debug information also +// when _DEBUG is not defined + +#define DBGPRINTF_ANYWAY + +#ifdef _DEBUG +#define DBG_PRINTF +#else +#ifdef DBGPRINTF_ANYWAY +#define DBG_PRINTF +#endif +#endif + +#include + +// print functions for debug purposes +#ifdef DBG_PRINTF + +#ifdef __cplusplus +extern "C" { +#endif + + void dbgprintf_vector_uchar(const char* name, const unsigned char* v, int vsize); + void dbgprintf_vector_float(const char* name, const float* v, int vsize); + void dbgprintf_rows_float(const char* name, const float* v, int vsize, int nrows); + +#ifdef __cplusplus +} +#endif + +#else + +#define dbgprintf_vector_uchar(a,b) +#define dbgprintf_vector_float(a,b) +#define dbgprintf_rows_float(a, b, c, d) + +#endif + +#endif // _dbgprintf_h_ diff --git a/lib/superfox/qpc/nhash2.c b/lib/superfox/qpc/nhash2.c new file mode 100644 index 000000000..78817c7d9 --- /dev/null +++ b/lib/superfox/qpc/nhash2.c @@ -0,0 +1,383 @@ +/* + + File name: nhash2.c + + *------------------------------------------------------------------------------ + * + * This file is part of the WSPR application, Weak Signal Propogation Reporter + * + * File Name: nhash.c + * Description: Functions to produce 32-bit hashes for hash table lookup + * + * Copyright (C) 2008-2014 Joseph Taylor, K1JT + * License: GNU GPL v3+ + * + * This program is free software; you can redistribute it and/or modify it under + * the terms of the GNU General Public License as published by the Free Software + * Foundation; either version 3 of the License, or (at your option) any later + * version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * this program; if not, write to the Free Software Foundation, Inc., 51 Franklin + * Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * Files: lookup3.c + * Copyright: Copyright (C) 2006 Bob Jenkins + * License: public-domain + * You may use this code any way you wish, private, educational, or commercial. + * It's free. + * + *------------------------------------------------------------------------------- +*/ + +/* +These are functions for producing 32-bit hashes for hash table lookup. +hashword(), hashlittle(), hashlittle2(), hashbig(), mix(), and final() +are externally useful functions. Routines to test the hash are included +if SELF_TEST is defined. You can use this free for any purpose. It's in +the public domain. It has no warranty. + +You probably want to use hashlittle(). hashlittle() and hashbig() +hash byte arrays. hashlittle() is is faster than hashbig() on +little-endian machines. Intel and AMD are little-endian machines. +On second thought, you probably want hashlittle2(), which is identical to +hashlittle() except it returns two 32-bit hashes for the price of one. +You could implement hashbig2() if you wanted but I haven't bothered here. + +If you want to find a hash of, say, exactly 7 integers, do + a = i1; b = i2; c = i3; + mix(a,b,c); + a += i4; b += i5; c += i6; + mix(a,b,c); + a += i7; + final(a,b,c); +then use c as the hash value. If you have a variable length array of +4-byte integers to hash, use hashword(). If you have a byte array (like +a character string), use hashlittle(). If you have several byte arrays, or +a mix of things, see the comments above hashlittle(). + +Why is this so big? I read 12 bytes at a time into 3 4-byte integers, +then mix those integers. This is fast (you can do a lot more thorough +mixing with 12*3 instructions on 3 integers than you can with 3 instructions +on 1 byte), but shoehorning those bytes into integers efficiently is messy. +*/ + +#define SELF_TEST 1 + +#include /* defines printf for tests */ +#include /* defines time_t for timings in the test */ +#include "nhash2.h" +//#include /* attempt to define endianness */ +//#ifdef linux +//# include /* attempt to define endianness */ +//#endif + +#define HASH_LITTLE_ENDIAN 1 + +#define hashsize(n) ((uint32_t)1<<(n)) +#define hashmask(n) (hashsize(n)-1) +#define rot(x,k) (((x)<<(k)) | ((x)>>(32-(k)))) + +/* +------------------------------------------------------------------------------- +mix -- mix 3 32-bit values reversibly. + +This is reversible, so any information in (a,b,c) before mix() is +still in (a,b,c) after mix(). + +If four pairs of (a,b,c) inputs are run through mix(), or through +mix() in reverse, there are at least 32 bits of the output that +are sometimes the same for one pair and different for another pair. +This was tested for: +* pairs that differed by one bit, by two bits, in any combination + of top bits of (a,b,c), or in any combination of bottom bits of + (a,b,c). +* "differ" is defined as +, -, ^, or ~^. For + and -, I transformed + the output delta to a Gray code (a^(a>>1)) so a string of 1's (as + is commonly produced by subtraction) look like a single 1-bit + difference. +* the base values were pseudorandom, all zero but one bit set, or + all zero plus a counter that starts at zero. + +Some k values for my "a-=c; a^=rot(c,k); c+=b;" arrangement that +satisfy this are + 4 6 8 16 19 4 + 9 15 3 18 27 15 + 14 9 3 7 17 3 +Well, "9 15 3 18 27 15" didn't quite get 32 bits diffing +for "differ" defined as + with a one-bit base and a two-bit delta. I +used http://burtleburtle.net/bob/hash/avalanche.html to choose +the operations, constants, and arrangements of the variables. + +This does not achieve avalanche. There are input bits of (a,b,c) +that fail to affect some output bits of (a,b,c), especially of a. The +most thoroughly mixed value is c, but it doesn't really even achieve +avalanche in c. + +This allows some parallelism. Read-after-writes are good at doubling +the number of bits affected, so the goal of mixing pulls in the opposite +direction as the goal of parallelism. I did what I could. Rotates +seem to cost as much as shifts on every machine I could lay my hands +on, and rotates are much kinder to the top and bottom bits, so I used +rotates. +------------------------------------------------------------------------------- +*/ +#define mix(a,b,c) \ +{ \ + a -= c; a ^= rot(c, 4); c += b; \ + b -= a; b ^= rot(a, 6); a += c; \ + c -= b; c ^= rot(b, 8); b += a; \ + a -= c; a ^= rot(c,16); c += b; \ + b -= a; b ^= rot(a,19); a += c; \ + c -= b; c ^= rot(b, 4); b += a; \ +} + +/* +------------------------------------------------------------------------------- +final -- final mixing of 3 32-bit values (a,b,c) into c + +Pairs of (a,b,c) values differing in only a few bits will usually +produce values of c that look totally different. This was tested for +* pairs that differed by one bit, by two bits, in any combination + of top bits of (a,b,c), or in any combination of bottom bits of + (a,b,c). +* "differ" is defined as +, -, ^, or ~^. For + and -, I transformed + the output delta to a Gray code (a^(a>>1)) so a string of 1's (as + is commonly produced by subtraction) look like a single 1-bit + difference. +* the base values were pseudorandom, all zero but one bit set, or + all zero plus a counter that starts at zero. + +These constants passed: + 14 11 25 16 4 14 24 + 12 14 25 16 4 14 24 +and these came close: + 4 8 15 26 3 22 24 + 10 8 15 26 3 22 24 + 11 8 15 26 3 22 24 +------------------------------------------------------------------------------- +*/ +#define final(a,b,c) \ +{ \ + c ^= b; c -= rot(b,14); \ + a ^= c; a -= rot(c,11); \ + b ^= a; b -= rot(a,25); \ + c ^= b; c -= rot(b,16); \ + a ^= c; a -= rot(c,4); \ + b ^= a; b -= rot(a,14); \ + c ^= b; c -= rot(b,24); \ +} + +/* +------------------------------------------------------------------------------- +hashlittle() -- hash a variable-length key into a 32-bit value + k : the key (the unaligned variable-length array of bytes) + length : the length of the key, counting by bytes + initval : can be any 4-byte value +Returns a 32-bit value. Every bit of the key affects every bit of +the return value. Two keys differing by one or two bits will have +totally different hash values. + +The best hash table sizes are powers of 2. There is no need to do +mod a prime (mod is sooo slow!). If you need less than 32 bits, +use a bitmask. For example, if you need only 10 bits, do + h = (h & hashmask(10)); +In which case, the hash table should have hashsize(10) elements. + +If you are hashing n strings (uint8_t **)k, do it like this: + for (i=0, h=0; i 12) + { + a += k[0]; + b += k[1]; + c += k[2]; + mix(a,b,c); + length -= 12; + k += 3; + } + + // printf("AAA %x %d %x %d\n",c,c,c&32767,c&32767); + + /*----------------------------- 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; + /* fall through */ + case 11: c+=((uint32_t)k[10])<<16; + /* fall through */ + case 10: c+=((uint32_t)k[9])<<8; + /* fall through */ + case 9 : c+=k[8]; + /* fall through */ + case 8 : b+=((uint32_t)k[7])<<24; + /* fall through */ + case 7 : b+=((uint32_t)k[6])<<16; + /* fall through */ + case 6 : b+=((uint32_t)k[5])<<8; + /* fall through */ + case 5 : b+=k[4]; + /* fall through */ + case 4 : a+=((uint32_t)k[3])<<24; + /* fall through */ + case 3 : a+=((uint32_t)k[2])<<16; + /* fall through */ + case 2 : a+=((uint32_t)k[1])<<8; + /* fall through */ + case 1 : a+=k[0]; + break; + case 0 : return c; + } + } + + final(a,b,c); + return c; +} diff --git a/lib/superfox/qpc/nhash2.h b/lib/superfox/qpc/nhash2.h new file mode 100644 index 000000000..1864a93fc --- /dev/null +++ b/lib/superfox/qpc/nhash2.h @@ -0,0 +1,22 @@ +#ifndef NHASH_H_ +#define NHASH_H_ + +#ifdef Win32 +#include "win_stdint.h" /* defines uint32_t etc */ +#else +#include +#include /* defines uint32_t etc */ +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +uint32_t nhash( const void * key, uint64_t length, uint32_t initval); + //int nhash(const char *key, int length, int initval); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/lib/superfox/qpc/np_qpc.c b/lib/superfox/qpc/np_qpc.c new file mode 100644 index 000000000..d44e544b8 --- /dev/null +++ b/lib/superfox/qpc/np_qpc.c @@ -0,0 +1,329 @@ +// ------------------------------------------------------------------------------ +// np_qpc.h +// Q-Ary Polar Codes encoding/decoding functions +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . + +#include +#include +#include + +#include "dbgprintf.h" +#include "qpc_fwht.h" +#include "np_qpc.h" + +// static constant / functions + +static const float knorm = 1.0f / QPC_Q; + +static float* pdf_conv(float *dst, float* pdf1, float* pdf2) +{ + // convolution between two pdf + + float fwht_pdf1[QPC_Q]; + float fwht_pdf2[QPC_Q]; + int k; + + qpc_fwht(fwht_pdf1, pdf1); + qpc_fwht(fwht_pdf2, pdf2); + + for (k = 0; k < QPC_Q; k++) + fwht_pdf1[k] *= fwht_pdf2[k]; + + qpc_fwht(dst, fwht_pdf1); + + for (k = 0; k < QPC_Q; k++) + dst[k] *= knorm; + + return dst; +} +static void pdfarray_conv(float* dstarray, float* pdf1array, float* pdf2array, int numrows) +{ + int k; + + // convolutions between rows of pdfs + + for (k = 0; k < numrows; k++) { + pdf_conv(dstarray, pdf1array, pdf2array); + dstarray += QPC_Q; + pdf1array += QPC_Q; + pdf2array += QPC_Q; + } + +} + +static float _pdfuniform[QPC_Q]; +static const float* _pdf_uniform1(); +static const float* _pdf_uniform0(); +typedef const float*(*ptr_pdfuniform)(void); + +static ptr_pdfuniform _ptr_pdf_uniform = _pdf_uniform0; + +static const float* _pdf_uniform1() +{ + return _pdfuniform; +}; +static const float* _pdf_uniform0() +{ + // compute uniform pdf once for all + int k; + for (k = 0; k < QPC_Q; k++) + _pdfuniform[k] = knorm; + + // next call to _qpc_pdfuniform + // will be handled directly by _pqc_pdfuniform1 + _ptr_pdf_uniform = _pdf_uniform1; + + return _pdfuniform; +}; +static const float* pdf_uniform() +{ + return _ptr_pdf_uniform(); +} + +static float * pdf_mul(float *dst, float* pdf1, float* pdf2) +{ + int k; + float v; + float norm = 0; + for (k = 0; k < QPC_Q; k++) { + v = pdf1[k] * pdf2[k]; + dst[k] = v; + norm += v; + } + // if norm of the result is not positive + // return in dst a uniform distribution + if (norm <= 0) + memcpy(dst, pdf_uniform(), QPC_Q * sizeof(float)); + else { + norm = 1.0f / norm; + for (k = 0; k < QPC_Q; k++) + dst[k] = dst[k] * norm; + } + + + return dst; +} +static void pdfarray_mul(float* dstarray, float* pdf1array, float* pdf2array, int numrows) +{ + int k; + + // products between rows of pdfs + + for (k = 0; k < numrows; k++) { + pdf_mul(dstarray, pdf1array, pdf2array); + dstarray += QPC_Q; + pdf1array += QPC_Q; + pdf2array += QPC_Q; + } +} + +static float* pdf_convhard(float* dst, const float* pdf, unsigned char hd) +{ + // convolution between a pdf and a hard-decision feedback + + int k; + for (k=0;k pdfmax) { + pdfmax = pdf[k]; + imax = k; + } + + return imax; +} + +// local stack functions --------------------------------------- +static float _qpc_stack[QPC_N * QPC_Q * 2]; +static float* _qpc_stack_base = _qpc_stack; + +static float* _qpc_stack_push(int numfloats) +{ + float* addr = _qpc_stack_base; + _qpc_stack_base += numfloats; + return addr; +} +static void _qpc_stack_pop(int numfloats) +{ + _qpc_stack_base -= numfloats; +} + +// qpc encoder function (internal use) ---------------------------------------------------------- +unsigned char* _qpc_encode(unsigned char* y, unsigned char* x) +{ + // Non recursive polar encoder + // Same architecture of a fast fourier transform + // in which the fft butteflies are replaced by the polar transform + // butterflies + + int k, j, m; + int groups; + int bfypergroup; + int stepbfy; + int stepgroup; + int basegroup; + int basebfy; + + memcpy(y, x, QPC_N); + + for (k = 0; k < QPC_LOG2N; k++) { + groups = 1 << (QPC_LOG2N - 1 - k); + stepbfy = bfypergroup = 1 << k; + stepgroup = stepbfy << 1; + basegroup = 0; + for (j = 0; j < groups; j++) { + basebfy = basegroup; + for (m = 0; m < bfypergroup; m++) { + // polar transform + y[basebfy + stepbfy] = y[basebfy + stepbfy] ^ y[basebfy]; + basebfy = basebfy + 1; + } + basegroup = basegroup + stepgroup; + } + } + + return y; +} + +// qpc polar decoder (internal use )-------------------------------------------------- +void _qpc_decode(unsigned char* xdec, unsigned char* ydec, const float* py, const unsigned char* f, const unsigned char* fsize, const int numrows) +{ + + if (numrows == 1) { + if (fsize[0] == 0) { + // dbgprintf_vector_float("py", py, QPC_Q); + + // frozen symbol + xdec[0] = pdf_max(py); + ydec[0] = f[0]; + } + else { + // fsize = 1 => information symbol + xdec[0] = pdf_max(py); + ydec[0] = xdec[0]; + } + + } + else { + + int k; + int nextrows = numrows >> 1; + int size = nextrows << QPC_LOG2Q; + + // upper block variables + unsigned char* xdech = xdec + nextrows; + unsigned char* ydech = ydec + nextrows; + const unsigned char* fh = f + nextrows; + const unsigned char* fsizeh = fsize + nextrows; + + // Step 1. + // stack and init variables used in the recursion + + float* pyl = _qpc_stack_push(size); + memcpy(pyl, py, size * sizeof(float)); + + float* pyh = _qpc_stack_push(size); + memcpy(pyh, py + size, size * sizeof(float)); + + // Step 2. Recursion on upper block + // Forward pdf convolutions for the upper block + // (place in the lower part of py the convolution of lower and higher blocks) + +// float* pyh = py + size; +// pdfarray_conv(py, pyl, pyh, nextrows); // convolution overwriting the lower block of py which is not needed + + pdfarray_conv(pyh, pyl, pyh, nextrows); + _qpc_decode(xdech, ydech, pyh, fh, fsizeh, nextrows); + + // Step 3. compute pdfs in the lower block + pdfarray_convhard(pyh, py+size, ydech,nextrows); // dst ptr must be different form src ptr + pdfarray_mul(pyl, pyl, pyh, nextrows); + // we don't need pyh anymore + _qpc_stack_pop(size); + + // Step 4. Recursion on the lower block + _qpc_decode(xdec, ydec, pyl, f, fsize, nextrows); + // we don't need pyl anymore + _qpc_stack_pop(size); + + // Step 5. Update backward results + // xdec is already ok, we need just to update ydech + for (k = 0; k < nextrows; k++) + ydech[k] = ydech[k] ^ ydec[k]; + + } + + +} + +// Public encoding/decoding functions ------------------------------------------------ +void qpc_encode(unsigned char* y, const unsigned char* x) +{ + + // map information symbols + int kk, pos; + for (kk = 0; kk < QPC_K; kk++) { + pos = qpccode.xpos[kk]; + qpccode.f[pos] = x[kk]; + } + // do polar encoding + _qpc_encode(y, qpccode.f); +} +void qpc_decode(unsigned char* xdec, unsigned char* ydec, float* py) +{ + int k; + unsigned char x[QPC_N]; + + // set the first py row with know frozen (punctured) symbol + if (qpccode.np < qpccode.n) { + // assume that we punctured only the first output symbol + memset(py, 0, QPC_Q * sizeof(float)); + py[qpccode.f[0]] = 1.0f; + } + + // decode + _qpc_decode(x, ydec, py, qpccode.f, qpccode.fsize, QPC_N); + + // demap information symbols + for (k = 0; k < QPC_K; k++) + xdec[k] = x[qpccode.xpos[k]]; + +} diff --git a/lib/superfox/qpc/np_qpc.h b/lib/superfox/qpc/np_qpc.h new file mode 100644 index 000000000..0fa440e76 --- /dev/null +++ b/lib/superfox/qpc/np_qpc.h @@ -0,0 +1,62 @@ +// ------------------------------------------------------------------------------ +// np_qpc.h +// Q-Ary Polar Codes encoding/decoding functions +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . + + + +#ifndef _np_qpc_h_ +#define _np_qpc_h_ + +#define QPC_LOG2N 7 // log2(codeword length) (not punctured) +#define QPC_N (1< frozen) +} qpccode_ds; + +#ifdef __cplusplus +extern "C" +{ +#endif + + void qpc_encode(unsigned char* y, const unsigned char* x); + void qpc_decode(unsigned char* xdec, unsigned char* ydec, float* py); + + unsigned char* _qpc_encode(unsigned char* y, unsigned char* x); + void _qpc_decode(unsigned char* xdec, unsigned char* ydec, + const float* py, const unsigned char* f, const unsigned char* fsize, + const int numrows); + + extern qpccode_ds qpccode; + +#ifdef __cplusplus +} +#endif + +#endif // _np_qpc_h_ diff --git a/lib/superfox/qpc/np_rnd.c b/lib/superfox/qpc/np_rnd.c new file mode 100644 index 000000000..df2afacc3 --- /dev/null +++ b/lib/superfox/qpc/np_rnd.c @@ -0,0 +1,135 @@ +// ------------------------------------------------------------------------------ +// np_rnd.c +// Functions to generate random numbers with uniform/gaussian probability distributions +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . + +#include "np_rnd.h" + +#if _WIN32 // note the underscore: without it, it's not msdn official! + // Windows (x64 and x86) + #include // required only for GetTickCount(...) + #define K_RAND_MAX UINT_MAX +#elif _SVID_SOURCE || _XOPEN_SOURCE || __unix__ || (defined (__APPLE__) && defined(__MACH__)) /* POSIX or Unix or Apple */ + #include + #define rand_s(x) (*x)=(unsigned int)lrand48() // returns unsigned integers in the range 0..0x7FFFFFFF + #define K_RAND_MAX 0x7FFFFFFF // that's the max number + // generated by lrand48 +#else + #error "No good quality PRNG found" +#endif + + +void np_normrnd_real(float *dst, int nitems, float mean, float stdev) +{ + // for odd or even nitems + + unsigned int r; + float phi=0, u=0; + int set = 0; + float scalephi = (M_2PI / (1.0f + K_RAND_MAX)); + float scaleu = (1.0f / (1.0f + K_RAND_MAX)); + + + while (nitems--) + if (set==1) { + // generate a second normally distributed number + *dst++ = sinf(phi)*u*stdev+mean; + set = 0; + } + else { + // generate a uniform distributed phase phi in the range [0,2*PI) + rand_s((unsigned int*)&r); phi = scalephi * r; + + // generate a uniform distributed random number u in the range (0,1] + rand_s((unsigned int*)&r); u = scaleu * (1.0f + r); + + // generate normally distributed number + u = (float)sqrt(-2.0f * log(u)); + *dst++ = cosf(phi) * u * stdev + mean; + + set=1; + } +} + +void np_normrnd_cpx(float* dst, int nitems, float mean, float stdev) +{ + // as qpc_normrnd_real, but generates always nitems complex numbers (2*nitems reals) + + unsigned int r; + float phi = 0, u = 0; + // JHT int set = 0; + float scalephi = (M_2PI / (1.0f + K_RAND_MAX)); + float scaleu = (1.0f / (1.0f + K_RAND_MAX)); + + while (nitems--) { + // generate a uniform distributed phase phi in the range [0,2*PI) + rand_s((unsigned int*)&r); phi = scalephi * r; + + // generate a uniform distributed random number u in the range (0,1] + rand_s((unsigned int*)&r); u = scaleu * (1.0f + r); + + // generate normally distributed real and imaginary parts + u = (float)sqrt(-2.0f * log(u)); + *dst++ = cosf(phi) * u * stdev + mean; + *dst++ = sinf(phi) * u * stdev + mean; + } +} + +void np_unidrnd(unsigned int* dst, int nitems, unsigned int nsetsize) +{ + // generate uniform unsigned int random numbers in the range [0..nsetsize) + + unsigned int r; + float u; + float scaleu = (1.0f / (1.0f + K_RAND_MAX)); + + while (nitems--) { + rand_s((unsigned int*)&r); u = scaleu * nsetsize * r; + *dst++ = (unsigned int)floorf(u); + } + +} +void np_unidrnd_uc(unsigned char* dst, int nitems, unsigned char nsetsize) +{ + // generate uniform unsigned char random numbers in the range [0..nsetsize) + + unsigned int r; + float u; + float scaleu = (1.0f / (1.0f + K_RAND_MAX)); + + while (nitems--) { + rand_s((unsigned int*)&r); u = scaleu * nsetsize * r; + *dst++ = (unsigned char)floorf(u); + } + +} + +void np_unifrnd(float* dst, int nitems, float fmax) +{ + // generate uniform float random numbers in the range [0..fmax) + + unsigned int r; + float u; + float scaleu = (1.0f / (1.0f + K_RAND_MAX)); + + while (nitems--) { + rand_s((unsigned int*)&r); u = scaleu * fmax * r; + *dst++ = u; + } + +} diff --git a/lib/superfox/qpc/np_rnd.h b/lib/superfox/qpc/np_rnd.h new file mode 100644 index 000000000..47a71a3c4 --- /dev/null +++ b/lib/superfox/qpc/np_rnd.h @@ -0,0 +1,59 @@ +// ------------------------------------------------------------------------------ +// np_rnd.h +// Functions to generate random numbers with uniform/gaussian probability distributions +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . + +#ifndef _np_rnd_h_ +#define _np_rnd_h_ + +#define _CRT_RAND_S +#include + +#define _USE_MATH_DEFINES +#include +#define M_2PI (2.0f*(float)M_PI) + +#ifdef __cplusplus +extern "C" { +#endif + +// generate a random array of real numbers with a gaussian distribution of given mean and stdev +void np_normrnd_real(float *dst, int nitems, float mean, float stdev); + +// generate a random array of nitems complex numbers with a gaussian distribution of given mean and stdev +void np_normrnd_cpx(float* dst, int nitems, float mean, float stdev); + +// generate a random array of nitems unsigned int numbers with uniform distribution +// in the range [0 .. nsetsize) +void np_unidrnd(unsigned int* dst, int nitems, unsigned int nsetsize); + +// generate a random array of nitems unsigned chars numbers with uniform distribution +// in the range [0 .. nsetsize) +void np_unidrnd_uc(unsigned char* dst, int nitems, unsigned char nsetsize); + +// generate a random array of nitems float numbers with uniform distribution +// in the range [0 .. fmax) +void np_unifrnd(float* dst, int nitems, float fmax); + + +#ifdef __cplusplus +} +#endif + +#endif // _np_rnd_h_ + diff --git a/lib/superfox/qpc/qpc_fwht.c b/lib/superfox/qpc/qpc_fwht.c new file mode 100644 index 000000000..774112110 --- /dev/null +++ b/lib/superfox/qpc/qpc_fwht.c @@ -0,0 +1,233 @@ +// ------------------------------------------------------------------------------ +// qpc_fwht.c +// Fast Walsh-Hadamard Transforms for q-ary polar codes +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . + +#include "qpc_fwht.h" + +static void _qpc_sumdiff8_16(float* y, float* t) +{ + y[0] = t[0] + t[16]; + y[16] = t[0] - t[16]; + + y[1] = t[1] + t[17]; + y[17] = t[1] - t[17]; + + y[2] = t[2] + t[18]; + y[18] = t[2] - t[18]; + + y[3] = t[3] + t[19]; + y[19] = t[3] - t[19]; + + y[4] = t[4] + t[20]; + y[20] = t[4] - t[20]; + + y[5] = t[5] + t[21]; + y[21] = t[5] - t[21]; + + y[6] = t[6] + t[22]; + y[22] = t[6] - t[22]; + + y[7] = t[7] + t[23]; + y[23] = t[7] - t[23]; + +} +static void _qpc_sumdiff8_32(float* y, float* t) +{ + y[0] = t[0] + t[32]; + y[32] = t[0] - t[32]; + + y[1] = t[1] + t[33]; + y[33] = t[1] - t[33]; + + y[2] = t[2] + t[34]; + y[34] = t[2] - t[34]; + + y[3] = t[3] + t[35]; + y[35] = t[3] - t[35]; + + y[4] = t[4] + t[36]; + y[36] = t[4] - t[36]; + + y[5] = t[5] + t[37]; + y[37] = t[5] - t[37]; + + y[6] = t[6] + t[38]; + y[38] = t[6] - t[38]; + + y[7] = t[7] + t[39]; + y[39] = t[7] - t[39]; + +} +static void _qpc_sumdiff8_64(float* y, float* t) +{ + y[0] = t[0] + t[64]; + y[64] = t[0] - t[64]; + + y[1] = t[1] + t[65]; + y[65] = t[1] - t[65]; + + y[2] = t[2] + t[66]; + y[66] = t[2] - t[66]; + + y[3] = t[3] + t[67]; + y[67] = t[3] - t[67]; + + y[4] = t[4] + t[68]; + y[68] = t[4] - t[68]; + + y[5] = t[5] + t[69]; + y[69] = t[5] - t[69]; + + y[6] = t[6] + t[70]; + y[70] = t[6] - t[70]; + + y[7] = t[7] + t[71]; + y[71] = t[7] - t[71]; + +} + +float* qpc_fwht8(float* y, float* x) +{ + float t[8]; + + // first stage + y[0] = x[0] + x[1]; + y[1] = x[0] - x[1]; + + y[2] = x[2] + x[3]; + y[3] = x[2] - x[3]; + + y[4] = x[4] + x[5]; + y[5] = x[4] - x[5]; + + y[6] = x[6] + x[7]; + y[7] = x[6] - x[7]; + + // second stage + t[0] = y[0] + y[2]; + t[2] = y[0] - y[2]; + + t[1] = y[1] + y[3]; + t[3] = y[1] - y[3]; + + t[4] = y[4] + y[6]; + t[6] = y[4] - y[6]; + + t[5] = y[5] + y[7]; + t[7] = y[5] - y[7]; + + // third stage + y[0] = t[0] + t[4]; + y[4] = t[0] - t[4]; + + y[1] = t[1] + t[5]; + y[5] = t[1] - t[5]; + + y[2] = t[2] + t[6]; + y[6] = t[2] - t[6]; + + y[3] = t[3] + t[7]; + y[7] = t[3] - t[7]; + + return y; +} +float* qpc_fwht16(float* y, float* x) +{ + float t[16]; + + qpc_fwht8(t, x); + qpc_fwht8(t + 8, x + 8); + + y[0] = t[0] + t[8]; + y[8] = t[0] - t[8]; + + y[1] = t[1] + t[9]; + y[9] = t[1] - t[9]; + + y[2] = t[2] + t[10]; + y[10] = t[2] - t[10]; + + y[3] = t[3] + t[11]; + y[11] = t[3] - t[11]; + + y[4] = t[4] + t[12]; + y[12] = t[4] - t[12]; + + y[5] = t[5] + t[13]; + y[13] = t[5] - t[13]; + + y[6] = t[6] + t[14]; + y[14] = t[6] - t[14]; + + y[7] = t[7] + t[15]; + y[15] = t[7] - t[15]; + + return y; + +} +float* qpc_fwht32(float* y, float* x) +{ + float t[32]; + + qpc_fwht16(t, x); + qpc_fwht16(t + 16, x + 16); + + _qpc_sumdiff8_16(y, t); + _qpc_sumdiff8_16(y + 8, t + 8); + + return y; +} +float* qpc_fwht64(float* y, float* x) +{ + float t[64]; + + qpc_fwht32(t, x); + qpc_fwht32(t + 32, x + 32); + + _qpc_sumdiff8_32(y, t); + _qpc_sumdiff8_32(y + 8, t + 8); + _qpc_sumdiff8_32(y + 16, t + 16); + _qpc_sumdiff8_32(y + 24, t + 24); + + return y; +} +float* qpc_fwht128(float* y, float* x) +{ + float t[128]; + + qpc_fwht64(t, x); + qpc_fwht64(t + 64, x + 64); + + _qpc_sumdiff8_64(y, t); + _qpc_sumdiff8_64(y + 8, t + 8); + _qpc_sumdiff8_64(y + 16, t + 16); + _qpc_sumdiff8_64(y + 24, t + 24); + _qpc_sumdiff8_64(y + 32, t + 32); + _qpc_sumdiff8_64(y + 40, t + 40); + _qpc_sumdiff8_64(y + 48, t + 48); + _qpc_sumdiff8_64(y + 56, t + 56); + + return y; +} + +// functions over pdfs used by the decoder ----------------------------------------------------------- + + + + diff --git a/lib/superfox/qpc/qpc_fwht.h b/lib/superfox/qpc/qpc_fwht.h new file mode 100644 index 000000000..30477a331 --- /dev/null +++ b/lib/superfox/qpc/qpc_fwht.h @@ -0,0 +1,57 @@ +// ------------------------------------------------------------------------------ +// qpc_fwht.h +// Fast Walsh-Hadamard Transforms for q-ary polar codes +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// This source 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 file 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 source distribution. +// If not, see . + +#ifndef _qpc_fwht_h_ +#define _qpc_fwht_h_ + +#include "np_qpc.h" + +#if QPC_LOG2Q==7 +#define qpc_fwht(a,b) qpc_fwht128(a,b) +#endif +#if QPC_LOG2Q==6 +#define qpc_fwht(a,b) qpc_fwht64(a,b) +#endif +#if QPC_LOG2Q==5 +#define qpc_fwht(a,b) qpc_fwht32(a,b) +#endif +#if QPC_LOG2Q==4 +#define qpc_fwht(a,b) qpc_fwht16(a,b) +#endif +// Note that it makes no sense to use fast convolutions +// for transforms that are less than 16 symbols in size. +// For such cases direct convolutions are faster. +#if QPC_LOG2Q==3 +#define qpc_fwht(a,b) qpc_fwht8(a,b) +#endif + +#ifdef __cplusplus +extern "C" { +#endif + float* qpc_fwht8(float* y, float* x); + float* qpc_fwht16(float* y, float* x); + float* qpc_fwht32(float* y, float* x); + float* qpc_fwht64(float* y, float* x); + float* qpc_fwht128(float* y, float* x); +#ifdef __cplusplus +} +#endif + +#endif // _qpc_fwht_h_ diff --git a/lib/superfox/qpc/qpc_main.cpp b/lib/superfox/qpc/qpc_main.cpp new file mode 100644 index 000000000..8e6c59556 --- /dev/null +++ b/lib/superfox/qpc/qpc_main.cpp @@ -0,0 +1,198 @@ + +// ------------------------------------------------------------------------------ +// qpc_main.cpp +// +// Test the WER performance of the QPC (127,50) Q=128 code +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ + +#include +#include +#include + +#include "dbgprintf.h" +#include "qpc_fwht.h" +#include "np_qpc.h" + +// for random numbers generation and NCFSK Rayleigh channel simulation +#include "np_rnd.h" +void qpc_channel(float* yout, unsigned char* y, float EsNo) +{ + // compute channel outputs amplitudes + + int k; + + // generate noise samples ----------------------------------------- + float sigman = 1.0f / sqrtf(2.0f); + np_normrnd_cpx(yout, QPC_N * QPC_Q, 0, sigman); + + // dbgprintf_rows_float("yout", yout, QPC_Q*2, QPC_N); + + // generate rayleigh distributed signal amplitudes ----------------- + float symamps[QPC_N * 2]; + np_normrnd_cpx(symamps, QPC_N, 0, sigman); + + // dbgprintf_vector_float("symamps", symamps, QPC_N*2); + + // normalize signal amps to unity ---------------------------------- + float pwr = 0.0f; + float* psymamps = symamps; + // compute sig power + for (k = 0; k < QPC_N; k++) { + pwr += psymamps[0] * psymamps[0] + psymamps[1] * psymamps[1]; + psymamps += 2; + } + pwr = pwr / QPC_N; + + // normalize to avg EsNo + float norm = sqrtf(EsNo/pwr); + psymamps = symamps; + for (k = 0; k < QPC_N; k++) { + psymamps[0] *= norm; + psymamps[1] *= norm; + psymamps += 2; + } + // dbgprintf_vector_float("symamps norm", symamps, QPC_N * 2); + + // add signal amplitudes to noise ----------------------------------- + float *pyout = yout; + psymamps= symamps; + for (k = 0; k < QPC_N; k++) { + pyout[y[k]<<1] += psymamps[0]; + pyout[(y[k]<<1)+1] += psymamps[1]; + pyout += (QPC_Q*2); + psymamps += 2; + } + + // dbgprintf_rows_float("s+n out", yout, QPC_Q * 2, QPC_N); + + return; +} +void qpc_likelihoods(float* py, float* yout, float EsNo, float No) +{ + // compute symbols likelihoods + // (rayleigh channel) + + int k, j; + + float norm = EsNo / (EsNo + 1) / No; + + // compute likelihoods from energies ---------------------------- + float* pybase; + float* ppyout = yout; + float normpwr; + float normpwrmax; + float pynorm; + + for (k = 0; k < QPC_N; k++) { + + // compute loglikelihoods and largest one from energies + pybase = py + k * QPC_Q; + normpwrmax = 0.0f; + for (j = 0; j < QPC_Q; j++) { + normpwr = norm * (ppyout[0] * ppyout[0] + ppyout[1] * ppyout[1]); + pybase[j] = normpwr; + if (normpwr > normpwrmax) + normpwrmax = normpwr; + ppyout += 2; + } + // subtract largest exponent + pynorm = 0.0f; + for (j = 0; j < QPC_Q; j++) { + pybase[j] = expf(pybase[j] - normpwrmax); + pynorm += pybase[j]; + } + // normalize to probabilities + pynorm = 1.0f / pynorm; + for (j = 0; j < QPC_Q; j++) + pybase[j] = pybase[j] * pynorm; + + } + + return; +} + +int main() +{ + int k; + + int NT = 20000; // number of transmissions to simultate + float EbNodB = 3.7f; // Eb/No to simulate + + + int kc = qpccode.k; + int np = qpccode.np; + float Rc = 1.0f * kc / np; + float Tm = 10.0f; // message duration + float Rb = 1.0f / Tm * (QPC_K - 2) * QPC_LOG2Q; // Bit rate assuming two symbols for CRC + + float EsNodB = EbNodB + 10.0f * log10f(Rc * QPC_LOG2Q); + + static unsigned char xin[QPC_K]; // input information symbols + static unsigned char xdec[QPC_K]; // decoded information symbols + static unsigned char y[QPC_N]; // encoded codeword + static unsigned char ydec[QPC_N]; // decoded codeword + + static float yout[QPC_N * QPC_Q * 2]; // channel complex amplitutes output + static float py[QPC_N * QPC_Q]; // channel output probabilities + + float EsNo = powf(10.0f, EsNodB / 10.0f); + float No = 1.0f; + + //JHT static float we; + static float wer; + + int kk; + + printf("QPC Word Error Rate Test\n"); + printf("Polar Code (%d,%d) Q=%d for the Rayleigh channel\n", qpccode.np, qpccode.k, qpccode.q); + printf("Eb/No = %.1f dB\n", EbNodB); + printf("Es/No = %.1f dB\n", EsNodB); + printf("Tmsg = %.1f s\n", Tm); + printf("Bit Rate = %.1f bit/s\n", Rb); + printf("SNR(2500) = %.1f dB\n\n", EbNodB + 10.0f * log10f(Rb/ 2500)); + + + for (k = 0; k < NT; k++) { + + np_unidrnd_uc(xin, QPC_K, QPC_Q); // generate random information symbols + + qpc_encode(y, xin); // encode + + // compute channel outputs and probabilities + + // Note that if the code is punctured (i.e. N=127) + // the first codeword symbol must not be transmitted over the channel + // Here, in order to avoid useless copies, + // the vector of likelihoods py continues to be a QPC_N*QPC_Q vector of floats. + // The first received symbol likelihoods should start always at offset QPC_Q. + // The first QPC_Q positions of py will be ignored (and set) by the decoder. + qpc_channel(yout, y, EsNo); + + // The decoder can't estimate the EsNo easily + // so we assume that in any case it is 5 dB + float EsNoDec = 3.16; + qpc_likelihoods(py, yout, EsNoDec, No); + + qpc_decode(xdec, ydec, py); // decode + + // count words in errors + for (kk = 0; kk < QPC_K; kk++) + if (xin[kk] ^ xdec[kk]) { + wer += 1; + break; + } + + // show the decoder performance every 200 codewords transmissions + if ((k % 200) == 0 && k>0) { + printf("k=%6d/%6d wer = %8.2E\n", k, NT, wer / k); + fflush(stdout); + } + + } + // show the final result + printf("k=%6d/%6d wer = %8.2E\n", k, NT, wer / k); + printf("\nAll done.\n"); + return 0; +} diff --git a/lib/superfox/qpc/qpc_mod.f90 b/lib/superfox/qpc/qpc_mod.f90 new file mode 100644 index 000000000..931e5b3d1 --- /dev/null +++ b/lib/superfox/qpc/qpc_mod.f90 @@ -0,0 +1,38 @@ +module qpc_mod + interface + integer(c_int32_t) function nhash2 (xin, length, initval) bind(C, & + name="nhash2") + use iso_c_binding, only: c_signed_char, c_int64_t, c_int32_t + integer(c_int64_t), intent(in), value :: length + integer(c_signed_char), intent(in) :: xin(length) + integer(c_int32_t), intent(in), value :: initval + end function nhash2 + end interface + + interface + subroutine qpc_encode(y, xin) bind(C,name="qpc_encode") + use iso_c_binding, only: c_ptr, c_signed_char + integer(c_signed_char), intent(out) :: y(127) + integer(c_signed_char), intent(in) :: xin(50) + end subroutine qpc_encode + end interface + + interface + subroutine qpc_channel(yout, y, EsNo) bind(C,name="qpc_channel") + use iso_c_binding, only: c_float_complex, c_float, c_signed_char + complex(c_float_complex), intent(out) :: yout(128,128) + integer(c_signed_char), intent(out) :: y(127) + real(c_float), intent(in), value :: EsNo + end subroutine qpc_channel + end interface + + interface + subroutine qpc_decode(xdec, ydec, py) bind(C,name="qpc_decode") + use iso_c_binding, only: c_float, c_signed_char + real(c_float), intent(in) :: py(128,128) + integer(c_signed_char), intent(out) :: ydec(127) + integer(c_signed_char), intent(out) :: xdec(50) + end subroutine qpc_decode + end interface + +end module qpc_mod diff --git a/lib/superfox/qpc/qpc_n127k50q128.c b/lib/superfox/qpc/qpc_n127k50q128.c new file mode 100644 index 000000000..3897cc833 --- /dev/null +++ b/lib/superfox/qpc/qpc_n127k50q128.c @@ -0,0 +1,60 @@ +// ------------------------------------------------------------------------------ +// qpc_n127k50q128.c +// +// Defines the parameters of the q-ary polar code (127,50) Q=128 +// for WSJT-X +// ------------------------------------------------------------------------------ +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// +// WARNING: +// This source is NOT free software and it is licensed only for use with WSJT-X. +// No derived work is authorized to use this code without the written +// permission of his author (Nico Palermo / IV3NWV) who owns all the rights +// of this IP (intellectual property). +// +// This file is distributed ONLY for the purpose of documenting and making of +// public domain the encoding scheme used in WSJT-X so that the transmitted +// messages can be decoded by anybody. +// Anyway this does not imply that one could use the following tables in a +// derived work without an explicit and written authorization from his author. +// Any unauthorized use, as for any intellectual property, is simply +// illegal. +// ------------------------------------------------------------------------------- +#include "np_qpc.h" +qpccode_ds qpccode = { + 128, //n + 127, //np + 50, //k + 128, //q +{ + 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 16, 32, 17, 18, 64, 20, + 33, 34, 24, 7, 11, 36, 13, 19, 14, 65, 40, 21, 66, 22, 35, 68, + 25, 48, 37, 26, 72, 15, 38, 28, 41, 67, 23, 80, 42, 69, 49, 96, + 44, 27, 70, 50, 73, 39, 29, 52, 74, 30, 56, 81, 76, 43, 82, 84, + 97, 45, 71, 88, 98, 46, 100, 51, 104, 53, 75, 112, 54, 57, 99, 119, + 92, 77, 58, 117, 59, 83, 106, 31, 85, 108, 115, 116, 122, 125, 124, 91, + 61, 90, 89, 111, 78, 93, 94, 126, 86, 107, 110, 118, 121, 62, 120, 87, +105, 55, 114, 60, 127, 63, 103, 101, 123, 95, 102, 47, 109, 79, 113, 0 +}, +{ + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 +}, +{ + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, + 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 +} +}; diff --git a/lib/superfox/qpc/qpc_subs.c b/lib/superfox/qpc/qpc_subs.c new file mode 100644 index 000000000..ec5f8e446 --- /dev/null +++ b/lib/superfox/qpc/qpc_subs.c @@ -0,0 +1,115 @@ +// ------------------------------------------------------------------------------ +// qpc_test.cpp +// +// Test the WER performance of the QPC (127,50) Q=128 code +// +// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ + +#include +#include +#include +#include + +#include "dbgprintf.h" +#include "qpc_fwht.h" +#include "np_qpc.h" + +#include "nhash2.h" + +// for random numbers generation and NCFSK Rayleigh channel simulation +#include "np_rnd.h" +void qpc_channel(float* yout, unsigned char* y, float EsNo) +{ + // compute channel outputs amplitudes + + int k; + + // generate noise samples ----------------------------------------- + float sigman = 1.0f / sqrtf(2.0f); + np_normrnd_cpx(yout, QPC_N * QPC_Q, 0, sigman); + + // dbgprintf_rows_float("yout", yout, QPC_Q*2, QPC_N); + + // generate rayleigh distributed signal amplitudes ----------------- + float symamps[QPC_N * 2]; + np_normrnd_cpx(symamps, QPC_N, 0, sigman); + + // dbgprintf_vector_float("symamps", symamps, QPC_N*2); + + // normalize signal amps to unity ---------------------------------- + float pwr = 0.0f; + float* psymamps = symamps; + // compute sig power + for (k = 0; k < QPC_N; k++) { + pwr += psymamps[0] * psymamps[0] + psymamps[1] * psymamps[1]; + psymamps += 2; + } + pwr = pwr / QPC_N; + + // normalize to avg EsNo + float norm = sqrtf(EsNo/pwr); + psymamps = symamps; + for (k = 0; k < QPC_N; k++) { + psymamps[0] *= norm; + psymamps[1] *= norm; + psymamps += 2; + } + // dbgprintf_vector_float("symamps norm", symamps, QPC_N * 2); + + // add signal amplitudes to noise ----------------------------------- + float *pyout = yout; + psymamps= symamps; + for (k = 0; k < QPC_N; k++) { + pyout[y[k]<<1] += psymamps[0]; + pyout[(y[k]<<1)+1] += psymamps[1]; + pyout += (QPC_Q*2); + psymamps += 2; + } + + // dbgprintf_rows_float("s+n out", yout, QPC_Q * 2, QPC_N); + + return; +} +void qpc_likelihoods(float* py, float* yout, float EsNo, float No) +{ + // compute symbols likelihoods + // (rayleigh channel) + + int k, j; + + float norm = EsNo / (EsNo + 1) / No; + + // compute likelihoods from energies ---------------------------- + float* pybase; + float* ppyout = yout; + float normpwr; + float normpwrmax; + float pynorm; + + for (k = 0; k < QPC_N; k++) { + + // compute loglikelihoods and largest one from energies + pybase = py + k * QPC_Q; + normpwrmax = 0.0f; + for (j = 0; j < QPC_Q; j++) { + normpwr = norm * (ppyout[0] * ppyout[0] + ppyout[1] * ppyout[1]); + pybase[j] = normpwr; + if (normpwr > normpwrmax) + normpwrmax = normpwr; + ppyout += 2; + } + // subtract largest exponent + pynorm = 0.0f; + for (j = 0; j < QPC_Q; j++) { + pybase[j] = expf(pybase[j] - normpwrmax); + pynorm += pybase[j]; + } + // normalize to probabilities + pynorm = 1.0f / pynorm; + for (j = 0; j < QPC_Q; j++) + pybase[j] = pybase[j] * pynorm; + } + + return; +} diff --git a/lib/superfox/qpc_decode2.f90 b/lib/superfox/qpc_decode2.f90 new file mode 100644 index 000000000..2a260e021 --- /dev/null +++ b/lib/superfox/qpc_decode2.f90 @@ -0,0 +1,218 @@ +subroutine qpc_decode2(c0,fsync,ftol,xdec,ndepth,dth,damp,crc_ok, & + snrsync,fbest,tbest,snr) + + use qpc_mod + + parameter(NMAX=15*12000,NFT=365,NZ=100) + complex c0(NMAX) !Signal as received + complex c(NMAX) !Signal as received + real py(0:127,0:127) !Probabilities for received synbol values + real py0(0:127,0:127) !Probabilities for strong signal + real pyd(0:127,0:127) !Dithered values for py + real s2(0:127,0:151) !Symbol spectra, including sync + real s3(0:127,0:127) !Synchronized symbol spectra + real No + integer crc_chk,crc_sent + integer*8 n47 + integer idf(NZ),idt(NZ) + integer nseed(33) + integer*1 xdec(0:49) !Decoded message + integer*1 ydec(0:127) !Decoded symbols + logical crc_ok + integer maxdither(8) + integer isync(24) !Symbol numbers for sync tones + data isync/1,2,4,7,11,16,22,29,37,39,42,43,45,48,52,57,63,70,78,80,83, & + 84,86,89/ + data n47/47/,maxdither/20,50,100,200,500,1000,2000,5000/ + data nseed/ & + 321278106, -658879006, 1239150429, -941466001, -698554454, & + 1136210962, 1633585627, 1261915021, -1134191465, -487888229, & + 2131958895, -1429290834, -1802468092, 1801346659, 1966248904, & + 402671397, -1961400750, -1567227835, 1895670987, -286583128, & + -595933665, -1699285543, 1518291336, 1338407128, 838354404, & + -2081343776, -1449416716, 1236537391, -133197638, 337355509, & + -460640480, 1592689606, 0/ + + data idf/0, 0, -1, 0, -1, 1, 0, -1, 1, -2, 0, -1, 1, -2, 2, & + 0, -1, 1, -2, 2, -3, 0, -1, 1, -2, 2, -3, 3, 0, -1, & + 1, -2, 2, -3, 3, -4, 0, -1, 1, -2, 2, -3, 3, -4, 4, & + 0, -1, 1, -2, 2, -3, 3, -4, 4, -5, -1, 1, -2, 2, -3, & + 3, -4, 4, -5, 1, -2, 2, -3, 3, -4, 4, -5, -2, 2, -3, & + 3, -4, 4, -5, 2, -3, 3, -4, 4, -5, -3, 3, -4, 4, -5, & + 3, -4, 4, -5, -4, 4, -5, 4, -5, -5/ + data idt/0 , -1, 0, 1, -1, 0, -2, 1, -1, 0, 2, -2, 1, -1, 0, & + -3, 2, -2, 1, -1, 0, 3, -3, 2, -2, 1, -1, 0, -4, 3, & + -3, 2, -2, 1, -1, 0, 4, -4, 3, -3, 2, -2, 1, -1, 0, & + -5, 4, -4, 3, -3, 2, -2, 1, -1, 0, -5, 4, -4, 3, -3, & + 2, -2, 1, -1, -5, 4, -4, 3, -3, 2, -2, 1, -5, 4, -4, & + 3, -3, 2, -2, -5, 4, -4, 3, -3, 2, -5, 4, -4, 3, -3, & + -5, 4, -4, 3, -5, 4, -4, -5, 4, -5/ + + + fsample=12000.0 + baud=12000.0/1024.0 + nstype=1 + n47=47 + mask21=2**21 - 1 + crc_ok=.false. + + call qpc_sync(c0,fsample,isync,fsync,ftol,f2,t2,snrsync) + f00=1500.0 + f2 + t00=t2 + fbest=f00 + tbest=t00 + maxd=1 + if(ndepth.gt.0) maxd=maxdither(ndepth) + maxft=NZ + if(snrsync.lt.4.0 .or. ndepth.le.0) maxft=1 + do idith=1,maxft + if(idith.ge.2) maxd=1 + deltaf=idf(idith)*0.5 + deltat=idt(idith)*8.0/1024.0 + f=f00+deltaf + t=t00+deltat + fshift=1500.0 - (f+baud) !Shift frequencies down by f + 1 bin + call twkfreq2(c0,c,NMAX,fsample,fshift) + a=1.0 + b=0.0 + do kk=1,4 + if(kk.eq.2) b=0.4 + if(kk.eq.3) b=0.5 + if(kk.eq.4) b=0.6 + call sfox_demod(c,1500.0,t,isync,s2,s3) !Compute s2 and s3 + + if(b.gt.0.0) then + do j=0,127 + call smo121a(s3(:,j),128,a,b) + enddo + endif + call pctile(s3,128*128,50,base3) + s3=s3/base3 + + EsNoDec=3.16 + No=1. + py0=s3 + call qpc_likelihoods2(py,s3,EsNoDec,No) !For weak signals + + call random_seed(put=nseed) + do kkk=1,maxd + if(kkk.eq.1) then + pyd=py0 + else + pyd=0. + if(kkk.gt.2) then + call random_number(pyd) + pyd=2.0*(pyd-0.5) + endif + where(py.gt.dth) pyd=0. !Don't perturb large likelihoods + pyd=py*(1.0 + damp*pyd) !Compute dithered likelihood + endif + do j=0,127 + ss=sum(pyd(:,j)) + if(ss.gt.0.0) then + pyd(:,j)=pyd(:,j)/ss + else + pyd(:,j)=0.0 + endif + enddo + + call qpc_decode(xdec,ydec,pyd) + xdec=xdec(49:0:-1) + crc_chk=iand(nhash2(xdec,n47,571),mask21) !Compute crc_chk + crc_sent=128*128*xdec(47) + 128*xdec(48) + xdec(49) + crc_ok=crc_chk.eq.crc_sent + + if(crc_ok) then + call qpc_snr(s3,ydec,snr) + if(snr.lt.-16.5) crc_ok=.false. +! write(61,3061) idith,kk,kkk,idf(idith),idt(idith),a,b +!3061 format(5i5,2f8.3) + return + endif + enddo !kk: dither of smoothing weights + enddo !kkk: dither of probabilities + enddo !idith: dither of frequency and time + return +end subroutine qpc_decode2 + +subroutine smo121a(x,nz,a,b) + + real x(nz) + fac=1.0/(a+2*b) + x0=x(1) + do i=2,nz-1 + x1=x(i) + x(i)=fac*(a*x(i) + b*(x0+x(i+1))) + x0=x1 + enddo + + return +end subroutine smo121a + +subroutine remove_tone(c0,npts) + + parameter (NFILT=8000) + complex c0(npts) + complex cwindow(15*12000) + complex cref(npts) + complex cfilt(npts) + real window(-NFILT/2:NFILT/2) +! real endcorrection(NFILT/2+1) + real s(npts/4) + integer ipk(1) + logical first + data first/.true./ + save cwindow,first,pi + + if(first) then + pi=4.0*atan(1.0) + fac=1.0/float(npts) + sumw=0.0 + do j=-NFILT/2,NFILT/2 + window(j)=cos(pi*j/NFILT)**2 + sumw=sumw+window(j) + enddo + cwindow=0. + cwindow(1:NFILT+1)=window/sumw + cwindow=cshift(cwindow,NFILT/2+1) + call four2a(cwindow,npts,1,-1,1) + cwindow=cwindow*fac ! frequency domain smoothing filter + first=.false. + endif + + fsample=12000.0 + df=fsample/npts + fac=1.0/npts + cfilt=fac*c0 + call four2a(cfilt,npts,1,-1,1) ! fourier transform of input data + iz=npts/4 + do i=1,iz + s(i)=real(cfilt(i))**2 + aimag(cfilt(i))**2 + enddo + + ia=nint(700.0/df) + ib=nint(800.0/df) + ipk=maxloc(s(ia:ib)) + i0=ipk(1) + ia - 1 + f2=df*i0 + write(*,*) 'remove_tone - frequency: ',f2 + + dt=1.0/fsample + do i=1, npts + arg=2*pi*f2*i*dt + cref(i)=cmplx(cos(arg),sin(arg)) + enddo + cfilt=c0*conjg(cref) ! baseband to be filtered + call four2a(cfilt,npts,1,-1,1) + cfilt=cfilt*cwindow + call four2a(cfilt,npts,1,1,1) + + nframe=50*3456 + do i=1,nframe + cref(i)=cfilt(i)*cref(i) + c0(i)=c0(i)-cref(i) + enddo + + return + +end subroutine remove_tone diff --git a/lib/superfox/qpc_likelihoods2.f90 b/lib/superfox/qpc_likelihoods2.f90 new file mode 100644 index 000000000..70a631814 --- /dev/null +++ b/lib/superfox/qpc_likelihoods2.f90 @@ -0,0 +1,28 @@ +subroutine qpc_likelihoods2(py,s3,EsNo,No) + + integer QQ,QN + parameter(QQ=128,QN=128) + real py(0:QQ-1,0:QN-1) + real s3(0:QQ-1,0:QN-1) + real No,norm,normpwr,normpwrmax + + norm=(EsNo/(EsNo+1.0))/No + +! Compute likelihoods for symbol values, from the symbol power spectra + do k=0,QN-1 + normpwrmax=0. + do j=0,QQ-1 + normpwr=norm*s3(j,k) + py(j,k)=normpwr + normpwrmax=max(normpwr,normpwrmax) + enddo + pynorm=0. + do j=0,QQ-1 + py(j,k)=exp(py(j,k)-normpwrmax) + pynorm=pynorm + py(j,k) + enddo + py(0:QQ-1,k)=py(0:QQ-1,k)/pynorm !Normalize to probabilities + enddo + + return +end subroutine qpc_likelihoods2 diff --git a/lib/superfox/qpc_snr.f90 b/lib/superfox/qpc_snr.f90 new file mode 100644 index 000000000..b2c807f0e --- /dev/null +++ b/lib/superfox/qpc_snr.f90 @@ -0,0 +1,17 @@ +subroutine qpc_snr(s3,y,snr) + + use qpc_mod +! real s2(0:NQ-1,0:151) !Symbol spectra, including sync + real s3(0:127,0:127) !Synchronized symbol spectra + integer*1 y(0:127) !Encoded symbols +! integer isync(24) + + p=0. + do j=1,127 + i=y(j) + p=p + s3(i,j) + enddo + snr=db(p/127.0) - db(127.0) - 4.0 + + return +end subroutine qpc_snr diff --git a/lib/superfox/qpc_sync.f90 b/lib/superfox/qpc_sync.f90 new file mode 100644 index 000000000..ecccc692f --- /dev/null +++ b/lib/superfox/qpc_sync.f90 @@ -0,0 +1,101 @@ +subroutine qpc_sync(crcvd0,fsample,isync,fsync,ftol,f2,t2,snrsync) + + parameter(N9SEC=9*12000,NMAX=15*12000,NDOWN=16,NZ=N9SEC/NDOWN) + complex crcvd0(NMAX) !Signal as received + complex c0(0:N9SEC-1) !For long FFT + complex c1(0:NZ-1) + complex c1sum(0:NZ-1) + complex z + real s(N9SEC/4) + real p(-1125:1125) + integer ipk(1) + integer isync(24) + + baud=12000.0/1024.0 + df2=fsample/N9SEC + fac=1.0/N9SEC + c0=fac*crcvd0(1:N9SEC) + call four2a(c0,N9SEC,1,-1,1) !Forward c2c FFT + iz=N9SEC/4 + do i=1,iz + s(i)=real(c0(i))**2 + aimag(c0(i))**2 + enddo + + do i=1,4 !Smooth the spectrum a bit + call smo121(s,iz) + enddo + + ia=nint((fsync-ftol)/df2) + ib=nint((fsync+ftol)/df2) + ipk=maxloc(s(ia:ib)) + i0=ipk(1) + ia - 1 + f2=df2*i0-750.0 ! f2 is the offset from nominal 750 Hz. + ia=nint(i0-baud/df2) + ib=nint(i0+baud/df2) + s1=0.0 + s0=0.0 + do i=ia,ib + s0=s0+s(i) + s1=s1+(i-i0)*s(i) + enddo + delta=s1/s0 + i0=nint(i0+delta) + f2=i0*df2-750.0 + + c1=0. + ia=nint(i0-baud/df2) + ib=nint(i0+baud/df2) + do i=ia,ib + j=i-i0 + if(j.ge.0) c1(j)=c0(i) + if(j.lt.0) c1(j+NZ)=c0(i) + enddo + call four2a(c1,NZ,1,1,1) !Reverse c2c FFT: back to time domain + + c1sum(0)=c1(0) + do i=1,NZ-1 + c1sum(i)=c1sum(i-1) + c1(i) + enddo + + nspsd=1024/NDOWN + dt=NDOWN/12000.0 + lagmax=1.5/dt + i0=nint(0.5*fsample/NDOWN) !Nominal start time is 0.5 s + pmax=0. + lagpk=0 + do lag=-lagmax,lagmax + sp=0. + do j=1,24 + i1=i0 + (isync(j)-1)*nspsd + lag + i2=i1 + nspsd + if(i1.lt.0 .or. i1.gt.NZ-1) cycle + if(i2.lt.0 .or. i2.gt.NZ-1) cycle + z=c1sum(i2)-c1sum(i1) + sp=sp + real(z)**2 + aimag(z)**2 + enddo + if(sp.gt.pmax) then + pmax=sp + lagpk=lag + endif + p(lag)=sp + enddo + + t2=lagpk*dt + snrsync=0. + sp=0. + sq=0. + nsum=0 + tsym=1024/12000.0 + do lag=-lagmax,lagmax + t=(lag-lagpk)*dt + if(abs(t).lt.tsym) cycle + nsum=nsum+1 + sp=sp + p(lag) + sq=sq + p(lag)*p(lag) + enddo + ave=sp/nsum + rms=sqrt(sq/nsum-ave*ave) + snrsync=(pmax-ave)/rms + + return +end subroutine qpc_sync diff --git a/lib/superfox/sfox_ana.f90 b/lib/superfox/sfox_ana.f90 new file mode 100644 index 000000000..bc9ea6cf1 --- /dev/null +++ b/lib/superfox/sfox_ana.f90 @@ -0,0 +1,17 @@ +subroutine sfox_ana(dd,npts,c0,npts2) + + real dd(npts) !Raw data at 12000 Hz + complex c0(0:npts2-1) !Complex data at 12000 Hz + save + + nfft1=npts + nfft2=nfft1 + fac=2.0/(32767.0*nfft1) + c0(0:npts-1)=fac*dd(1:npts) + call four2a(c0,nfft1,1,-1,1) !Forward c2c FFT + c0(nfft2/2+1:nfft2-1)=0. !Remove negative frequencies + c0(0)=0.5*c0(0) !Scale the DC term to 1/2 + call four2a(c0,nfft2,1,1,1) !Inverse c2c FFT; c0 is analytic sig + + return +end subroutine sfox_ana diff --git a/lib/superfox/sfox_demod.f90 b/lib/superfox/sfox_demod.f90 new file mode 100644 index 000000000..cbb11df71 --- /dev/null +++ b/lib/superfox/sfox_demod.f90 @@ -0,0 +1,79 @@ +subroutine sfox_demod(crcvd,f,t,isync,s2,s3) + + use sfox_mod + complex crcvd(NMAX) !Signal as received + complex c(0:NSPS-1) !Work array, one symbol long + real s2(0:NQ-1,0:151) !Symbol spectra, including sync + real s3(0:NQ-1,0:NN) !Synchronized symbol spectra + integer isync(24) + integer ipk(1) + integer hist1(0:NQ-1),hist2(0:NQ-1) + + j0=nint(12000.0*(t+0.5)) + df=12000.0/NSPS + i0=nint(f/df)-NQ/2 + k2=0 + s2(:,0)=0. !The punctured symbol + s3(:,0)=0. !The punctured symbol + + do n=1,NDS !Loop over all symbols + jb=n*NSPS + j0 + ja=jb-NSPS+1 + k2=k2+1 + if(ja.lt.1 .or. jb.gt.NMAX) cycle + c=crcvd(ja:jb) + call four2a(c,NSPS,1,-1,1) !Compute symbol spectrum + do i=0,NQ-1 + s2(i,k2)=real(c(i0+i))**2 + aimag(c(i0+i))**2 + enddo + enddo + + + call pctile(s2,NQ*151,50,base2) + s2=s2/base2 + + hist1=0 + hist2=0 + do j=0,151 + ipk=maxloc(s2(1:NQ-1,j)) !Find the spectral peak + i=ipk(1)-1 + hist1(i)=hist1(i)+1 + enddo + + hist1(0)=0 !Don't treat sync tone as a birdie + do i=0,123 + hist2(i)=sum(hist1(i:i+3)) + enddo + + ipk=maxloc(hist1) + i1=ipk(1)-1 + m1=maxval(hist1) + ipk=maxloc(hist2) + i2=ipk(1)-1 + m2=maxval(hist2) + if(m1.gt.12) then + do i=0,127 + if(hist1(i).gt.12) then + s2(i,:)=1.0 + endif + enddo + endif + + if(m2.gt.20) then + if(i2.ge.1) i2=i2-1 + if(i2.gt.120) i2=120 + s2(i2:i2+7,:)=1.0 + endif + + k3=0 + do n=1,NDS !Copy symbol spectra from s2 into s3 + if(any(isync(1:NS).eq.n)) cycle !Skip the sync symbols + k3=k3+1 + s3(:,k3)=s2(:,n) + enddo + + call pctile(s3,NQ*NN,50,base3) + s3=s3/base3 + + return +end subroutine sfox_demod diff --git a/lib/superfox/sfox_mod.f90 b/lib/superfox/sfox_mod.f90 new file mode 100644 index 000000000..01fc06f26 --- /dev/null +++ b/lib/superfox/sfox_mod.f90 @@ -0,0 +1,73 @@ +module sfox_mod + + parameter (NMAX=15*12000) !Samples in iwave (180,000) + integer MM,NQ,NN,KK,NS,NDS,NFZ,NSPS,NSYNC,NZ,NFFT1 + real baud,tsym,bw + +contains + subroutine sfox_init(mm0,nn0,kk0,itu,fspread,delay,fsample,ns0) + + character*2 itu + integer isps(54) + integer iloc(1) + data isps/ 896, 960, 972, 980,1000,1008,1024,1029,1050,1080, & + 1120,1125,1134,1152,1176,1200,1215,1225,1250,1260, & + 1280,1296,1323,1344,1350,1372,1400,1440,1458,1470, & + 1500,1512,1536,1568,1575,1600,1620,1680,1701,1715, & + 1728,1750,1764,1792,1800,1875,1890,1920,1944,1960, & + 2000,2016,2025,2048/ + + MM=mm0 !Bits per symbol + NQ=2**MM !Q, number of MFSK tones + NN=nn0 !Codeword length + KK=kk0 !Number of information symbols + NS=ns0 !Number of sync symbols + NDS=NN+NS !Total number of channel symbols + NFZ=3 !First zero + + jsps=nint(12.8*fsample/NDS) + iloc=minloc(abs(isps-jsps)) + NSPS=isps(iloc(1)) !Samples per symbol + NSYNC=NS*NSPS !Samples in sync waveform + NZ=NSPS*NDS !Samples in full Tx waveform + NFFT1=2*NSPS !Length of FFTs for symbol spectra + + baud=fsample/NSPS + tsym=1.0/baud + bw=NQ*baud + + fspread=0.0 + delay=0.0 + if(itu.eq.'LQ') then + fspread=0.5 + delay=0.5 + else if(itu.eq.'LM') then + fspread=1.5 + delay=2.0 + else if(itu.eq.'LD') then + fspread=10.0 + delay=6.0 + else if(itu.eq.'MQ') then + fspread=0.1 + delay=0.5 + else if(itu.eq.'MM') then + fspread=0.5 + delay=1.0 + else if(itu.eq.'MD') then + fspread=1.0 + delay=2.0 + else if(itu.eq.'HQ') then + fspread=0.5 + delay=1.0 + else if(itu.eq.'HM') then + fspread=10.0 + delay=3.0 + else if(itu.eq.'HD') then + fspread=30.0 + delay=7.0 + endif + + return + end subroutine sfox_init + +end module sfox_mod diff --git a/lib/superfox/sfox_prob.f90 b/lib/superfox/sfox_prob.f90 new file mode 100644 index 000000000..c8e9f25c9 --- /dev/null +++ b/lib/superfox/sfox_prob.f90 @@ -0,0 +1,55 @@ +subroutine sfox_prob(s3,rxdat,rxprob,rxdat2,rxprob2) + +! Demodulate the 64-bin spectra for each of 63 symbols in a frame. + +! Parameters +! rxdat most reliable symbol value +! rxdat2 second most likely symbol value +! rxprob probability that rxdat was the transmitted value +! rxprob2 probability that rxdat2 was the transmitted value + + use sfox_mod + implicit real*8 (a-h,o-z) + real*4 s3(0:NQ-1,0:NN-1) + integer rxdat(0:NN-1),rxprob(0:NN-1),rxdat2(0:NN-1),rxprob2(0:NN-1) + + afac=1.1 +! scale=255.999 + scale=2047.999 + +! Compute average spectral value + ave=sum(s3)/(NQ*ND) + i1=1 !Silence warning + i2=1 + +! Compute probabilities for most reliable symbol values + do j=0,NN-1 !Loop over all symbols + s1=-1.e30 + psum=0. + do i=0,NQ-1 !Loop over frequency bins + x=min(afac*s3(i,j)/ave,50.d0) + psum=psum+s3(i,j) + if(s3(i,j).gt.s1) then + s1=s3(i,j) !Find max signal+noise power + i1=i !Find most reliable symbol value + endif + enddo + if(psum.eq.0.0) psum=1.e-6 !Guard against zero signal+noise + + s2=-1.e30 + do i=0,NQ-1 + if(i.ne.i1 .and. s3(i,j).gt.s2) then + s2=s3(i,j) !Second largest signal+noise power + i2=i !Bin number for second largest power + endif + enddo + p1=s1/psum !p1, p2 are symbol metrics for ftrsd + p2=s2/psum + rxdat(j)=i1 + rxdat2(j)=i2 + rxprob(j)=scale*p1 !Scaled probabilities, 0 - 255 + rxprob2(j)=scale*p2 + enddo + + return +end subroutine sfox_prob diff --git a/lib/superfox/sfox_remove_ft8.f90 b/lib/superfox/sfox_remove_ft8.f90 new file mode 100644 index 000000000..ce47ac008 --- /dev/null +++ b/lib/superfox/sfox_remove_ft8.f90 @@ -0,0 +1,229 @@ +subroutine sfox_remove_ft8(dd,npts) + use packjt77 + include 'ft8_params.f90' + parameter (MAXCAND=100) + parameter (NP2=2812) + integer itone(NN) + integer iloc(1), ip(1) + integer icos7(0:6) + integer graymap(0:7) + integer*1 cw(174),apmask(174) + integer*1 message91(91) + integer*1 message77(77) + + character*77 c77 + character*37 msg37 + + real s8(0:7,NN) + real s2(0:7) + real ss(9) + real candidate(3,MAXCAND) + real sbase(NH1) + real dd(npts) + real bmeta(174),bmetd(174) + real llra(174),llrd(174) + complex cd0(0:3199) + complex csymb(32) + complex ctwk(32) + complex cs(0:7,NN) + data icos7/3,1,4,0,6,5,2/ + logical first,newdat + logical one(0:7,0:2) + logical unpk77_success + data first/.true./ + data graymap/0,1,3,2,5,6,4,7/ + save one,first,twopi + + if(first) then + one=.false. + do i=0,7 + do j=0,2 + if(iand(i,2**j).ne.0) one(i,j)=.true. + enddo + enddo + first=.false. + twopi=8.0*atan(1.0) + endif + + fs2=12000.0/NDOWN + dt2=1.0/fs2 + nfa=500 + nfb=2500 + syncmin=1.5 + nfqso=750 + call sync8(dd,npts,nfa,nfb,syncmin,nfqso,MAXCAND,candidate,ncand,sbase) + newdat=.true. + do icand=1,ncand + f1=candidate(1,icand) + xdt=candidate(2,icand) + xbase=10.0**(0.1*(sbase(nint(f1/3.125))-40.0)) + call ft8_downsample(dd,newdat,f1,cd0) + newdat=.false. + i0=nint((xdt+0.5)*fs2) + smax=0.0 + do idt=i0-10,i0+10 + call sync8d(cd0,idt,ctwk,0,sync) + if(sync.gt.smax) then + smax=sync + ibest=idt + endif + enddo + + smax=0.0 + delfbest=0.0 + do ifr=-5,5 !Search over +/- 2.5 Hz + delf=ifr*0.5 + dphi=twopi*delf*dt2 + phi=0.0 + do i=1,32 + ctwk(i)=cmplx(cos(phi),sin(phi)) + phi=mod(phi+dphi,twopi) + enddo + call sync8d(cd0,ibest,ctwk,1,sync) + if( sync .gt. smax ) then + smax=sync + delfbest=delf + endif + enddo + + f1=f1+delfbest !Improved estimate of DF + + call ft8_downsample(dd,newdat,f1,cd0) !Mix f1 to baseband and downsample + + smax=0.0 + do idt=-4,4 !Search over +/- one quarter symbol + call sync8d(cd0,ibest+idt,ctwk,0,sync) + ss(idt+5)=sync + enddo + smax=maxval(ss) + iloc=maxloc(ss) + ibest=iloc(1)-5+ibest + xdt=(ibest-1)*dt2 + sync=smax + + do k=1,NN + i1=ibest+(k-1)*32 + csymb=cmplx(0.0,0.0) + if( i1.ge.0 .and. i1+31 .le. NP2-1 ) csymb=cd0(i1:i1+31) + call four2a(csymb,32,1,-1,1) + cs(0:7,k)=csymb(1:8)/1e3 + s8(0:7,k)=abs(csymb(1:8)) + enddo + +! sync quality check + is1=0 + is2=0 + is3=0 + do k=1,7 + ip=maxloc(s8(:,k)) + if(icos7(k-1).eq.(ip(1)-1)) is1=is1+1 + ip=maxloc(s8(:,k+36)) + if(icos7(k-1).eq.(ip(1)-1)) is2=is2+1 + ip=maxloc(s8(:,k+72)) + if(icos7(k-1).eq.(ip(1)-1)) is3=is3+1 + enddo +! hard sync sum - max is 21 + nsync=is1+is2+is3 + if(nsync .le. 6) then ! bail out + nbadcrc=1 + cycle + endif + + nsym=1 + nt=2**(3*nsym) + do ihalf=1,2 + do k=1,29,nsym + if(ihalf.eq.1) ks=k+7 + if(ihalf.eq.2) ks=k+43 + amax=-1.0 + do i=0,nt-1 + i3=iand(i,7) + s2(i)=abs(cs(graymap(i3),ks)) + enddo + i32=1+(k-1)*3+(ihalf-1)*87 + ibmax=2 + do ib=0,ibmax + bm=maxval(s2(0:nt-1),one(0:nt-1,ibmax-ib)) - & + maxval(s2(0:nt-1),.not.one(0:nt-1,ibmax-ib)) + if(i32+ib .gt.174) cycle + bmeta(i32+ib)=bm + den=max(maxval(s2(0:nt-1),one(0:nt-1,ibmax-ib)), & + maxval(s2(0:nt-1),.not.one(0:nt-1,ibmax-ib))) + if(den.gt.0.0) then + cm=bm/den + else ! erase it + cm=0.0 + endif + bmetd(i32+ib)=cm + enddo + enddo + enddo + + call normalizebmet(bmeta,174) + call normalizebmet(bmetd,174) + + scalefac=2.83 + llra=scalefac*bmeta + llrd=scalefac*bmetd + + cw=0 + dmin=0 + norder=2 + maxosd=-1 + Keff=91 + apmask=0 + message91=0 + cw=0 + call decode174_91(llra,Keff,maxosd,norder,apmask,message91,cw, & + ntype,nharderrors,dmin) + + if(nharderrors.ge.0) then + message77=message91(1:77) + else + cycle + endif + + if(count(cw.eq.0).eq.174) cycle + + write(c77,'(77i1)') message77 + read(c77(72:74),'(b3)') n3 + read(c77(75:77),'(b3)') i3 + if(i3.gt.5 .or. (i3.eq.0.and.n3.gt.6)) cycle + if(i3.eq.0 .and. n3.eq.2) cycle + + call unpack77(c77,1,msg37,unpk77_success) +! write(77,*) 'FT8 interference: ',msg37 + if(.not.unpk77_success) cycle +! Message structure: S7 D29 S7 D29 S7 + itone(1:7)=icos7 + itone(36+1:36+7)=icos7 + itone(NN-6:NN)=icos7 + k=7 + do j=1,ND + i=3*j -2 + k=k+1 + if(j.eq.30) k=k+7 + indx=cw(i)*4 + cw(i+1)*2 + cw(i+2) + itone(k)=graymap(indx) + enddo + call subtractft8(dd,itone,f1,xdt,.true.) +! return + enddo + return +end subroutine sfox_remove_ft8 + + +subroutine normalizebmet(bmet,n) + real bmet(n) + + bmetav=sum(bmet)/real(n) + bmet2av=sum(bmet*bmet)/real(n) + var=bmet2av-bmetav*bmetav + if( var .gt. 0.0 ) then + bmetsig=sqrt(var) + else + bmetsig=sqrt(bmet2av) + endif + bmet=bmet/bmetsig + return +end subroutine normalizebmet diff --git a/lib/superfox/sfox_unpack.f90 b/lib/superfox/sfox_unpack.f90 new file mode 100644 index 000000000..903e92edf --- /dev/null +++ b/lib/superfox/sfox_unpack.f90 @@ -0,0 +1,125 @@ +subroutine sfox_unpack(nutc,x,nsnr,f0,dt0,foxcall,nsignature) + + use packjt77 + integer*1 x(0:49) + integer*8 n58 + logical success + character*336 msgbits + character*22 msg(10) !### only msg(1) is used ??? ### + character*13 foxcall,c13 + character*10 ssignature + character*4 crpt(5),grid4 + character*26 freeTextMsg + character*38 c + logical use_otp + data c/' 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ/'/ + + ncq=0 + if (nsignature.eq.0) then + use_otp = .FALSE. + else + use_otp = .TRUE. + endif + write(msgbits,1000) x(0:46) +1000 format(47b7.7) + read(msgbits(327:329),'(b6)') i3 !Message type + read(msgbits(1:28),'(b28)') n28 !Standard Fox call + call unpack28(n28,foxcall,success) +! print*,'aa',i3,foxcall + + if(i3.eq.1) then !Compound Fox callsign +! read(msgbits(87:101),'(b15)') n15 +! call unpackgrid(n15,grid4) +! msg(1)='CQ '//trim(foxcall)//' '//grid4 +! write(*,1100) nutc,nsnr,dt0,nint(f0),trim(msg(1)) +! go to 100 + else if(i3.eq.2) then !Up to 4 Hound calls and free text + call unpacktext77(msgbits(161:231),freeTextMsg(1:13)) + call unpacktext77(msgbits(232:302),freeTextMsg(14:26)) + do i=26,1,-1 + if(freeTextMsg(i:i).ne.'.') exit + freeTextMsg(i:i)=' ' + enddo +! print*,'aa1 ',freeTextMsg + write(*,1100) nutc,nsnr,dt0,nint(f0),freeTextMsg +1100 format(i6.6,i4,f5.1,i5,1x,"~",2x,a) + else if(i3.eq.3) then !CQ FoxCall Grid + read(msgbits(1:58),'(b58)') n58 !FoxCall + do i=11,1,-1 + j=mod(n58,38)+1 + foxcall(i:i)=c(j:j) + n58=n58/38 + enddo + foxcall(12:13)=' ' + read(msgbits(59:73),'(b15)') n15 + call unpackgrid(n15,grid4) + msg(1)='CQ '//trim(foxcall)//' '//grid4 + write(*,1100) nutc,nsnr,dt0,nint(f0),trim(msg(1)) + call unpacktext77(msgbits(74:144),freeTextMsg(1:13)) + call unpacktext77(msgbits(145:215),freeTextMsg(14:26)) + do i=26,1,-1 + if(freeTextMsg(i:i).ne.'.') exit + freeTextMsg(i:i)=' ' + enddo + if(len(trim(freeTextMsg)).gt.0) write(*,1100) nutc,nsnr,dt0,& + nint(f0),freeTextMsg + go to 100 + endif + + j=281 + iz=4 !Max number of reports + if(i3.eq.2) j=141 +! print*,'aa2',j,iz + do i=1,iz !Extract the reports + read(msgbits(j:j+4),'(b5)') n + if(n.eq.31) then + crpt(i)='RR73' + else + write(crpt(i),1006) n-18 +1006 format(i3.2) + if(crpt(i)(1:1).eq.' ') crpt(i)(1:1)='+' + endif +! print*,'aa3',i,n,crpt(i) + j=j+5 + enddo + +! Unpack Hound callsigns and format user-level messages: + iz=9 !Max number of hound calls + if(i3.eq.2 .or. i3.eq.3) iz=4 +! print*,'bb',i3,iz +! print*,msgbits + do i=1,iz + j=28*i + 1 + read(msgbits(j:j+27),'(b28)') n28 + call unpack28(n28,c13,success) +! print*,'cc',i,j,n28 + if(n28.eq.0) cycle + msg(i)=trim(c13)//' '//trim(foxcall) + if(msg(i)(1:3).eq.'CQ ') then + ncq=ncq+1 + else + if(i3.eq.2) then +! print*,'dd',i,crpt(i) + msg(i)=trim(msg(i))//' '//crpt(i) + else + if(i.le.5) msg(i)=trim(msg(i))//' RR73' + if(i.gt.5) msg(i)=trim(msg(i))//' '//crpt(i-5) + endif + endif + if(ncq.le.1 .or. msg(i)(1:3).ne.'CQ ') then + write(*,1100) nutc,nsnr,dt0,nint(f0),trim(msg(i)) + endif + enddo + + if(msgbits(306:306).eq.'1' .and. ncq.lt.1) then + write(*,1100) nutc,nsnr,dt0,nint(f0),'CQ '//foxcall + endif + +100 read(msgbits(307:326),'(b20)') nsignature +! print*,'i3:',i3 + if (use_otp) then + write(ssignature,'(I0)') nsignature + write(*,1100) nutc,nsnr,dt0,nint(f0),'$VERIFY$ '//trim(foxcall)//' '//trim(ssignature) + endif + return +end subroutine sfox_unpack diff --git a/lib/superfox/sfrx.f90 b/lib/superfox/sfrx.f90 new file mode 100644 index 000000000..f73ad0d3b --- /dev/null +++ b/lib/superfox/sfrx.f90 @@ -0,0 +1,137 @@ +program sfrx + + use sfox_mod + use julian + use popen_module, only: get_command_as_string + + integer*2 iwave(NMAX) + integer ihdr(11) + integer*8 secday,ntime8,ntime8_chk + integer*1 xdec(0:49) + character*120 fname,cmnd,cresult + character*13 foxcall,foxcall_chk + character*256 ppath + complex c0(NMAX) !Complex form of signal as received + real dd(NMAX) + logical crc_ok + logical use_otp + data secday/86400/ + include 'gtag.f90' + + use_otp = .FALSE. + fsync=750.0 + ftol=100.0 + narg=iargc() + + if(narg.lt.1) then +! print*,'Usage: sfrx [fsync [ftol]] infile [...]' +! print*,'Examples: sfrx 230305_011230.wav' +! print*,' sfrx 775 10 240811_102400.wav' +! print*,'Reads one or more .wav files and calls SuperFox decoder on each.' +! print*,'Defaults: fsync=750, ftol=100' + print '(" Git tag: ",z9)',ntag + go to 999 + endif + + ifile1=1 + call getarg(1,fname) + read(fname,*,err=1) fsync + ifile1=2 + call getarg(2,fname) + read(fname,*,err=1) ftol + ifile1=3 + +1 nf=0 + nd=0 + nv=0 + + fsample=12000.0 + call sfox_init(7,127,50,'no',fspread,delay,fsample,24) + + npts=15*12000 + + do ifile=ifile1,narg + call getarg(ifile,fname) + if(fname.eq.'OTP') then + use_otp = .TRUE. + cycle + endif + open(10,file=trim(fname),status='old',access='stream',err=4) + + go to 5 +4 print*,'Cannot open file ',trim(fname) + go to 999 +5 read(10) ihdr,iwave + close(10) + + nz=len(trim(fname)) + nyymmdd=ihdr(1) + nutc=ihdr(2) + if(fname(nz-3:nz).eq.'.wav') then + read(fname(nz-16:nz-11),*) nyymmdd + read(fname(nz-9:nz-4),*) nutc + endif + if(nyymmdd.eq.-1) then + ntime8=itime8()/30 + ntime8=30*ntime8 + else + iyr=2000+nyymmdd/10000 + imo=mod(nyymmdd/100,100) + iday=mod(nyymmdd,100) + ih=nutc/10000 + im=mod(nutc/100,100) + is=mod(nutc,100) + ntime8=secday*(JD(iyr,imo,iday)-2440588) + 3600*ih + 60*im + is + endif + + dd=iwave + call sfox_remove_ft8(dd,npts) + + call sfox_ana(dd,npts,c0,npts) + + ndepth=3 + dth=0.5 + damp=1.0 + + call qpc_decode2(c0,fsync,ftol, xdec,ndepth,dth,damp,crc_ok, & + snrsync,fbest,tbest,snr) + if(crc_ok) then + nsnr=nint(snr) + if (use_otp) then + nsignature = 1 + else + nsignature = 0 + endif + + call sfox_unpack(nutc,xdec,nsnr,fbest-750.0,tbest,foxcall,nsignature) +! Execute 'foxchk ' to get the correct signature value. + write(cmnd,1102) trim(foxcall),ntime8 +1102 format('foxchk ',a,i12) + call getarg(0,ppath) + lindex=index(ppath,'sfrx')-1 + if (.not.use_otp) then + cmnd=ppath(1:lindex)//cmnd + cresult=get_command_as_string(trim(cmnd)) + read(cresult,1104) foxcall_chk,ntime8_chk,nsignature_chk +1104 format(a11,i13,i10.7) + if(nsignature.eq.nsignature_chk) write(*,1110) trim(foxcall) +1110 format(a,' verified') + if(nsignature.eq.nsignature_chk) nv=nv+1 + endif + nd=nd+1 +! else +! i0=index(fname,'.wav') +! write(60,3060) trim(fname) +!3060 format('cp ',a,' nodecode') + endif + nf=nf+1 + enddo + ncarg=narg + if (use_otp) then + ncarg=ncarg-1 + endif + + if(ncarg.gt.ifile1) write(*,1999) nf,nd,nv +1999 format('nfiles:',i5,' ndecodes:',i5,' nverified:',i5) + +999 end program sfrx diff --git a/lib/superfox/twkfreq2.f90 b/lib/superfox/twkfreq2.f90 new file mode 100644 index 000000000..3cb26f930 --- /dev/null +++ b/lib/superfox/twkfreq2.f90 @@ -0,0 +1,19 @@ +subroutine twkfreq2(c3,c4,npts,fsample,fshift) + +! Adjust frequency of complex waveform + + complex c3(npts) + complex c4(npts) + complex w,wstep + data twopi/6.283185307/ + + w=1.0 + dphi=fshift*twopi/fsample + wstep=cmplx(cos(dphi),sin(dphi)) + do i=1,npts + w=w*wstep + c4(i)=w*c3(i) + enddo + + return +end subroutine twkfreq2