mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2025-03-23 12:38:53 -04:00
Remove the .../lib/superfox directory and its contents.
This commit is contained in:
Joe Taylor
parent
ce08913e3a
commit
ca18e58c77
@ -342,7 +342,6 @@ set (wsjt_FSRCS
|
||||
lib/wavhdr.f90
|
||||
lib/qra/q65/q65_encoding_modules.f90
|
||||
lib/ft8/ft8_a7.f90
|
||||
lib/superfox/sfox_mod.f90
|
||||
|
||||
# remaining non-module sources
|
||||
lib/addit.f90
|
||||
@ -428,9 +427,6 @@ set (wsjt_FSRCS
|
||||
lib/fspread_lorentz.f90
|
||||
lib/ft8/foxfilt.f90
|
||||
lib/ft8/foxgen.f90
|
||||
lib/superfox/foxgen2.f90
|
||||
lib/superfox/sfox_assemble.f90
|
||||
lib/superfox/sfox_wave.f90
|
||||
lib/ft8/foxgen_wrap.f90
|
||||
lib/freqcal.f90
|
||||
lib/ft8/ft8apset.f90
|
||||
@ -525,6 +521,7 @@ set (wsjt_FSRCS
|
||||
lib/sec0.f90
|
||||
lib/sec_midn.f90
|
||||
lib/setup65.f90
|
||||
lib/sfox_wave.f90
|
||||
lib/sh65.f90
|
||||
lib/sh65snr.f90
|
||||
lib/slasubs.f
|
||||
@ -590,13 +587,6 @@ set (wsjt_FSRCS
|
||||
lib/fst4/get_crc24.f90
|
||||
lib/fst4/fst4_baseline.f90
|
||||
lib/77bit/hash22calc.f90
|
||||
lib/superfox/sfox_gen.f90
|
||||
lib/superfox/sfox_sync.f90
|
||||
lib/superfox/sfox_demod.f90
|
||||
lib/superfox/sym_prob.f90
|
||||
lib/superfox/getpp3.f90
|
||||
lib/superfox/ftrsd3.f90
|
||||
lib/superfox/ran1.f90
|
||||
)
|
||||
|
||||
# temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit
|
||||
@ -637,10 +627,6 @@ set (wsjt_CSRCS
|
||||
lib/vit213.c
|
||||
lib/wisdom.c
|
||||
lib/wrapkarn.c
|
||||
lib/superfox/init_rs.c
|
||||
lib/superfox/encode_rs.c
|
||||
lib/superfox/decode_rs.c
|
||||
lib/superfox/rs_sf.c
|
||||
${ldpc_CSRCS}
|
||||
${qra_CSRCS}
|
||||
)
|
||||
@ -1169,12 +1155,6 @@ target_link_libraries (test_snr wsjt_fort)
|
||||
add_executable (q65sim lib/qra/q65/q65sim.f90)
|
||||
target_link_libraries (q65sim wsjt_fort wsjt_cxx)
|
||||
|
||||
#add_executable (rstest lib/superfox/rstest.f90)
|
||||
#target_link_libraries (rstest wsjt_fort wsjt_cxx)
|
||||
|
||||
add_executable (sfoxtest lib/superfox/sfoxtest.f90)
|
||||
target_link_libraries (sfoxtest wsjt_fort wsjt_cxx)
|
||||
|
||||
add_executable (q65code lib/qra/q65/q65code.f90)
|
||||
target_link_libraries (q65code wsjt_fort wsjt_cxx)
|
||||
|
||||
|
||||
@ -1,27 +0,0 @@
|
||||
CC = gcc
|
||||
FC = gfortran
|
||||
FFLAGS = -O2 -Wall -fbounds-check
|
||||
CFLAGS= -O9 -Wall
|
||||
|
||||
# Default rules
|
||||
%.o: %.c
|
||||
${CC} ${CFLAGS} -c $<
|
||||
%.o: %.f
|
||||
${FC} ${FFLAGS} -c $<
|
||||
%.o: %.F
|
||||
${FC} ${FFLAGS} -c $<
|
||||
%.o: %.f90
|
||||
${FC} ${FFLAGS} -c $<
|
||||
%.o: %.F90
|
||||
${FC} ${FFLAGS} -c $<
|
||||
|
||||
all: sfox_tx
|
||||
|
||||
OBJS1 = sfox_tx.o foxgen2.o sfox_assemble.o
|
||||
|
||||
sfox_tx: $(OBJS1)
|
||||
$(FC) -o sfox_tx $(OBJS1) libwsjt_fort.a
|
||||
|
||||
.PHONY : clean
|
||||
clean:
|
||||
-rm -f *.o *.exe
|
||||
@ -1,263 +0,0 @@
|
||||
/* Reed-Solomon decoder
|
||||
* Copyright 2002 Phil Karn, KA9Q
|
||||
* May be used under the terms of the GNU General Public License (GPL)
|
||||
*/
|
||||
|
||||
#ifdef DEBUG
|
||||
#include <stdio.h>
|
||||
#endif
|
||||
|
||||
#include <string.h>
|
||||
|
||||
#define NULL ((void *)0)
|
||||
#define min(a,b) ((a) < (b) ? (a) : (b))
|
||||
|
||||
#ifdef FIXED
|
||||
#include "fixed.h"
|
||||
#elif defined(BIGSYM)
|
||||
#include "int_sf.h"
|
||||
#else
|
||||
#include "char.h"
|
||||
#endif
|
||||
|
||||
int DECODE_RS(
|
||||
#ifdef FIXED
|
||||
DTYPE *data, int *eras_pos, int no_eras,int pad){
|
||||
#else
|
||||
void *p,DTYPE *data, int *eras_pos, int no_eras){
|
||||
struct rs *rs = (struct rs *)p;
|
||||
#endif
|
||||
int deg_lambda, el, deg_omega;
|
||||
int i, j, r,k;
|
||||
DTYPE u,q,tmp,num1,num2,den,discr_r;
|
||||
DTYPE lambda[NROOTS+1], s[NROOTS]; /* Err+Eras Locator poly
|
||||
* and syndrome poly */
|
||||
DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
|
||||
DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS];
|
||||
int syn_error, count;
|
||||
|
||||
#ifdef FIXED
|
||||
/* Check pad parameter for validity */
|
||||
if(pad < 0 || pad >= NN)
|
||||
return -1;
|
||||
#endif
|
||||
|
||||
/* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
|
||||
for(i=0;i<NROOTS;i++)
|
||||
s[i] = data[0];
|
||||
|
||||
for(j=1;j<NN-PAD;j++){
|
||||
for(i=0;i<NROOTS;i++){
|
||||
if(s[i] == 0){
|
||||
s[i] = data[j];
|
||||
} else {
|
||||
s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Convert syndromes to index form, checking for nonzero condition */
|
||||
syn_error = 0;
|
||||
for(i=0;i<NROOTS;i++){
|
||||
syn_error |= s[i];
|
||||
s[i] = INDEX_OF[s[i]];
|
||||
}
|
||||
|
||||
if (!syn_error) {
|
||||
/* if syndrome is zero, data[] is a codeword and there are no
|
||||
* errors to correct. So return data[] unmodified
|
||||
*/
|
||||
count = 0;
|
||||
goto finish;
|
||||
}
|
||||
memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
|
||||
lambda[0] = 1;
|
||||
|
||||
if (no_eras > 0) {
|
||||
/* Init lambda to be the erasure locator polynomial */
|
||||
lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
|
||||
for (i = 1; i < no_eras; i++) {
|
||||
u = MODNN(PRIM*(NN-1-eras_pos[i]));
|
||||
for (j = i+1; j > 0; j--) {
|
||||
tmp = INDEX_OF[lambda[j - 1]];
|
||||
if(tmp != A0)
|
||||
lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
|
||||
}
|
||||
}
|
||||
|
||||
#if DEBUG >= 1
|
||||
/* Test code that verifies the erasure locator polynomial just constructed
|
||||
Needed only for decoder debugging. */
|
||||
|
||||
/* find roots of the erasure location polynomial */
|
||||
for(i=1;i<=no_eras;i++)
|
||||
reg[i] = INDEX_OF[lambda[i]];
|
||||
|
||||
count = 0;
|
||||
for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
|
||||
q = 1;
|
||||
for (j = 1; j <= no_eras; j++)
|
||||
if (reg[j] != A0) {
|
||||
reg[j] = MODNN(reg[j] + j);
|
||||
q ^= ALPHA_TO[reg[j]];
|
||||
}
|
||||
if (q != 0)
|
||||
continue;
|
||||
/* store root and error location number indices */
|
||||
root[count] = i;
|
||||
loc[count] = k;
|
||||
count++;
|
||||
}
|
||||
if (count != no_eras) {
|
||||
printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
|
||||
count = -1;
|
||||
goto finish;
|
||||
}
|
||||
#if DEBUG >= 2
|
||||
printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
|
||||
for (i = 0; i < count; i++)
|
||||
printf("%d ", loc[i]);
|
||||
printf("\n");
|
||||
#endif
|
||||
#endif
|
||||
}
|
||||
for(i=0;i<NROOTS+1;i++)
|
||||
// printf("%d %d %d\n",i,lambda[i],INDEX_OF[lambda[i]]);
|
||||
b[i] = INDEX_OF[lambda[i]];
|
||||
|
||||
/*
|
||||
* Begin Berlekamp-Massey algorithm to determine error+erasure
|
||||
* locator polynomial
|
||||
*/
|
||||
r = no_eras;
|
||||
el = no_eras;
|
||||
while (++r <= NROOTS) { /* r is the step number */
|
||||
/* Compute discrepancy at the r-th step in poly-form */
|
||||
discr_r = 0;
|
||||
for (i = 0; i < r; i++){
|
||||
if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
|
||||
discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
|
||||
}
|
||||
}
|
||||
discr_r = INDEX_OF[discr_r]; /* Index form */
|
||||
if (discr_r == A0) {
|
||||
/* 2 lines below: B(x) <-- x*B(x) */
|
||||
memmove(&b[1],b,NROOTS*sizeof(b[0]));
|
||||
b[0] = A0;
|
||||
} else {
|
||||
/* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
|
||||
t[0] = lambda[0];
|
||||
for (i = 0 ; i < NROOTS; i++) {
|
||||
if(b[i] != A0)
|
||||
t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
|
||||
else
|
||||
t[i+1] = lambda[i+1];
|
||||
}
|
||||
if (2 * el <= r + no_eras - 1) {
|
||||
el = r + no_eras - el;
|
||||
/*
|
||||
* 2 lines below: B(x) <-- inv(discr_r) *
|
||||
* lambda(x)
|
||||
*/
|
||||
for (i = 0; i <= NROOTS; i++)
|
||||
b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
|
||||
} else {
|
||||
/* 2 lines below: B(x) <-- x*B(x) */
|
||||
memmove(&b[1],b,NROOTS*sizeof(b[0]));
|
||||
b[0] = A0;
|
||||
}
|
||||
memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
|
||||
}
|
||||
}
|
||||
|
||||
/* Convert lambda to index form and compute deg(lambda(x)) */
|
||||
deg_lambda = 0;
|
||||
for(i=0;i<NROOTS+1;i++){
|
||||
lambda[i] = INDEX_OF[lambda[i]];
|
||||
if(lambda[i] != A0)
|
||||
deg_lambda = i;
|
||||
}
|
||||
/* Find roots of the error+erasure locator polynomial by Chien search */
|
||||
memcpy(®[1],&lambda[1],NROOTS*sizeof(reg[0]));
|
||||
count = 0; /* Number of roots of lambda(x) */
|
||||
for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
|
||||
q = 1; /* lambda[0] is always 0 */
|
||||
for (j = deg_lambda; j > 0; j--){
|
||||
if (reg[j] != A0) {
|
||||
reg[j] = MODNN(reg[j] + j);
|
||||
q ^= ALPHA_TO[reg[j]];
|
||||
}
|
||||
}
|
||||
if (q != 0)
|
||||
continue; /* Not a root */
|
||||
/* store root (index-form) and error location number */
|
||||
#if DEBUG>=2
|
||||
printf("count %d root %d loc %d\n",count,i,k);
|
||||
#endif
|
||||
root[count] = i;
|
||||
loc[count] = k;
|
||||
/* If we've already found max possible roots,
|
||||
* abort the search to save time
|
||||
*/
|
||||
if(++count == deg_lambda)
|
||||
break;
|
||||
}
|
||||
if (deg_lambda != count) {
|
||||
/*
|
||||
* deg(lambda) unequal to number of roots => uncorrectable
|
||||
* error detected
|
||||
*/
|
||||
count = -1;
|
||||
goto finish;
|
||||
}
|
||||
/*
|
||||
* Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
|
||||
* x**NROOTS). in index form. Also find deg(omega).
|
||||
*/
|
||||
deg_omega = deg_lambda-1;
|
||||
for (i = 0; i <= deg_omega;i++){
|
||||
tmp = 0;
|
||||
for(j=i;j >= 0; j--){
|
||||
if ((s[i - j] != A0) && (lambda[j] != A0))
|
||||
tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
|
||||
}
|
||||
omega[i] = INDEX_OF[tmp];
|
||||
}
|
||||
|
||||
/*
|
||||
* Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
|
||||
* inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
|
||||
*/
|
||||
for (j = count-1; j >=0; j--) {
|
||||
num1 = 0;
|
||||
for (i = deg_omega; i >= 0; i--) {
|
||||
if (omega[i] != A0)
|
||||
num1 ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
|
||||
}
|
||||
num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
|
||||
den = 0;
|
||||
|
||||
/* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
|
||||
for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
|
||||
if(lambda[i+1] != A0)
|
||||
den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
|
||||
}
|
||||
#if DEBUG >= 1
|
||||
if (den == 0) {
|
||||
printf("\n ERROR: denominator = 0\n");
|
||||
count = -1;
|
||||
goto finish;
|
||||
}
|
||||
#endif
|
||||
/* Apply error to data */
|
||||
if (num1 != 0 && loc[j] >= PAD) {
|
||||
data[loc[j]-PAD] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
|
||||
}
|
||||
}
|
||||
finish:
|
||||
if(eras_pos != NULL){
|
||||
for(i=0;i<count;i++)
|
||||
eras_pos[i] = loc[i];
|
||||
}
|
||||
return count;
|
||||
}
|
||||
@ -1,52 +0,0 @@
|
||||
/* Reed-Solomon encoder
|
||||
* Copyright 2002, Phil Karn, KA9Q
|
||||
* May be used under the terms of the GNU General Public License (GPL)
|
||||
*/
|
||||
#include <string.h>
|
||||
|
||||
#ifdef FIXED
|
||||
#include "fixed.h"
|
||||
#elif defined(BIGSYM)
|
||||
#include "int_sf.h"
|
||||
#else
|
||||
#include "char.h"
|
||||
#endif
|
||||
|
||||
void ENCODE_RS(
|
||||
#ifdef FIXED
|
||||
DTYPE *data, DTYPE *bb,int pad){
|
||||
#else
|
||||
void *p,DTYPE *data, DTYPE *bb){
|
||||
struct rs *rs = (struct rs *)p;
|
||||
#endif
|
||||
int i, j;
|
||||
DTYPE feedback;
|
||||
|
||||
#ifdef FIXED
|
||||
/* Check pad parameter for validity */
|
||||
if(pad < 0 || pad >= NN)
|
||||
return;
|
||||
#endif
|
||||
|
||||
memset(bb,0,NROOTS*sizeof(DTYPE));
|
||||
|
||||
for(i=0;i<NN-NROOTS-PAD;i++){
|
||||
feedback = INDEX_OF[data[i] ^ bb[0]];
|
||||
if(feedback != A0){ /* feedback term is non-zero */
|
||||
#ifdef UNNORMALIZED
|
||||
/* This line is unnecessary when GENPOLY[NROOTS] is unity, as it must
|
||||
* always be for the polynomials constructed by init_rs()
|
||||
*/
|
||||
feedback = MODNN(NN - GENPOLY[NROOTS] + feedback);
|
||||
#endif
|
||||
for(j=1;j<NROOTS;j++)
|
||||
bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])];
|
||||
}
|
||||
/* Shift */
|
||||
memmove(&bb[0],&bb[1],sizeof(DTYPE)*(NROOTS-1));
|
||||
if(feedback != A0)
|
||||
bb[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
|
||||
else
|
||||
bb[NROOTS-1] = 0;
|
||||
}
|
||||
}
|
||||
@ -1,128 +0,0 @@
|
||||
subroutine foxgen2(nslots,cmsg,line)
|
||||
|
||||
! Called from foxgen() when it's time to encode a SuperFox message and
|
||||
! generate the waveform to be transmitted. We need to parse the old-style
|
||||
! Fox messages and extract the necessary pieces.
|
||||
|
||||
! use packjt77
|
||||
character*120 line
|
||||
character*40 cmsg(5) !Old-style Fox messages are here
|
||||
character*37 msg
|
||||
character*26 sfmsg
|
||||
character*13 mycall
|
||||
character*4 mygrid
|
||||
character*6 hiscall_1,hiscall_2
|
||||
character*4 rpt1,rpt2
|
||||
character*13 w(19)
|
||||
integer nw(19)
|
||||
integer ntype !Message type: 0 Free Text
|
||||
! 1 CQ MyCall MyGrid
|
||||
! 2 Call_1 MyCall RR73
|
||||
! 3 Call_1 MyCall rpt1
|
||||
! 4 Call_1 RR73; Call_2 <MyCall> rpt2
|
||||
|
||||
if(nslots.lt.1 .or. nslots.gt.5) return
|
||||
k=0
|
||||
do i=1,nslots
|
||||
hiscall_1=''
|
||||
hiscall_2=''
|
||||
mycall=''
|
||||
mygrid=''
|
||||
rpt1=''
|
||||
rpt2=''
|
||||
msg=cmsg(i)(1:37)
|
||||
call split77(msg,nwords,nw,w)
|
||||
ntype=0
|
||||
if(msg(1:3).eq.'CQ ') then
|
||||
ntype=1
|
||||
mycall=w(2)(1:12)
|
||||
mygrid=w(3)(1:4)
|
||||
else if(index(msg,';').gt.0) then
|
||||
ntype=4
|
||||
hiscall_1=w(1)(1:6)
|
||||
hiscall_2=w(3)(1:6)
|
||||
rpt1='RR73'
|
||||
rpt2=w(5)(1:4)
|
||||
mycall=w(4)(2:nw(4)-1)
|
||||
else if(index(msg,' RR73').gt.0) then
|
||||
ntype=2
|
||||
hiscall_1=w(1)(1:6)
|
||||
mycall=w(2)(1:12)
|
||||
rpt1='RR73'
|
||||
else if(nwords.eq.3 .and. nw(3).eq.3 .and. &
|
||||
(w(3)(1:1).eq.'-' .or. w(3)(1:1).eq.'+')) then
|
||||
ntype=3
|
||||
hiscall_1=w(1)(1:6)
|
||||
mycall=w(2)(1:12)
|
||||
rpt1=w(3)(1:4)
|
||||
endif
|
||||
! write(*,3001) ntype,cmsg(i),hiscall_1,rpt1,hiscall_2,rpt2, &
|
||||
! mycall(1:6),mygrid
|
||||
!3001 format(i1,2x,a37,1x,a6,1x,a4,1x,a6,1x,a4,1x,a6,1x,a4)
|
||||
|
||||
k=k+1
|
||||
if(ntype.le.3) call sfox_assemble(ntype,k,msg(1:26),mycall,mygrid,line)
|
||||
if(ntype.eq.4) then
|
||||
sfmsg=w(1)(1:nw(1))//' '//mycall(1:len(trim(mycall))+1)//'RR73'
|
||||
call sfox_assemble(2,k,sfmsg,mycall,mygrid,line)
|
||||
sfmsg=w(3)(1:nw(3))//' '//mycall(1:len(trim(mycall))+1)//w(5)(1:3)
|
||||
k=k+1
|
||||
call sfox_assemble(3,k,sfmsg,mycall,mygrid,line)
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
call sfox_assemble(ntype,11,msg(1:26),mycall,mygrid,line) !k=11 to finish up
|
||||
|
||||
return
|
||||
end subroutine foxgen2
|
||||
|
||||
subroutine split77(msg,nwords,nw,w)
|
||||
|
||||
! Convert msg to upper case; collapse multiple blanks; parse into words.
|
||||
|
||||
character*37 msg
|
||||
character*13 w(19)
|
||||
character*1 c,c0
|
||||
character*6 bcall_1
|
||||
logical ok1
|
||||
integer nw(19)
|
||||
|
||||
iz=len(trim(msg))
|
||||
j=0
|
||||
k=0
|
||||
n=0
|
||||
c0=' '
|
||||
w=' '
|
||||
do i=1,iz
|
||||
if(ichar(msg(i:i)).eq.0) msg(i:i)=' '
|
||||
c=msg(i:i) !Single character
|
||||
if(c.eq.' ' .and. c0.eq.' ') cycle !Skip leading/repeated blanks
|
||||
if(c.ne.' ' .and. c0.eq.' ') then
|
||||
k=k+1 !New word
|
||||
n=0
|
||||
endif
|
||||
j=j+1 !Index in msg
|
||||
n=n+1 !Index in word
|
||||
if(c.ge.'a' .and. c.le.'z') c=char(ichar(c)-32) !Force upper case
|
||||
msg(j:j)=c
|
||||
if(n.le.13) w(k)(n:n)=c !Copy character c into word
|
||||
c0=c
|
||||
enddo
|
||||
iz=j !Message length
|
||||
nwords=k !Number of words in msg
|
||||
if(nwords.le.0) go to 900
|
||||
do i=1,nwords
|
||||
nw(i)=len(trim(w(i)))
|
||||
enddo
|
||||
msg(iz+1:)=' '
|
||||
if(nwords.lt.3) go to 900
|
||||
call chkcall(w(3),bcall_1,ok1)
|
||||
if(ok1 .and. w(1)(1:3).eq.'CQ ') then
|
||||
w(1)='CQ_'//w(2)(1:10) !Make "CQ " into "CQ_"
|
||||
w(2:12)=w(3:13) !Move all remaining words down by one
|
||||
nwords=nwords-1
|
||||
endif
|
||||
|
||||
900 return
|
||||
end subroutine split77
|
||||
@ -1,217 +0,0 @@
|
||||
/*
|
||||
ftrsd3.c
|
||||
|
||||
A soft-decision decoder for Reed-Solomon codes.
|
||||
|
||||
This decoding scheme is built around Phil Karn's Berlekamp-Massey
|
||||
errors and erasures decoder. The approach is inspired by a number of
|
||||
publications, including the stochastic Chase decoder described
|
||||
in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
|
||||
IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
|
||||
"Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
|
||||
and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.
|
||||
|
||||
Steve Franke K9AN and Joe Taylor K1JT
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <unistd.h>
|
||||
#include <time.h>
|
||||
#include <string.h>
|
||||
#include "rs2.h"
|
||||
|
||||
static void *rs;
|
||||
void getpp3_(int workdat[], float *pp);
|
||||
|
||||
void ftrsd3_(int rxdat[], int rxprob[], int rxdat2[], int rxprob2[],
|
||||
int* ntrials0, int correct[], int param[], int ntry[])
|
||||
{
|
||||
// int rxdat[127], rxprob[127], rxdat2[127], rxprob2[127];
|
||||
int workdat[127];
|
||||
int indexes[127];
|
||||
int era_pos[79];
|
||||
int i, j, numera, nerr, nn=127;
|
||||
int ntrials = *ntrials0;
|
||||
int nhard=0,nhard_min=32768,nsoft=0,nsoft_min=32768;
|
||||
int ntotal=0,ntotal_min=32768,ncandidates;
|
||||
int nera_best=0;
|
||||
float pp,pp1,pp2;
|
||||
static unsigned int nseed;
|
||||
|
||||
// Power-percentage symbol metrics - composite gnnf/hf
|
||||
int perr[8][8] = {
|
||||
{ 4, 9, 11, 13, 14, 14, 15, 15},
|
||||
{ 2, 20, 20, 30, 40, 50, 50, 50},
|
||||
{ 7, 24, 27, 40, 50, 50, 50, 50},
|
||||
{13, 25, 35, 46, 52, 70, 50, 50},
|
||||
{17, 30, 42, 54, 55, 64, 71, 70},
|
||||
{25, 39, 48, 57, 64, 66, 77, 77},
|
||||
{32, 45, 54, 63, 66, 75, 78, 83},
|
||||
{51, 58, 57, 66, 72, 77, 82, 86}};
|
||||
|
||||
printf("A\n");
|
||||
// Initialize the KA9Q Reed-Solomon encoder/decoder
|
||||
unsigned int symsize=7, gfpoly=0x89, fcr=3, prim=1, nroots=79;
|
||||
rs=init_rs_int(symsize, gfpoly, fcr, prim, nroots, 0);
|
||||
printf("B\n");
|
||||
|
||||
// Sort rxprob to find indexes of the least reliable symbols
|
||||
int k, pass, tmp, nsym=127;
|
||||
int probs[127];
|
||||
for (i=0; i<127; i++) {
|
||||
indexes[i]=i;
|
||||
probs[i]=rxprob[i];
|
||||
}
|
||||
printf("C\n");
|
||||
|
||||
for (pass = 1; pass <= nsym-1; pass++) {
|
||||
for (k = 0; k < nsym - pass; k++) {
|
||||
if( probs[k] < probs[k+1] ) {
|
||||
tmp = probs[k];
|
||||
probs[k] = probs[k+1];
|
||||
probs[k+1] = tmp;
|
||||
tmp = indexes[k];
|
||||
indexes[k] = indexes[k+1];
|
||||
indexes[k+1] = tmp;
|
||||
}
|
||||
}
|
||||
}
|
||||
printf("D\n");
|
||||
|
||||
// See if we can decode using BM HDD, and calculate the syndrome vector.
|
||||
memset(era_pos,0,79*sizeof(int));
|
||||
numera=0;
|
||||
memcpy(workdat,rxdat,127*sizeof(int));
|
||||
nerr=decode_rs_int(rs,workdat,era_pos,numera,1);
|
||||
printf("E\n");
|
||||
|
||||
if( nerr >= 0 ) {
|
||||
// Hard-decision decoding succeeded. Save codeword and some parameters.
|
||||
nhard=0;
|
||||
for (i=0; i<127; i++) {
|
||||
if( workdat[i] != rxdat[i] ) nhard=nhard+1;
|
||||
}
|
||||
memcpy(correct,workdat,127*sizeof(int));
|
||||
param[0]=0;
|
||||
param[1]=nhard;
|
||||
param[2]=0;
|
||||
param[3]=0;
|
||||
param[4]=0;
|
||||
param[5]=0;
|
||||
param[7]=1000*1000;
|
||||
ntry[0]=0;
|
||||
return;
|
||||
}
|
||||
|
||||
if(nerr==-1) {
|
||||
printf("nerr:', %d\n",nerr);
|
||||
return;
|
||||
}
|
||||
|
||||
/*
|
||||
Hard-decision decoding failed. Try the FT soft-decision method.
|
||||
Generate random erasure-locator vectors and see if any of them
|
||||
decode. This will generate a list of "candidate" codewords. The
|
||||
soft distance between each candidate codeword and the received
|
||||
word is estimated by finding the largest (pp1) and second-largest
|
||||
(pp2) outputs from a synchronized filter-bank operating on the
|
||||
symbol spectra, and using these to decide which candidate
|
||||
codeword is "best".
|
||||
*/
|
||||
|
||||
nseed=1; //Seed for random numbers
|
||||
float ratio;
|
||||
int thresh, nsum;
|
||||
int thresh0[127];
|
||||
ncandidates=0;
|
||||
nsum=0;
|
||||
int ii,jj;
|
||||
for (i=0; i<nn; i++) {
|
||||
nsum=nsum+rxprob[i];
|
||||
j = indexes[126-i];
|
||||
ratio = (float)rxprob2[j]/((float)rxprob[j]+0.01);
|
||||
ii = 7.999*ratio;
|
||||
jj = (126-i)/8;
|
||||
thresh0[i] = 1.3*perr[ii][jj];
|
||||
}
|
||||
if(nsum<=0) return;
|
||||
|
||||
pp1=0.0;
|
||||
pp2=0.0;
|
||||
for (k=1; k<=ntrials; k++) {
|
||||
memset(era_pos,0,79*sizeof(int));
|
||||
memcpy(workdat,rxdat,127*sizeof(int));
|
||||
|
||||
/*
|
||||
Mark a subset of the symbols as erasures.
|
||||
Run through the ranked symbols, starting with the worst, i=0.
|
||||
NB: j is the symbol-vector index of the symbol with rank i.
|
||||
*/
|
||||
numera=0;
|
||||
for (i=0; i<nn; i++) {
|
||||
j = indexes[126-i];
|
||||
thresh=thresh0[i];
|
||||
long int ir;
|
||||
|
||||
// Generate a random number ir, 0 <= ir < 100 (see POSIX.1-2001 example).
|
||||
nseed = nseed * 1103515245 + 12345;
|
||||
ir = (unsigned)(nseed/65536) % 32768;
|
||||
ir = (100*ir)/32768;
|
||||
|
||||
if((ir < thresh ) && numera < 79) {
|
||||
era_pos[numera]=j;
|
||||
numera=numera+1;
|
||||
}
|
||||
}
|
||||
|
||||
nerr=decode_rs_int(rs,workdat,era_pos,numera,0);
|
||||
if( nerr >= 0 ) {
|
||||
// We have a candidate codeword. Find its hard and soft distance from
|
||||
// the received word. Also find pp1 and pp2 from the full array
|
||||
// s3(64,127) of synchronized symbol spectra.
|
||||
ncandidates=ncandidates+1;
|
||||
nhard=0;
|
||||
nsoft=0;
|
||||
for (i=0; i<127; i++) {
|
||||
if(workdat[i] != rxdat[i]) {
|
||||
nhard=nhard+1;
|
||||
if(workdat[i] != rxdat2[i]) {
|
||||
nsoft=nsoft+rxprob[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
nsoft=127*nsoft/nsum;
|
||||
ntotal=nsoft+nhard;
|
||||
|
||||
pp=0.;
|
||||
// getpp3_(workdat,&pp);
|
||||
if(pp>pp1) {
|
||||
pp2=pp1;
|
||||
pp1=pp;
|
||||
nsoft_min=nsoft;
|
||||
nhard_min=nhard;
|
||||
ntotal_min=ntotal;
|
||||
memcpy(correct,workdat,127*sizeof(int));
|
||||
nera_best=numera;
|
||||
ntry[0]=k;
|
||||
} else {
|
||||
if(pp>pp2 && pp!=pp1) pp2=pp;
|
||||
}
|
||||
if(nhard_min <= 41 && ntotal_min <= 71) break;
|
||||
}
|
||||
if(k == ntrials) ntry[0]=k;
|
||||
}
|
||||
|
||||
param[0]=ncandidates;
|
||||
param[1]=nhard_min;
|
||||
param[2]=nsoft_min;
|
||||
param[3]=nera_best;
|
||||
param[4]= pp1 > 0 ? 1000.0*pp2/pp1 : 1000.0;
|
||||
param[5]=ntotal_min;
|
||||
param[6]=ntry[0];
|
||||
param[7]=1000.0*pp2;
|
||||
param[8]=1000.0*pp1;
|
||||
if(param[0]==0) param[2]=-1;
|
||||
return;
|
||||
}
|
||||
@ -1,192 +0,0 @@
|
||||
subroutine ftrsd3(s3,chansym0,rxdat,rxprob,rxdat2,rxprob2,ntrials0, &
|
||||
correct,param,ntry)
|
||||
|
||||
! Soft-decision decoder for Reed-Solomon codes.
|
||||
|
||||
! This decoding scheme is built around Phil Karn's Berlekamp-Massey
|
||||
! errors and erasures decoder. The approach is inspired by a number of
|
||||
! publications, including the stochastic Chase decoder described
|
||||
! in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
|
||||
! IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
|
||||
! "Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
|
||||
! and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.
|
||||
|
||||
! Steve Franke K9AN and Joe Taylor K1JT
|
||||
|
||||
use sfox_mod
|
||||
|
||||
real s3(0:NQ-1,0:NN-1) !Symbol spectra
|
||||
integer chansym0(0:NN-1) !Transmitted codeword
|
||||
integer rxdat(0:NN-1) !Hard-decision symbol values
|
||||
integer rxprob(0:NN-1) !Probabilities that rxdat values are correct
|
||||
integer rxdat2(0:NN-1) !Second most probable symbol values
|
||||
integer rxprob2(0:NN-1) !Probabilities that rxdat2 values are correct
|
||||
integer workdat(0:NN-1) !Work array
|
||||
integer correct(0:NN-1) !Corrected codeword
|
||||
integer indexes(0:NN-1) !For sorting probabilities
|
||||
integer probs(0:NN-1) !Temp array for sorting probabilities
|
||||
integer thresh0(0:NN-1) !Temp array for thresholds
|
||||
integer era_pos(0:NN-KK-1) !Index values for erasures
|
||||
integer param(0:8)
|
||||
integer*8 nseed,ir !No unsigned int in Fortran
|
||||
integer pass,tmp,thresh
|
||||
|
||||
integer perr(0:7,0:7)
|
||||
data perr/ 4, 9,11,13,14,14,15,15, &
|
||||
2,20,20,30,40,50,50,50, &
|
||||
7,24,27,40,50,50,50,50, &
|
||||
13,25,35,46,52,70,50,50, &
|
||||
17,30,42,54,55,64,71,70, &
|
||||
25,39,48,57,64,66,77,77, &
|
||||
32,45,54,63,66,75,78,83, &
|
||||
51,58,57,66,72,77,82,86/
|
||||
|
||||
ntrials=ntrials0
|
||||
nhard=0
|
||||
nhard_min=32768
|
||||
nsoft=0
|
||||
nsoft_min=32768
|
||||
ntotal=0
|
||||
ntotal_min=32768
|
||||
nera_best=0
|
||||
nsym=nn
|
||||
|
||||
do i=0,NN-1
|
||||
indexes(i)=i
|
||||
probs(i)=rxprob(i)
|
||||
enddo
|
||||
|
||||
do pass=1,nsym-1
|
||||
do k=0,nsym-pass-1
|
||||
if(probs(k).lt.probs(k+1)) then
|
||||
tmp=probs(k)
|
||||
probs(k)=probs(k+1)
|
||||
probs(k+1)=tmp
|
||||
tmp=indexes(k)
|
||||
indexes(k)=indexes(k+1)
|
||||
indexes(k+1)=tmp
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
correct=-1
|
||||
era_pos=0
|
||||
numera=0
|
||||
workdat=rxdat
|
||||
call rs_decode_sf(workdat,era_pos,numera,nerr) !Call the decoder
|
||||
nerr=-1
|
||||
|
||||
if(nerr.ge.0) then
|
||||
! Hard-decision decoding succeeded. Save codeword and some parameters.
|
||||
nhard=count(workdat.ne.rxdat)
|
||||
correct=workdat
|
||||
param(0)=0
|
||||
param(1)=nhard
|
||||
param(2)=0
|
||||
param(3)=0
|
||||
param(4)=0
|
||||
param(5)=0
|
||||
param(7)=1000*1000 !???
|
||||
ntry=0
|
||||
go to 900
|
||||
endif
|
||||
|
||||
! Hard-decision decoding failed. Try the FT soft-decision method.
|
||||
! Generate random erasure-locator vectors and see if any of them
|
||||
! decode. This will generate a list of "candidate" codewords. The
|
||||
! soft distance between each candidate codeword and the received
|
||||
! word is estimated by finding the largest (pp1) and second-largest
|
||||
! (pp2) outputs from a synchronized filter-bank operating on the
|
||||
! symbol spectra, and using these to decide which candidate
|
||||
! codeword is "best".
|
||||
|
||||
nseed=1 !Seed for random numbers
|
||||
ncandidates=0
|
||||
nsum=0
|
||||
do i=0,NN-1
|
||||
nsum=nsum+rxprob(i)
|
||||
j=indexes(NN-1-i)
|
||||
ratio=float(rxprob2(j))/(float(rxprob(j))+0.01)
|
||||
ii=7.999*ratio
|
||||
jj=int((7.999/NN)*(NN-1-i))
|
||||
thresh0(i)=0.90*perr(jj,ii)
|
||||
enddo
|
||||
if(nsum.le.0) return
|
||||
|
||||
pp1=0.
|
||||
pp2=0.
|
||||
do k=1,ntrials
|
||||
era_pos=0
|
||||
workdat=rxdat
|
||||
|
||||
! Mark a subset of the symbols as erasures.
|
||||
! Run through the ranked symbols, starting with the worst, i=0.
|
||||
! NB: j is the symbol-vector index of the symbol with rank i.
|
||||
|
||||
ncaught=0
|
||||
numera=0
|
||||
do i=0,NN-1
|
||||
j=indexes(NN-1-i)
|
||||
thresh=thresh0(i)
|
||||
! Generate a random number ir, 0 <= ir <= 100 (see POSIX.1-2001 example).
|
||||
ir=100.0*ran1(nseed)
|
||||
if((ir.lt.thresh) .and. numera.lt. 0.69*(NN-KK)) then
|
||||
era_pos(numera)=j
|
||||
numera=numera+1
|
||||
if(rxdat(j).ne.chansym0(j)) then
|
||||
ncaught=ncaught+1
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
call rs_decode_sf(workdat,era_pos,numera,nerr) !Call the decoder
|
||||
if( nerr.ge.0) then
|
||||
! We have a candidate codeword. Find its hard and soft distance from
|
||||
! the received word. Also find pp1 and pp2 from the full array
|
||||
! s3(NQ,NN) of synchronized symbol spectra.
|
||||
ncandidates=ncandidates+1
|
||||
nhard=0
|
||||
nsoft=0
|
||||
do i=0,NN-1
|
||||
if(workdat(i).ne. rxdat(i)) then
|
||||
nhard=nhard+1;
|
||||
if(workdat(i) .ne. rxdat2(i)) nsoft=nsoft+rxprob(i)
|
||||
endif
|
||||
enddo
|
||||
nsoft=NN*nsoft/nsum
|
||||
ntotal=nsoft+nhard
|
||||
|
||||
pp=0.
|
||||
call getpp3(s3,workdat,pp)
|
||||
! write(*,5001) ncandidates,nhard,nsoft,ntotal,pp,pp1,pp2
|
||||
!5001 format(4i8,3f7.3)
|
||||
if(pp.gt.pp1) then
|
||||
pp2=pp1
|
||||
pp1=pp
|
||||
nsoft_min=nsoft
|
||||
nhard_min=nhard
|
||||
ntotal_min=ntotal
|
||||
correct=workdat
|
||||
nera_best=numera
|
||||
ntry=k
|
||||
else
|
||||
if(pp.gt.pp2 .and. pp.ne.pp1) pp2=pp
|
||||
endif
|
||||
if(nhard_min.le.60 .and. ntotal_min.le.90) exit !### Needs tuning
|
||||
endif
|
||||
if(k.eq.ntrials) ntry=k
|
||||
enddo
|
||||
|
||||
param(0)=ncandidates
|
||||
param(1)=nhard_min
|
||||
param(2)=nsoft_min
|
||||
param(3)=nera_best
|
||||
param(4)=1000
|
||||
if(pp1.gt.0.0) param(4)=1000.0*pp2/pp1
|
||||
param(5)=ntotal_min
|
||||
param(6)=ntry
|
||||
param(7)=1000.0*pp2
|
||||
param(8)=1000.0*pp1
|
||||
if(param(0).eq.0) param(2)=-1
|
||||
|
||||
900 return
|
||||
end subroutine ftrsd3
|
||||
@ -1,25 +0,0 @@
|
||||
subroutine get_crc14(mc,len,ncrc)
|
||||
!
|
||||
! 1. To calculate 14-bit CRC, mc(1:len-14) is the message and mc(len-13:len) are zero.
|
||||
! 2. To check a received CRC, mc(1:len is the received message plus CRC.
|
||||
! ncrc will be zero if the received message/CRC are consistent
|
||||
!
|
||||
character c14*14
|
||||
integer*1 mc(len)
|
||||
integer*1 r(15),p(15)
|
||||
integer ncrc
|
||||
! polynomial for 14-bit CRC 0x6757
|
||||
data p/1,1,0,0,1,1,1,0,1,0,1,0,1,1,1/
|
||||
|
||||
! divide by polynomial
|
||||
r=mc(1:15)
|
||||
do i=0,len-15
|
||||
r(15)=mc(i+15)
|
||||
r=mod(r+r(1)*p,2)
|
||||
r=cshift(r,1)
|
||||
enddo
|
||||
|
||||
write(c14,'(14b1)') r(1:14)
|
||||
read(c14,'(b14.14)') ncrc
|
||||
|
||||
end subroutine get_crc14
|
||||
@ -1,22 +0,0 @@
|
||||
subroutine getpp3(s3,workdat,p)
|
||||
|
||||
use sfox_mod
|
||||
real s3(NQ,NN)
|
||||
integer workdat(NN)
|
||||
integer a(NN)
|
||||
|
||||
! a(1:NN)=workdat(NN:1:-1)
|
||||
a=workdat
|
||||
|
||||
psum=0.
|
||||
do j=1,NN
|
||||
i=a(j)+1
|
||||
x=s3(i,j)
|
||||
s3(i,j)=0.
|
||||
psum=psum + x
|
||||
s3(i,j)=x
|
||||
enddo
|
||||
p=psum/NN
|
||||
|
||||
return
|
||||
end subroutine getpp3
|
||||
@ -1,126 +0,0 @@
|
||||
/* Initialize a RS codec
|
||||
*
|
||||
* Copyright 2002 Phil Karn, KA9Q
|
||||
* May be used under the terms of the GNU General Public License (GPL)
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
|
||||
#ifdef CCSDS
|
||||
#include "ccsds.h"
|
||||
#elif defined(BIGSYM)
|
||||
#include "int_sf.h"
|
||||
#else
|
||||
#include "char.h"
|
||||
#endif
|
||||
|
||||
#define NULL ((void *)0)
|
||||
|
||||
void FREE_RS(void *p){
|
||||
struct rs *rs = (struct rs *)p;
|
||||
|
||||
free(rs->alpha_to);
|
||||
free(rs->index_of);
|
||||
free(rs->genpoly);
|
||||
free(rs);
|
||||
}
|
||||
|
||||
/* Initialize a Reed-Solomon codec
|
||||
* symsize = symbol size, bits (1-8)
|
||||
* gfpoly = Field generator polynomial coefficients
|
||||
* fcr = first root of RS code generator polynomial, index form
|
||||
* prim = primitive element to generate polynomial roots
|
||||
* nroots = RS code generator polynomial degree (number of roots)
|
||||
* pad = padding bytes at front of shortened block
|
||||
*/
|
||||
void *INIT_RS(int symsize,int gfpoly,int fcr,int prim,
|
||||
int nroots,int pad){
|
||||
struct rs *rs;
|
||||
int i, j, sr,root,iprim;
|
||||
|
||||
/* Check parameter ranges */
|
||||
if(symsize < 0 || symsize > 8*sizeof(DTYPE))
|
||||
return NULL; /* Need version with ints rather than chars */
|
||||
|
||||
if(fcr < 0 || fcr >= (1<<symsize))
|
||||
return NULL;
|
||||
if(prim <= 0 || prim >= (1<<symsize))
|
||||
return NULL;
|
||||
if(nroots < 0 || nroots >= (1<<symsize))
|
||||
return NULL; /* Can't have more roots than symbol values! */
|
||||
if(pad < 0 || pad >= ((1<<symsize) -1 - nroots))
|
||||
return NULL; /* Too much padding */
|
||||
|
||||
rs = (struct rs *)calloc(1,sizeof(struct rs));
|
||||
rs->mm = symsize;
|
||||
rs->nn = (1<<symsize)-1;
|
||||
rs->pad = pad;
|
||||
|
||||
rs->alpha_to = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
|
||||
if(rs->alpha_to == NULL){
|
||||
free(rs);
|
||||
return NULL;
|
||||
}
|
||||
rs->index_of = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
|
||||
if(rs->index_of == NULL){
|
||||
free(rs->alpha_to);
|
||||
free(rs);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* Generate Galois field lookup tables */
|
||||
rs->index_of[0] = A0; /* log(zero) = -inf */
|
||||
rs->alpha_to[A0] = 0; /* alpha**-inf = 0 */
|
||||
sr = 1;
|
||||
for(i=0;i<rs->nn;i++){
|
||||
rs->index_of[sr] = i;
|
||||
rs->alpha_to[i] = sr;
|
||||
sr <<= 1;
|
||||
if(sr & (1<<symsize))
|
||||
sr ^= gfpoly;
|
||||
sr &= rs->nn;
|
||||
}
|
||||
if(sr != 1){
|
||||
/* field generator polynomial is not primitive! */
|
||||
free(rs->alpha_to);
|
||||
free(rs->index_of);
|
||||
free(rs);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* Form RS code generator polynomial from its roots */
|
||||
rs->genpoly = (DTYPE *)malloc(sizeof(DTYPE)*(nroots+1));
|
||||
if(rs->genpoly == NULL){
|
||||
free(rs->alpha_to);
|
||||
free(rs->index_of);
|
||||
free(rs);
|
||||
return NULL;
|
||||
}
|
||||
rs->fcr = fcr;
|
||||
rs->prim = prim;
|
||||
rs->nroots = nroots;
|
||||
|
||||
/* Find prim-th root of 1, used in decoding */
|
||||
for(iprim=1;(iprim % prim) != 0;iprim += rs->nn)
|
||||
;
|
||||
rs->iprim = iprim / prim;
|
||||
|
||||
rs->genpoly[0] = 1;
|
||||
for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) {
|
||||
rs->genpoly[i+1] = 1;
|
||||
|
||||
/* Multiply rs->genpoly[] by @**(root + x) */
|
||||
for (j = i; j > 0; j--){
|
||||
if (rs->genpoly[j] != 0)
|
||||
rs->genpoly[j] = rs->genpoly[j-1] ^ rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[j]] + root)];
|
||||
else
|
||||
rs->genpoly[j] = rs->genpoly[j-1];
|
||||
}
|
||||
/* rs->genpoly[0] can never be zero */
|
||||
rs->genpoly[0] = rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[0]] + root)];
|
||||
}
|
||||
/* convert rs->genpoly[] to index form for quicker encoding */
|
||||
for (i = 0; i <= nroots; i++)
|
||||
rs->genpoly[i] = rs->index_of[rs->genpoly[i]];
|
||||
|
||||
return rs;
|
||||
}
|
||||
@ -1,56 +0,0 @@
|
||||
/* Include file to configure the RS codec for integer symbols
|
||||
*
|
||||
* Copyright 2002, Phil Karn, KA9Q
|
||||
* May be used under the terms of the GNU General Public License (GPL)
|
||||
*/
|
||||
#define DTYPE int
|
||||
|
||||
/* Reed-Solomon codec control block */
|
||||
struct rs {
|
||||
int mm; /* Bits per symbol */
|
||||
int nn; /* Symbols per block (= (1<<mm)-1) */
|
||||
DTYPE *alpha_to; /* log lookup table */
|
||||
DTYPE *index_of; /* Antilog lookup table */
|
||||
DTYPE *genpoly; /* Generator polynomial */
|
||||
int nroots; /* Number of generator roots = number of parity symbols */
|
||||
int fcr; /* First consecutive root, index form */
|
||||
int prim; /* Primitive element, index form */
|
||||
int iprim; /* prim-th root of 1, index form */
|
||||
int pad; /* Padding bytes in shortened block */
|
||||
};
|
||||
|
||||
static inline int modnn(struct rs *rs,int x){
|
||||
while (x >= rs->nn) {
|
||||
x -= rs->nn;
|
||||
x = (x >> rs->mm) + (x & rs->nn);
|
||||
}
|
||||
return x;
|
||||
}
|
||||
#define MODNN(x) modnn(rs,x)
|
||||
|
||||
#define MM (rs->mm)
|
||||
#define NN (rs->nn)
|
||||
#define ALPHA_TO (rs->alpha_to)
|
||||
#define INDEX_OF (rs->index_of)
|
||||
#define GENPOLY (rs->genpoly)
|
||||
#define NROOTS (rs->nroots)
|
||||
#define FCR (rs->fcr)
|
||||
#define PRIM (rs->prim)
|
||||
#define IPRIM (rs->iprim)
|
||||
#define PAD (rs->pad)
|
||||
#define A0 (NN)
|
||||
|
||||
#define ENCODE_RS encode_rs_sf
|
||||
#define DECODE_RS decode_rs_sf
|
||||
#define INIT_RS init_rs_sf
|
||||
#define FREE_RS free_rs_sf
|
||||
|
||||
void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity);
|
||||
int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras);
|
||||
void *INIT_RS(int symsize,int gfpoly,int fcr,
|
||||
int prim,int nroots,int pad);
|
||||
void FREE_RS(void *p);
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,39 +0,0 @@
|
||||
subroutine plotspec(dat)
|
||||
|
||||
use sfox_mod
|
||||
real dat(NMAX)
|
||||
real s(0:NSPS/2)
|
||||
complex c(0:NSPS-1)
|
||||
integer ipk(1)
|
||||
|
||||
nblks=NZ/NSPS
|
||||
s=0.
|
||||
fac=1.0/NSPS
|
||||
do j=1,nblks
|
||||
ib=j*NSPS
|
||||
ia=ib-NSPS+1
|
||||
c=fac*dat(ia:ib)
|
||||
call four2a(c,NSPS,1,-1,1)
|
||||
do i=0,NSPS/2
|
||||
s(i)=s(i) + real(c(i))**2 + aimag(c(i))**2
|
||||
enddo
|
||||
enddo
|
||||
|
||||
df=12000.0/NSPS
|
||||
ipk=maxloc(s)
|
||||
f0=df*(ipk(1)-1)
|
||||
p_sig_plus_noise=maxval(s)
|
||||
p_noise=0.
|
||||
do i=0,NSPS/2
|
||||
f=i*df
|
||||
if(f.le.2500+df .and. abs(f-f0).gt.0.5*df) p_noise=p_noise + s(i)
|
||||
write(40,1000) f,s(i)
|
||||
1000 format(2f10.3)
|
||||
enddo
|
||||
p_sig=p_sig_plus_noise - p_noise*df/2500.0
|
||||
snr=p_sig/p_noise
|
||||
snrdb=db(snr)
|
||||
write(*,1100) snrdb
|
||||
1100 format('Measured SNR:',f7.2)
|
||||
|
||||
end subroutine plotspec
|
||||
@ -1,28 +0,0 @@
|
||||
FUNCTION ran1(idum)
|
||||
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
|
||||
REAL ran1,AM,EPS,RNMX
|
||||
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, &
|
||||
NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
|
||||
INTEGER j,k,iv(NTAB),iy
|
||||
SAVE iv,iy
|
||||
DATA iv /NTAB*0/, iy /0/
|
||||
if (idum.le.0.or.iy.eq.0) then
|
||||
idum=max(-idum,1)
|
||||
do j=NTAB+8,1,-1
|
||||
k=idum/IQ
|
||||
idum=IA*(idum-k*IQ)-IR*k
|
||||
if (idum.lt.0) idum=idum+IM
|
||||
if (j.le.NTAB) iv(j)=idum
|
||||
enddo
|
||||
iy=iv(1)
|
||||
endif
|
||||
k=idum/IQ
|
||||
idum=IA*(idum-k*IQ)-IR*k
|
||||
if (idum.lt.0) idum=idum+IM
|
||||
j=1+iy/NDIV
|
||||
iy=iv(j)
|
||||
iv(j)=idum
|
||||
ran1=min(AM*iy,RNMX)
|
||||
|
||||
return
|
||||
END FUNCTION ran1
|
||||
@ -1,33 +0,0 @@
|
||||
#include <stdio.h>
|
||||
#include "rs_sf.h"
|
||||
|
||||
static void *rs_sf;
|
||||
static int first=1;
|
||||
static int nn,kk,nroots,npad;
|
||||
|
||||
void rs_init_sf_(int *mm, int *nq, int *nn0, int *kk0, int *nfz)
|
||||
// Initialize the RS decoder.
|
||||
{
|
||||
// Save parameters nn, kk, nroots, npad for global access
|
||||
nn=*nn0;
|
||||
kk=*kk0;
|
||||
nroots=nn-kk;
|
||||
npad=*nq-1-nn;
|
||||
|
||||
int gfpoly=0x43; //For *mm=6
|
||||
if(*mm==7) gfpoly=0x89;
|
||||
if(*mm==8) gfpoly=0x11d;
|
||||
rs_sf=init_rs_sf(*mm,gfpoly,*nfz,1,nroots,npad);
|
||||
first=0;
|
||||
}
|
||||
|
||||
void rs_encode_sf_(int *msg, int *parsym)
|
||||
// Encode information symbols msg[KK], producing parity symbols parsym[nroots].
|
||||
{
|
||||
encode_rs_sf(rs_sf,msg,parsym); //Compute the parity symbols
|
||||
}
|
||||
|
||||
void rs_decode_sf_(int *recd, int *era_pos, int *numera, int *nerr)
|
||||
{
|
||||
*nerr=decode_rs_sf(rs_sf,recd,era_pos,*numera);
|
||||
}
|
||||
@ -1,35 +0,0 @@
|
||||
/* User include file for the Reed-Solomon codec
|
||||
* Copyright 2002, Phil Karn KA9Q
|
||||
* May be used under the terms of the GNU General Public License (GPL)
|
||||
*/
|
||||
|
||||
/* General purpose RS codec, 8-bit symbols */
|
||||
void encode_rs_char(void *rs,unsigned char *data,unsigned char *parity);
|
||||
int decode_rs_char(void *rs,unsigned char *data,int *eras_pos,
|
||||
int no_eras);
|
||||
void *init_rs_char(int symsize,int gfpoly,
|
||||
int fcr,int prim,int nroots,
|
||||
int pad);
|
||||
void free_rs_char(void *rs);
|
||||
|
||||
/* General purpose RS codec, integer symbols */
|
||||
void encode_rs_sf(void *rs,int *data,int *parity);
|
||||
int decode_rs_sf(void *rs,int *data,int *eras_pos,int no_eras);
|
||||
void *init_rs_sf(int symsize,int gfpoly,int fcr,
|
||||
int prim,int nroots,int pad);
|
||||
void free_rs_sf(void *rs);
|
||||
|
||||
/* CCSDS standard (255,223) RS codec with conventional (*not* dual-basis)
|
||||
* symbol representation
|
||||
*/
|
||||
void encode_rs_8(unsigned char *data,unsigned char *parity,int pad);
|
||||
int decode_rs_8(unsigned char *data,int *eras_pos,int no_eras,int pad);
|
||||
|
||||
/* CCSDS standard (255,223) RS codec with dual-basis symbol representation */
|
||||
void encode_rs_ccsds(unsigned char *data,unsigned char *parity,int pad);
|
||||
int decode_rs_ccsds(unsigned char *data,int *eras_pos,int no_eras,int pad);
|
||||
|
||||
/* Tables to map from conventional->dual (Taltab) and
|
||||
* dual->conventional (Tal1tab) bases
|
||||
*/
|
||||
extern unsigned char Taltab[],Tal1tab[];
|
||||
@ -1,71 +0,0 @@
|
||||
program rst8
|
||||
|
||||
character arg*8
|
||||
integer*1 dat0(223) !Generated data
|
||||
integer*1 parsym(32) !Parity symbols
|
||||
integer*1 cword0(255) !Generated codeword
|
||||
integer*1 cword(255) !Rcvd codeword with errors; will be corrected in place
|
||||
integer iera(0:200) !Positions of additional erasures
|
||||
integer decode_rs_8
|
||||
|
||||
nargs=iargc()
|
||||
if(nargs.ne.3) then
|
||||
print*,'Usage: rst8 npad nera nerr'
|
||||
print*,'Example: rst8 178 0 16'
|
||||
go to 999
|
||||
endif
|
||||
nkv=0
|
||||
call getarg(1,arg)
|
||||
read(arg,*) npad
|
||||
call getarg(2,arg)
|
||||
read(arg,*) nera
|
||||
call getarg(3,arg)
|
||||
read(arg,*) nerr
|
||||
|
||||
! The basic code RS(255,223) is punctured with npad leading zeros.
|
||||
nn=255-npad
|
||||
kk=223-npad
|
||||
write(*,1000) nn,kk
|
||||
1000 format('Basic code is RS(255,223). npad:',i4,' N:',i4,' K:',i4)
|
||||
|
||||
! Generate a message, kk symbols with values 1 to kk.
|
||||
do i=1,kk
|
||||
dat0(i)=i
|
||||
enddo
|
||||
|
||||
write(*,1002)
|
||||
1002 format('Generated message symbols:')
|
||||
write(*,1004) dat0(1:kk)
|
||||
1004 format(20i4)
|
||||
|
||||
call encode_rs_8(dat0,parsym,npad) !Get parity symbols
|
||||
cword0(1:kk)=dat0(1:kk) !Genetated codeword
|
||||
cword0(kk+1:nn)=parsym(1:nn-kk)
|
||||
write(*,1006)
|
||||
1006 format(/'Encoded channel symbols')
|
||||
write(*,1004) cword0(1:nn)
|
||||
|
||||
cword=cword0
|
||||
do i=1,nerr !Introduce errors
|
||||
j=nn+1-i
|
||||
cword(j)=mod(cword(j)+1,256)
|
||||
enddo
|
||||
write(*,1008) nera
|
||||
1008 format(/'Received channel symbols, with',i4,' errors at the end:')
|
||||
write(*,1004) cword(1:nn)
|
||||
|
||||
do i=0,nera-1
|
||||
iera(i)=i
|
||||
enddo
|
||||
|
||||
nfixed=decode_rs_8(cword,iera,nera,npad)
|
||||
ibad=count(cword(1:kk).ne.cword0(1:kk))
|
||||
|
||||
write(*,1010)
|
||||
1010 format(/'Decoded result:')
|
||||
write(*,1004) cword(1:kk)
|
||||
maxfix=(nn-kk)/2 + nera/2
|
||||
write(*,1100) nerr,nera,nfixed,maxfix
|
||||
1100 format(/'nerr:',i3,' nera:',i3,' nfixed:',i3,' maxfix:',i3)
|
||||
|
||||
999 end program rst8
|
||||
@ -1,84 +0,0 @@
|
||||
program rstest
|
||||
|
||||
character arg*8
|
||||
integer dat0(255) !Message symbols
|
||||
integer parsym(255) !Parity symbols
|
||||
integer chansym0(255) !Encoded data, Karn
|
||||
integer chansym(255) !Encoded data with errors
|
||||
integer dat(235) !Decoded data, i*4
|
||||
! integer, target :: parsym(255)
|
||||
integer iera(0:200) !Positions of erasures
|
||||
integer gfpoly
|
||||
! type(c_ptr) :: rs
|
||||
|
||||
data gfpoly/z'11d'/
|
||||
|
||||
nargs=iargc()
|
||||
if(nargs.ne.5) then
|
||||
print*,'Usage: rstest M N K nera nerr'
|
||||
print*,'Examples: rstest 6 63 12 0 25'
|
||||
print*,' rstest 7 127 51 0 38'
|
||||
print*,' rstest 8 255 51 0 102'
|
||||
print*,' rstest 8 255 223 0 16'
|
||||
go to 999
|
||||
endif
|
||||
nkv=0
|
||||
call getarg(1,arg)
|
||||
read(arg,*) mm
|
||||
call getarg(2,arg)
|
||||
read(arg,*) nn
|
||||
call getarg(3,arg)
|
||||
read(arg,*) kk
|
||||
call getarg(4,arg)
|
||||
read(arg,*) nera
|
||||
call getarg(5,arg)
|
||||
read(arg,*) nerr
|
||||
|
||||
! Initialize the Karn codec
|
||||
nq=2**mm
|
||||
nfz=3
|
||||
call rs_init_sf(mm,nq,nn,kk,nfz) !Initialize the Karn RS codec
|
||||
|
||||
! Generate kk message symbols. (Values must be in range 0 to nq-1.)
|
||||
do i=1,kk
|
||||
dat0(i)=i
|
||||
enddo
|
||||
|
||||
write(*,1000) mm,nn,kk,nera,nerr
|
||||
1000 format('M:',i2,' N:',i4,' K:',i4,' nera:',i4,' nerr:',i4/ &
|
||||
'Generated data symbols')
|
||||
write(*,1002) dat0(1:kk)
|
||||
1002 format(20i4)
|
||||
|
||||
call rs_encode_sf(dat0,parsym) !Compute parity symbols
|
||||
chansym0(1:kk)=dat0(1:kk)
|
||||
chansym0(kk+1:nn)=parsym(1:nn-kk)
|
||||
|
||||
write(*,1004)
|
||||
1004 format(/'Encoded channel symbols')
|
||||
write(*,1002) chansym0(1:nn)
|
||||
|
||||
chansym=chansym0
|
||||
do i=1,nerr !Introduce errors
|
||||
chansym(i)=mod(chansym(i)+1,nq)
|
||||
enddo
|
||||
write(*,1006) nera
|
||||
1006 format(/'Recovered channel symbols, with',i4,' errors at the start:')
|
||||
write(*,1002) chansym(1:nn)
|
||||
|
||||
do i=0,nera-1
|
||||
iera(i)=i
|
||||
enddo
|
||||
|
||||
call rs_decode_sf(chansym,iera,nera,nfixed)
|
||||
dat(1:kk)=chansym(1:kk)
|
||||
ibad=count(dat(1:kk).ne.dat0(1:kk))
|
||||
write(*,1008)
|
||||
1008 format(/'Decoded result:')
|
||||
write(*,1002) dat(1:kk)
|
||||
maxfix=(nn-kk)/2 + nera/2
|
||||
write(*,1100) nerr,nera,nfixed,maxfix
|
||||
1100 format(/'nerr:',i3,' nera:',i3,' nfixed:',i3,' maxfix:',i3)
|
||||
|
||||
999 end program rstest
|
||||
|
||||
@ -1,92 +0,0 @@
|
||||
subroutine sfox_assemble(ntype,k,msg,mycall0,mygrid0,line)
|
||||
|
||||
! In subsequent calls, assemble all necessary information for a SuperFox
|
||||
! transmission.
|
||||
|
||||
character*120 line
|
||||
character*26 msg
|
||||
character*26 msg0,msg1,msg2(5),msg3(5)
|
||||
character*4 rpt2(5)
|
||||
character*6 hiscall(10)
|
||||
character*13 mycall0,mycall
|
||||
character*4 mygrid0,mygrid
|
||||
integer ntype !Message type: 0 Free Text
|
||||
! 1 CQ MyCall MyGrid
|
||||
! 2 Call_1 MyCall RR73
|
||||
! 3 Call_1 MyCall rpt
|
||||
integer nmsg(0:3) !Number of messages of type ntype
|
||||
data nmsg/0,0,0,0/,nbits/0/,ntx/0/,nb_mycall/0/
|
||||
save
|
||||
|
||||
if(mycall0(1:1).ne.' ') mycall=mycall0
|
||||
if(mygrid0(1:1).ne.' ') mygrid=mygrid0
|
||||
if(ntype.ge.1) nb_mycall=28 !### Allow for nonstandard MyCall ###
|
||||
if(sum(nmsg).eq.0) then
|
||||
hiscall=' '
|
||||
rpt2=' '
|
||||
endif
|
||||
|
||||
if(k.le.10) then
|
||||
if(ntype.eq.0) then
|
||||
if(nbits+nb_mycall.le.191) then !Enough room for a free text message?
|
||||
nmsg(ntype)=nmsg(ntype)+1
|
||||
nbits=nbits+142
|
||||
msg0=msg
|
||||
endif
|
||||
else if(ntype.eq.1) then
|
||||
if(nbits+nb_mycall.le.318) then !Enough room for a CQ ?
|
||||
nmsg(ntype)=nmsg(ntype)+1
|
||||
nbits=nbits+15
|
||||
msg1=msg
|
||||
endif
|
||||
else if(ntype.eq.2) then
|
||||
if(nbits+nb_mycall.le.305) then !Enough room for a RR73 message?
|
||||
nmsg(ntype)=nmsg(ntype)+1
|
||||
nbits=nbits+28
|
||||
j=nmsg(ntype)
|
||||
msg2(j)=msg
|
||||
i1=index(msg,' ')
|
||||
hiscall(j+5)=msg(1:i1-1)
|
||||
endif
|
||||
else if(ntype.eq.3) then
|
||||
if(nbits+nb_mycall.le.300) then !Enough room for a message with report?
|
||||
nmsg(ntype)=nmsg(ntype)+1
|
||||
nbits=nbits+33
|
||||
j=nmsg(ntype)
|
||||
msg3(j)=msg
|
||||
i1=index(msg,' ')
|
||||
hiscall(j)=msg(1:i1-1)
|
||||
i1=max(index(msg,'-'),index(msg,'+'))
|
||||
rpt2(j)=msg(i1:i1+3)
|
||||
endif
|
||||
endif
|
||||
return
|
||||
endif
|
||||
|
||||
if(k.ge.11) then
|
||||
! All pieces are now available. Put them into a command line for external
|
||||
! program sfox_tx.
|
||||
ntx=ntx+1 !Transmission number
|
||||
nbits=nbits+nb_mycall !Add bits for MyCall
|
||||
|
||||
if(nmsg(1).ge.1) then
|
||||
line=msg1
|
||||
else
|
||||
line=trim(mycall)
|
||||
do i=1,nmsg(3)
|
||||
line=trim(line)//' '//trim(hiscall(i))//' '//rpt2(i)
|
||||
enddo
|
||||
do i=1,nmsg(2)
|
||||
line=trim(line)//' '//trim(hiscall(i+5))
|
||||
enddo
|
||||
endif
|
||||
|
||||
nmsg=0
|
||||
nbits=0
|
||||
nb_mycall=0
|
||||
hiscall=' '
|
||||
rpt2=' '
|
||||
endif
|
||||
|
||||
return
|
||||
end subroutine sfox_assemble
|
||||
@ -1,33 +0,0 @@
|
||||
subroutine sfox_demod(crcvd,f,t,isync,s3)
|
||||
|
||||
use sfox_mod
|
||||
complex crcvd(NMAX) !Signal as received
|
||||
complex c(0:NSPS-1) !Work array, one symbol long
|
||||
real s3(0:NQ-1,0:NN-1) !Synchronized symbol spectra
|
||||
integer isync(44)
|
||||
! integer ipk(1)
|
||||
|
||||
j0=nint(12000.0*(t+0.5))
|
||||
df=12000.0/NSPS
|
||||
i0=nint(f/df)-NQ/2
|
||||
k=-1
|
||||
do n=1,NDS !Loop over all symbols
|
||||
if(any(isync(1:NS).eq.n)) cycle
|
||||
jb=n*NSPS + j0
|
||||
ja=jb-NSPS+1
|
||||
if(ja.lt.1 .or. jb.gt.NMAX) cycle
|
||||
k=k+1
|
||||
c=crcvd(ja:jb)
|
||||
call four2a(c,NSPS,1,-1,1) !Compute symbol spectrum
|
||||
do i=0,NQ-1
|
||||
s3(i,k)=real(c(i0+i))**2 + aimag(c(i0+i))**2
|
||||
enddo
|
||||
! ipk=maxloc(s3(0:NQ-1,k))
|
||||
! if(k.lt.10) print*,'AAA',k,ipk(1)-1
|
||||
enddo
|
||||
|
||||
call pctile(s3,NQ*NN,50,base)
|
||||
s3=s3/base
|
||||
|
||||
return
|
||||
end subroutine sfox_demod
|
||||
@ -1,41 +0,0 @@
|
||||
subroutine sfox_gen(idat,f0,fsample,isync,itone,cdat)
|
||||
|
||||
use sfox_mod
|
||||
complex cdat(NMAX) !Generated complex waveform
|
||||
complex w,wstep
|
||||
integer idat(NN)
|
||||
integer isync(44)
|
||||
integer itone(171)
|
||||
|
||||
twopi=8.0*atan(1.0)
|
||||
|
||||
! Create the itone sequence: data symbols and interspersed sync symbols
|
||||
j=1
|
||||
k=0
|
||||
do i=1,NDS
|
||||
if(j.le.NS .and. i.eq.isync(j)) then
|
||||
if(j.lt.NS) j=j+1 !Index for next sync symbol
|
||||
itone(i)=0 !Insert sync symbol at tone 0
|
||||
else
|
||||
k=k+1
|
||||
itone(i)=idat(k) + 1 !Symbol value 0 is transmitted at tone 1, etc.
|
||||
endif
|
||||
enddo
|
||||
|
||||
df=fsample/NSPS
|
||||
w=1.0
|
||||
j=0
|
||||
i0=NQ/2
|
||||
! Generate the waveform
|
||||
do k=1,NDS !Loop over all symbols
|
||||
dphi=(f0 + (itone(k)-i0)*df)*(twopi/fsample)
|
||||
wstep=cmplx(cos(dphi),sin(dphi))
|
||||
do i=1,NSPS !NSPS samples per symbol
|
||||
j=j+1
|
||||
w=w*wstep
|
||||
cdat(j)=w
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
end subroutine sfox_gen
|
||||
@ -1,73 +0,0 @@
|
||||
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
|
||||
@ -1,153 +0,0 @@
|
||||
subroutine sfox_sync(iwave,fsample,isync,f,t,fwidth)
|
||||
|
||||
use sfox_mod
|
||||
parameter (NSTEP=8)
|
||||
integer*2 iwave(0:NMAX-1)
|
||||
integer isync(44)
|
||||
integer ipeak(2)
|
||||
integer ipeak2(1)
|
||||
complex, allocatable :: c(:) !Work array
|
||||
real, allocatable :: s(:,:) !Symbol spectra, stepped by NSTEP
|
||||
real, allocatable :: savg(:) !Average spectrum
|
||||
real, allocatable :: ccf(:,:)
|
||||
real, allocatable :: s2(:) !Fine spectrum of sync tone
|
||||
|
||||
nfft=nsps
|
||||
nh=nfft/2
|
||||
istep=NSPS/NSTEP
|
||||
jz=(13.5*fsample)/istep
|
||||
df=fsample/nfft
|
||||
dtstep=istep/fsample
|
||||
fsync=1500.0-bw/2
|
||||
ftol=50.0
|
||||
ia=nint((fsync-ftol)/df)
|
||||
ib=nint((fsync+ftol)/df)
|
||||
lagmax=1.5/dtstep
|
||||
lag1=-lagmax
|
||||
lag2=lagmax
|
||||
|
||||
allocate(s(0:nh/2,jz))
|
||||
allocate(savg(0:nh/2))
|
||||
allocate(c(0:nfft-1))
|
||||
allocate(ccf(ia:ib,lag1:lag2))
|
||||
|
||||
s=0.
|
||||
savg=0.
|
||||
fac=1.0/nfft
|
||||
|
||||
! Compute symbol spectra with df=baud/2 and NSTEP steps per symbol.
|
||||
do j=1,jz
|
||||
i1=(j-1)*istep
|
||||
i2=i1+nsps-1
|
||||
k=-1
|
||||
do i=i1,i2,2 !Load iwave data into complex array c0, for r2c FFT
|
||||
xx=iwave(i)
|
||||
yy=iwave(i+1)
|
||||
k=k+1
|
||||
c(k)=fac*cmplx(xx,yy)
|
||||
enddo
|
||||
c(k+1:)=0.
|
||||
call four2a(c,nfft,1,-1,0) !r2c FFT
|
||||
do i=1,nh/2
|
||||
s(i,j)=real(c(i))**2 + aimag(c(i))**2
|
||||
savg(i)=savg(i) + s(i,j)
|
||||
enddo
|
||||
enddo
|
||||
savg=savg/jz
|
||||
|
||||
ccfbest=0.
|
||||
ibest=0
|
||||
lagpk=0
|
||||
lagbest=0
|
||||
j0=0.5/dtstep !Nominal start-signal index
|
||||
|
||||
do i=ia,ib
|
||||
ccfmax=0.
|
||||
do lag=lag1,lag2
|
||||
ccft=0.
|
||||
do m=1,NS
|
||||
k=isync(m)
|
||||
n=NSTEP*(k-1) + 1
|
||||
j=n+lag+j0
|
||||
if(j.ge.1 .and. j.le.jz) ccft=ccft + s(i,j)
|
||||
enddo ! m
|
||||
ccft=ccft - NS*savg(i)
|
||||
ccf(i,lag)=ccft
|
||||
if(ccft.gt.ccfmax) then
|
||||
ccfmax=ccft
|
||||
lagpk=lag
|
||||
endif
|
||||
enddo ! lag
|
||||
|
||||
if(ccfmax.gt.ccfbest) then
|
||||
ccfbest=ccfmax
|
||||
ibest=i
|
||||
lagbest=lagpk
|
||||
endif
|
||||
enddo ! i
|
||||
|
||||
ipeak=maxloc(ccf)
|
||||
ipk=ipeak(1)-1+ia
|
||||
jpk=ipeak(2)-1+lag1
|
||||
|
||||
dxj=0.
|
||||
if(jpk.gt.lag1 .and. jpk.lt.lag2) then
|
||||
call peakup(ccf(ipk,jpk-1),ccf(ipk,jpk),ccf(ipk,jpk+1),dxj)
|
||||
endif
|
||||
|
||||
f=ibest*df + bw/2 + dxi*df
|
||||
t=(lagbest+dxj)*dtstep
|
||||
t=t-0.01 !### Why is this needed? ###
|
||||
|
||||
nfft2=4*NSPS
|
||||
deallocate(c)
|
||||
allocate(c(0:nfft2-1))
|
||||
allocate(s2(0:nfft2-1))
|
||||
|
||||
i0=(t+0.5)*fsample
|
||||
s2=0.
|
||||
df2=fsample/nfft2
|
||||
do m=1,NS
|
||||
i1=i0+(isync(m)-1)*NSPS
|
||||
i2=i1+NSPS-1
|
||||
k=-1
|
||||
do i=i1,i2,2 !Load iwave data into complex array c0, for r2c FFT
|
||||
if(i.gt.0) then
|
||||
xx=iwave(i)
|
||||
yy=iwave(i+1)
|
||||
else
|
||||
xx=0.
|
||||
yy=0.
|
||||
endif
|
||||
k=k+1
|
||||
c(k)=fac*cmplx(xx,yy)
|
||||
enddo
|
||||
c(k+1:)=0.
|
||||
call four2a(c,nfft2,1,-1,0) !r2c FFT
|
||||
do i=1,nfft2/4
|
||||
s2(i)=s2(i) + real(c(i))**2 + aimag(c(i))**2
|
||||
enddo
|
||||
enddo
|
||||
|
||||
ia=nint((fsync-ftol)/df2)
|
||||
ib=nint((fsync+ftol)/df2)
|
||||
ipeak2=maxloc(s2(ia:ib))
|
||||
ipk=ipeak2(1)-1+ia
|
||||
|
||||
dxi=0.
|
||||
if(ipk.gt.1 .and. ipk.lt.nfft/4) then
|
||||
call peakup(s2(ipk-1),s2(ipk),s2(ipk+1),dxi)
|
||||
endif
|
||||
f=(ipk+dxi)*df2 + bw/2.0
|
||||
fwidth=0.
|
||||
|
||||
if(ipk.gt.100 .and. ipk.lt.nfft2/4-100) then
|
||||
call pctile(s2(ipk-100:ipk+100),201,48,base)
|
||||
s2=s2-base
|
||||
smax=maxval(s2(ipk-10:ipk+10))
|
||||
w=count(s2(ipk-10:ipk+10).gt.0.5*smax)
|
||||
if(w.gt.4.0) fwidth=sqrt(w*w - 4*4)*df2
|
||||
endif
|
||||
|
||||
return
|
||||
end subroutine sfox_sync
|
||||
@ -1,46 +0,0 @@
|
||||
program sfox_tx
|
||||
|
||||
! This program is required in order to create a SuperFox transmission.
|
||||
|
||||
! The present version goes through the following steps:
|
||||
! 1. Read old-style Fox messages from file 'sfox_1.dat' in the WSJT-X
|
||||
! writable data directory.
|
||||
! 2. Parse up to NSlots=5 messages to extract MyCall, up to 10 Hound
|
||||
! calls, and the report or RR73 to be sent to each Hound.
|
||||
! 3. Assemble and encode a single SuperFox message to produce itone(1:151),
|
||||
! the array of channel symbol values.
|
||||
! 4. Write the contents of array itone to file 'sfox_2.dat'.
|
||||
|
||||
character*120 fname !Full path for sfox.dat
|
||||
character*120 line !List of SuperFox message pieces
|
||||
character*40 cmsg(5) !Old-style Fox messages
|
||||
integer itone(151) !SuperFox channel-symbol values
|
||||
|
||||
open(70,file='fort.70',status='unknown',position='append')
|
||||
call getarg(1,fname)
|
||||
open(25,file=trim(fname),status='unknown')
|
||||
do i=1,5
|
||||
read(25,1000,end=10) cmsg(i)
|
||||
1000 format(a40)
|
||||
write(70,*) 'AAA',i,cmsg(i)
|
||||
enddo
|
||||
i=6
|
||||
|
||||
10 close(25)
|
||||
nslots=i-1
|
||||
write(70,*) 'BBB',nslots
|
||||
call foxgen2(nslots,cmsg,line)
|
||||
write(70,*) 'CCC ',trim(line)
|
||||
|
||||
do i=1,151 !Dummy loop to populate itone during tests
|
||||
itone(i)=i-1
|
||||
enddo
|
||||
|
||||
i1=index(fname,'sfox_1.dat')
|
||||
fname(i1:i1+9)='sfox_2.dat'
|
||||
open(25,file=trim(fname),status='unknown')
|
||||
write(25,1100) itone
|
||||
1100 format(20i4)
|
||||
close(25)
|
||||
|
||||
end program sfox_tx
|
||||
@ -1,116 +0,0 @@
|
||||
program sfoxsim
|
||||
|
||||
! Generate a SuperFox waveform with specified SNR and channel parameters.
|
||||
! Output is saved to a *.wav file.
|
||||
! SuperFox uses a (127,51) code with 7-bit symbols, punctured to (125,49).
|
||||
! The puncured symbols contain a 14-bit CRC.
|
||||
! First tests use RS(127,51) code and Berlekamp-Massey decoder.
|
||||
|
||||
use wavhdr
|
||||
! use packjt77
|
||||
parameter (NMAX=15*12000)
|
||||
parameter (NSPS=1024,NSYNC=2*12000)
|
||||
parameter (NWAVE=125*NSPS+NSYNC)
|
||||
type(hdr) h !Header for .wav file
|
||||
character arg*12,fname*17
|
||||
! character msg37*37,msgsent37*37
|
||||
complex c0(0:NMAX-1)
|
||||
complex c(0:NMAX-1)
|
||||
complex cwave(0:NWAVE-1)
|
||||
real wave(NMAX)
|
||||
real xjunk(NWAVE)
|
||||
real xdat(51)
|
||||
integer*1 idat(51)
|
||||
integer itone(125)
|
||||
integer*1 msgbits(77)
|
||||
integer*2 iwave(NMAX) !Generated full-length waveform
|
||||
|
||||
! Get command-line argument(s)
|
||||
nargs=iargc()
|
||||
if(nargs.ne.6) then
|
||||
print*,'Usage: sfoxsim f0 DT fSpread del nfiles snr'
|
||||
print*,'Example: sfoxsim 1500.0 0.0 0.1 1.0 10 -15'
|
||||
go to 999
|
||||
endif
|
||||
call getarg(1,arg)
|
||||
read(arg,*) f0 !Frequency (only used for single-signal
|
||||
call getarg(2,arg)
|
||||
read(arg,*) xdt !Time offset from nominal (s)
|
||||
call getarg(3,arg)
|
||||
read(arg,*) fspread !Watterson frequency spread (Hz)
|
||||
call getarg(4,arg)
|
||||
read(arg,*) delay !Watterson delay (ms)
|
||||
call getarg(5,arg)
|
||||
read(arg,*) nfiles !Number of files
|
||||
call getarg(6,arg)
|
||||
read(arg,*) snrdb !SNR_2500
|
||||
|
||||
twopi=8.0*atan(1.0)
|
||||
fs=12000.0 !Sample rate (Hz)
|
||||
dt=1.0/fs !Sample interval (s)
|
||||
tt=NSPS*dt !Duration of symbols (s)
|
||||
baud=1.0/tt !Keying rate (baud)
|
||||
bw=128*baud !Occupied bandwidth (Hz)
|
||||
tsync=NSYNC*dt !Duration of analog sync function
|
||||
txt=tsync + 125*NSPS*dt !Overall transmission length (s)
|
||||
bandwidth_ratio=2500.0/(fs/2.0)
|
||||
sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
|
||||
if(snrdb.gt.90.0) sig=1.0
|
||||
|
||||
write(*,1000) f0,xdt,fspread,delay,tsync,txt,bw,snrdb
|
||||
1000 format('f0:',f7.1,' DT:',f6.2,' fSpread:',f5.1,' delay:',f4.1/ &
|
||||
'Tsync:',f4.1,' TxT:',f6.1,' BW:',f7.1,' SNR:',f6.1)
|
||||
write(*,*)
|
||||
|
||||
! Source-encode, then get itone()
|
||||
|
||||
call random_number(xdat)
|
||||
idat=int(128*xdat)
|
||||
itone=0
|
||||
itone(1:49)=idat(1:49)
|
||||
|
||||
write(*,'(20i4)') idat
|
||||
write(*,*)
|
||||
write(*,'(a17)') 'Channel symbols: '
|
||||
write(*,'(20i4)') itone
|
||||
write(*,*)
|
||||
if(nsps.ne.-99) go to 999
|
||||
|
||||
do ifile=1,nfiles
|
||||
c0=0.
|
||||
c0(0:NWAVE-1)=cwave
|
||||
c0=cshift(c0,-nint((xdt+0.5)/dt))
|
||||
! if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c0,NMAX,NWAVE,fs,delay,fspread)
|
||||
c=sig*c0
|
||||
|
||||
wave=imag(c)
|
||||
peak=maxval(abs(wave))
|
||||
nslots=1
|
||||
|
||||
if(snrdb.lt.90) then
|
||||
do i=1,NMAX !Add gaussian noise at specified SNR
|
||||
xnoise=gran()
|
||||
wave(i)=wave(i) + xnoise
|
||||
enddo
|
||||
endif
|
||||
|
||||
gain=100.0
|
||||
if(snrdb.lt.90.0) then
|
||||
wave=gain*wave
|
||||
else
|
||||
datpk=maxval(abs(wave))
|
||||
fac=32766.9/datpk
|
||||
wave=fac*wave
|
||||
endif
|
||||
if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped."
|
||||
iwave=nint(wave)
|
||||
h=default_header(12000,NMAX)
|
||||
write(fname,1102) ifile
|
||||
1102 format('000000_',i6.6,'.wav')
|
||||
open(10,file=fname,status='unknown',access='stream')
|
||||
write(10) h,iwave !Save to *.wav file
|
||||
close(10)
|
||||
write(*,1110) ifile,xdt,f0,snrdb,fname
|
||||
1110 format(i4,f7.2,f8.2,f7.1,2x,a17)
|
||||
enddo
|
||||
999 end program sfoxsim
|
||||
@ -1,285 +0,0 @@
|
||||
program sfoxtest
|
||||
|
||||
! Generate and test possible waveforms for SuperFox signal.
|
||||
|
||||
use wavhdr
|
||||
use sfox_mod
|
||||
use timer_module, only: timer
|
||||
use timer_impl, only: init_timer !, limtrace
|
||||
|
||||
type(hdr) h !Header for .wav file
|
||||
integer*2 iwave(NMAX) !Generated i*2 waveform
|
||||
integer param(0:8)
|
||||
integer isync(44)
|
||||
integer jsync(171)
|
||||
integer itone(171)
|
||||
integer nsb(10)
|
||||
real*4 xnoise(NMAX) !Random noise
|
||||
real*4 dat(NMAX) !Generated real data
|
||||
complex cdat(NMAX) !Generated complex waveform
|
||||
complex cnoise(NMAX) !Complex noise
|
||||
complex crcvd(NMAX) !Signal as received
|
||||
real a(3)
|
||||
real, allocatable :: s3(:,:) !Symbol spectra: will be s3(NQ,NN)
|
||||
integer, allocatable :: msg0(:) !Information symbols
|
||||
integer, allocatable :: parsym(:) !Parity symbols
|
||||
integer, allocatable :: chansym0(:) !Encoded data, 7-bit integers
|
||||
integer, allocatable :: chansym(:) !Recovered hard-decision symbols
|
||||
integer, allocatable :: iera(:) !Positions of erasures
|
||||
integer, allocatable :: rxdat(:)
|
||||
integer, allocatable :: rxprob(:)
|
||||
integer, allocatable :: rxdat2(:)
|
||||
integer, allocatable :: rxprob2(:)
|
||||
integer, allocatable :: correct(:)
|
||||
character fname*17,arg*12,itu*2
|
||||
|
||||
data isync/ 1, 2, 5, 11, 19, 24, 26, 28, 29, 35, &
|
||||
39, 48, 51, 53, 55, 56, 66, 71, 74, 78, &
|
||||
80, 82, 84, 85, 92, 98, 103, 107, 109, 111, &
|
||||
116, 122, 130, 131, 134, 136, 137, 140, 146, 154, &
|
||||
159, 161, 163, 165/
|
||||
data nsb/1,2,4,7,11,16,22,29,37,39/
|
||||
|
||||
nargs=iargc()
|
||||
if(nargs.ne.11) then
|
||||
print*,'Usage: sfoxtest f0 DT ITU M N K NS v st nfiles snr'
|
||||
print*,'Example: sfoxtest 1500 0.15 MM 7 127 48 33 0 3 10 -10'
|
||||
print*,' f0=0 means f0, DT will assume suitable random values'
|
||||
print*,' LQ: Low Latitude Quiet'
|
||||
print*,' MM: Mid Latitude Moderate'
|
||||
print*,' HD: High Latitude Disturbed'
|
||||
print*,' ... and similarly for LM LD MQ MD HQ HM'
|
||||
print*,' NS: number of sync symbols'
|
||||
print*,' v=1 for .wav files, 2 for verbose output, 3 for both'
|
||||
print*,' st: Sync type, 0 for hard-wired, otherwise 1-3'
|
||||
print*,' snr=0 means loop over SNRs'
|
||||
go to 999
|
||||
endif
|
||||
call getarg(1,arg)
|
||||
read(arg,*) f0
|
||||
call getarg(2,arg)
|
||||
read(arg,*) xdt
|
||||
call getarg(3,itu)
|
||||
call getarg(4,arg)
|
||||
read(arg,*) mm0
|
||||
call getarg(5,arg)
|
||||
read(arg,*) nn0
|
||||
call getarg(6,arg)
|
||||
read(arg,*) kk0
|
||||
call getarg(7,arg)
|
||||
read(arg,*) ns0
|
||||
call getarg(8,arg)
|
||||
read(arg,*) nv
|
||||
call getarg(9,arg)
|
||||
read(arg,*) nstype
|
||||
call getarg(10,arg)
|
||||
read(arg,*) nfiles
|
||||
call getarg(11,arg)
|
||||
read(arg,*) snrdb
|
||||
|
||||
call init_timer ('timer.out')
|
||||
call timer('sfoxtest',0)
|
||||
|
||||
fsample=12000.0 !Sample rate (Hz)
|
||||
call sfox_init(mm0,nn0,kk0,itu,fspread,delay,fsample,ns0)
|
||||
tsync=NSYNC/fsample
|
||||
txt=(NN+NS)*NSPS/fsample
|
||||
|
||||
write(*,1000) MM,NN,KK,NSPS,baud,bw,itu,NS,txt,nstype
|
||||
1000 format('M:',i2,' Base code: (',i3,',',i3,') NSPS:',i5, &
|
||||
' Baud:',f7.3,' BW:',f9.3/ &
|
||||
'Channel: ',a2,' NS:',i3,' TxT:',f5.1,' SyncType:',i2/)
|
||||
|
||||
! Allocate storage for arrays that depend on code parameters.
|
||||
allocate(s3(0:NQ-1,0:NN-1))
|
||||
allocate(msg0(1:KK))
|
||||
allocate(parsym(1:NN-KK))
|
||||
allocate(chansym0(0:NN-1))
|
||||
allocate(chansym(0:NN-1))
|
||||
allocate(iera(0:NN-1))
|
||||
allocate(rxdat(0:NN-1))
|
||||
allocate(rxprob(0:NN-1))
|
||||
allocate(rxdat2(0:NN-1))
|
||||
allocate(rxprob2(0:NN-1))
|
||||
allocate(correct(0:NN-1))
|
||||
|
||||
idum=-1
|
||||
if(nstype.eq.2) then
|
||||
jsync=0
|
||||
jsync(1)=1
|
||||
jsync(NDS)=1
|
||||
ms=2
|
||||
do i=1,100000
|
||||
j=1 + (NDS-1)*ran1(idum)
|
||||
if(jsync(j).eq.0) then
|
||||
jsync(j)=1
|
||||
ms=ms+1
|
||||
if(ms.eq.NS) exit
|
||||
endif
|
||||
enddo
|
||||
j=0
|
||||
do i=1,NDS
|
||||
if(jsync(i).eq.1) then
|
||||
j=j+1
|
||||
isync(j)=i
|
||||
endif
|
||||
enddo
|
||||
else if(nstype.eq.3) then
|
||||
isync(1:10)=nsb
|
||||
isync(11:20)=nsb + isync(10) + 2
|
||||
isync(21:30)=nsb + isync(20) + 2
|
||||
isync(31:40)=nsb + isync(30) + 2
|
||||
isync(41:44)=nsb(1:4) + isync(40) + 2
|
||||
endif
|
||||
|
||||
rms=100.
|
||||
baud=fsample/nsps !Keying rate, 11.719 baud for nsps=1024
|
||||
bandwidth_ratio=2500.0/fsample
|
||||
fgood0=1.0
|
||||
|
||||
! Generate a message
|
||||
msg0=0
|
||||
do i=1,KK
|
||||
msg0(i)=int(NQ*ran1(idum))
|
||||
enddo
|
||||
|
||||
call rs_init_sf(MM,NQ,NN,KK,NFZ) !Initialize the Karn codec
|
||||
call rs_encode_sf(msg0,parsym) !Compute parity symbols
|
||||
chansym0(0:kk-1)=msg0(1:kk)
|
||||
chansym0(kk:nn-1)=parsym(1:nn-kk)
|
||||
|
||||
isnr0=-8
|
||||
|
||||
do isnr=isnr0,-20,-1
|
||||
snr=isnr
|
||||
if(snrdb.ne.0.0) snr=snrdb
|
||||
sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snr)
|
||||
sigr=sqrt(2.)*sig
|
||||
if(snr.gt.90.0) sig=1.0
|
||||
ngoodsync=0
|
||||
ngood=0
|
||||
ntot=0
|
||||
nworst=0
|
||||
sqt=0.
|
||||
sqf=0.
|
||||
sumw=0.
|
||||
|
||||
do ifile=1,nfiles
|
||||
xnoise=0.
|
||||
cnoise=0.
|
||||
if(snr.lt.90) then
|
||||
do i=1,NMAX !Generate Gaussian noise
|
||||
x=gran()
|
||||
y=gran()
|
||||
xnoise(i)=x
|
||||
cnoise(i)=cmplx(x,y)
|
||||
enddo
|
||||
endif
|
||||
|
||||
f1=f0
|
||||
if(f0.eq.0.0) then
|
||||
f1=1500.0 + 20.0*(ran1(idum)-0.5)
|
||||
xdt=ran1(idum)-0.5
|
||||
endif
|
||||
call timer('gen ',0)
|
||||
! Generate cdat, the SuperFox waveform
|
||||
call sfox_gen(chansym0,f1,fsample,isync,itone,cdat)
|
||||
call timer('gen ',1)
|
||||
|
||||
crcvd=0.
|
||||
crcvd(1:NMAX)=cshift(cdat(1:NMAX),-nint((0.5+xdt)*fsample))
|
||||
call timer('watterso',0)
|
||||
if(fspread.ne.0 .or. delay.ne.0) call watterson(crcvd,NMAX,NZ,fsample,&
|
||||
delay,fspread)
|
||||
call timer('watterso',1)
|
||||
|
||||
dat=aimag(sigr*crcvd(1:NMAX)) + xnoise !Add generated AWGN noise
|
||||
fac=32767.0
|
||||
if(snr.ge.90.0) iwave(1:NMAX)=nint(fac*dat(1:NMAX))
|
||||
if(snr.lt.90.0) iwave(1:NMAX)=nint(rms*dat(1:NMAX))
|
||||
|
||||
crcvd=sig*crcvd+cnoise
|
||||
if(nstype.eq.0) then
|
||||
f=f1 !Hard-wired sync
|
||||
t=xdt
|
||||
else
|
||||
call timer('sync ',0)
|
||||
call sfox_sync(iwave,fsample,isync,f,t,fwidth) !Find freq, DT, width
|
||||
call timer('sync ',1)
|
||||
endif
|
||||
ferr=f-f1
|
||||
terr=t-xdt
|
||||
|
||||
igoodsync=0
|
||||
if(abs(ferr).lt.baud/2.0 .and. abs(terr).lt.tsym/4.0) then
|
||||
igoodsync=1
|
||||
ngoodsync=ngoodsync+1
|
||||
sqt=sqt + terr*terr
|
||||
sqf=sqf + ferr*ferr
|
||||
sumw=sumw+fwidth
|
||||
endif
|
||||
|
||||
a=0.
|
||||
a(1)=1500.0-f - baud !Shift frequencies down by one bin
|
||||
call timer('twkfreq ',0)
|
||||
call twkfreq(crcvd,crcvd,NMAX,fsample,a)
|
||||
call timer('twkfreq ',1)
|
||||
f=1500.0
|
||||
call timer('demod ',0)
|
||||
call sfox_demod(crcvd,f,t,isync,s3) !Get s3(0:NQ-1,0:127)
|
||||
call timer('demod ',1)
|
||||
|
||||
call timer('prob ',0)
|
||||
call sym_prob(s3,rxdat,rxprob,rxdat2,rxprob2)
|
||||
call timer('prob ',1)
|
||||
|
||||
nera=0
|
||||
chansym=mod(chansym,nq) !Enforce 0 to nq-1
|
||||
nharderr=count(rxdat.ne.chansym0) !Count hard errors
|
||||
ntot=ntot+nharderr
|
||||
nworst=max(nworst,nharderr)
|
||||
|
||||
ntrials=1000
|
||||
call timer('ftrsd3 ',0)
|
||||
call ftrsd3(s3,chansym0,rxdat,rxprob,rxdat2,rxprob2,ntrials, &
|
||||
correct,param,ntry)
|
||||
call timer('ftrsd3 ',1)
|
||||
|
||||
if(iand(nv,1).ne.0) then
|
||||
h=default_header(12000,NMAX)
|
||||
fname='000000_000001.wav'
|
||||
write(fname(8:13),'(i6.6)') ifile
|
||||
open(10,file=trim(fname),access='stream',status='unknown')
|
||||
write(10) h,iwave(1:NMAX) !Save the .wav file
|
||||
close(10)
|
||||
endif
|
||||
|
||||
if(count(correct.ne.chansym0).eq.0) ngood=ngood+1
|
||||
enddo ! ifile
|
||||
fgoodsync=float(ngoodsync)/nfiles
|
||||
fgood=float(ngood)/nfiles
|
||||
if(isnr.eq.isnr0) write(*,1300)
|
||||
1300 format(' SNR Eb/No iters fsync fgood averr rmsf rmst avew'/ &
|
||||
'-------------------------------------------------------------')
|
||||
ave_harderr=float(ntot)/nfiles
|
||||
rmst=sqrt(sqt/ngoodsync)
|
||||
rmsf=sqrt(sqf/ngoodsync)
|
||||
ebno=snr-10*log10(baud/2500*mm0*KK/NN)
|
||||
avew=sumw/ngoodsync
|
||||
write(*,1310) snr,ebno,nfiles,fgoodsync,fgood,ave_harderr,rmsf,rmst,avew
|
||||
1310 format(f7.2,f7.2 i6,2f7.4,f7.1,f7.2,f7.4,f6.1)
|
||||
if(fgood.le.0.5 .and. fgood0.gt.0.5) then
|
||||
threshold=isnr + 1 - (fgood0-0.50)/(fgood0-fgood+0.000001)
|
||||
endif
|
||||
fgood0=fgood
|
||||
if(snrdb.ne.0.0) exit
|
||||
! if(fgood.eq.0.0) exit
|
||||
if(fgoodsync.lt.0.5) exit
|
||||
enddo ! isnr
|
||||
if(snrdb.eq.0.0) write(*,1320) threshold
|
||||
1320 format(/'Threshold sensitivity (50% decoding):',f6.1,' dB')
|
||||
call timer('sfoxtest',1)
|
||||
|
||||
999 call timer('sfoxtest',101)
|
||||
end program sfoxtest
|
||||
|
||||
@ -1,55 +0,0 @@
|
||||
subroutine sym_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 sym_prob
|
||||
@ -1,18 +0,0 @@
|
||||
#include <stdio.h>
|
||||
|
||||
void encode_rs_8(unsigned char *data, unsigned char *parity, int pad);
|
||||
|
||||
int decode_rs_8(unsigned char *data, int *eras_pos, int no_eras, int pad);
|
||||
|
||||
void encode_rs_8_(unsigned char data[], unsigned char parity[], int *npad)
|
||||
{
|
||||
encode_rs_8(data,parity,*npad); //Compute the parity symbols
|
||||
}
|
||||
|
||||
|
||||
int decode_rs_8_(unsigned char *data, int *era_pos, int *numera, int *npad)
|
||||
{
|
||||
int nerr;
|
||||
nerr=decode_rs_8(data,era_pos,*numera,*npad);
|
||||
return nerr;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user