Add new file ftrsd3.f90.

CMakeLists.txt

@@ -592,6 +592,7 @@ set (wsjt_FSRCS
   lib/superfox/sfox_demod.f90
   lib/superfox/sfox_clo.f90
   lib/superfox/sym_prob.f90
+  lib/superfox/getpp3.f90
   lib/superfox/ran1.f90
   )
 
@@ -637,6 +638,7 @@ set (wsjt_CSRCS
   lib/superfox/encode_rs.c
   lib/superfox/decode_rs.c
   lib/superfox/rs_sf.c
+  lib/superfox/ftrsd3.c
   ${ldpc_CSRCS}
   ${qra_CSRCS}
   )

lib/superfox/ftrsd3.f90

@@ -0,0 +1,189 @@
+subroutine ftrsd3(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
+
+  integer, dimension(0:NN-1) :: rxdat,rxprob,rxdat2,rxprob2,workdat,   &
+       correct,indexes
+  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*8 nseed,ir              !No unsigned int in Fortran
+  integer pass,tmp
+  
+  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
+           probs(k)=probs(k+1)
+           probs(k+1)=tmp
+           tmp=indexes(k)
+           indexes(k)=indexes(k+1)
+           indexes(k+1)=tmp
+        enddo
+     enddo
+  enddo
+
+  era_pos=0
+  numera=0
+  workdat=rxdat
+  call rs_decode_sf(workdat,era_pos,numera,nerr)    !Call the decoder
+
+  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
+     return
+  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=(NN-1-i)/8
+     thresh0(i)=1.3*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.
+
+     numera=0
+     do i=0,NN-1
+        j=indexes(126-i)
+        thresh=thresh0(i)
+
+! Generate a random number ir, 0 <= ir <= 100 (see POSIX.1-2001 example).
+        nseed=nseed*1103515245 + 12345
+        ir=mod(nseed/65536),32768)
+        ir=(100*ir)/32768
+        nseed=iand(ir,4294967295)
+
+        if((ir.lt.thresh ) .and. numera.lt.(NN-KK)) then
+           era_pos(numera)=j
+           numera=numera+1
+        endif
+     enddo
+
+!     nerr=decode_rs_int(rs,workdat,era_pos,numera,0);
+     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(64,127) 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-1)*nsoft/nsum
+        ntotal=nsoft+nhard
+
+        pp=0.
+        call getpp3(s3,workdat,pp)
+        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 <= 41 && ntotal_min <= 71) exit   !### New values ###
+     enddo
+     if(k.eq.ntrials) ntry=k
+  enddo
+
+  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
+  param(7)=1000.0*pp2
+  param(8)=1000.0*pp1
+  if(param(0).eq.0) param(2)=-1
+
+  return
+end subroutine ftrsd3
+

lib/superfox/sfox_demod.f90

@@ -4,15 +4,14 @@ subroutine sfox_demod(crcvd,f,t,s3,chansym)
   complex crcvd(NMAX)                    !Signal as received
   complex c(0:NSPS-1)                    !Work array, one symbol long
   real s(0:NQ-1)                         !Power spectrum
-  real s3(0:NQ-1,0:ND-1)                 !Symbol spectra
+  real s3(0:NQ-1,0:NN-1)                 !Symbol spectra
   integer chansym(NN)                    !Hard-decision symbol values
   integer ipk(1)
 
   i0=nint(12000.0*t)
   df=12000.0/NSPS
   j0=nint(f/df)-NQ/2
-  chansym=0
+  do n=1,NN                             !Loop over all symbols
-  do n=1,ND                             !Loop over all symbols
      ib=n*NSPS + i0
      if(n.gt.ND1) ib=(NS+n)*NSPS + i0
      ia=ib-NSPS+1

lib/superfox/sfox_gen.f90

@@ -4,7 +4,7 @@ subroutine sfox_gen(idat,f0,fsample,syncwidth,cdat)
 !  include "sfox_params.f90"
   complex cdat(NMAX)                     !Generated complex waveform
   complex w,wstep
-  integer idat(ND)
+  integer idat(NN)
 
   twopi=8.0*atan(1.0)
   tsync=NS*NSPS/fsample

lib/superfox/sfox_mod.f90

@@ -1,7 +1,7 @@
 module sfox_mod
   
   parameter (NMAX=15*12000)       !Samples in iwave (180,000)
-  integer MM,NQ,NN,KK,ND1,ND2,ND,NFZ,NSPS,NS,NSYNC,NZ,NFFT,NFFT1
+  integer MM,NQ,NN,KK,ND1,ND2,NFZ,NSPS,NS,NSYNC,NZ,NFFT,NFFT1
 
 contains
   subroutine sfox_init(mm0,nn0,kk0,itu,fspread,delay)
@@ -17,21 +17,20 @@ contains
 
     MM=mm0              !Bits per symbol
     NQ=2**MM            !Q, number of MFSK tones
-    NN=nn0              !Channel symbols, before puncture
+    NN=nn0              !Number of channel symbols
-    KK=kk0              !Information symbols, before puncture
+    KK=kk0              !Information symbols
     ND1=25              !Data symbols before sync 
     ND2=NN-ND1          !Data symbols after sync 
-    ND=NN               !Total data symbols
     NFZ=3               !First zero
 
     tsync=2.0
-    jsps=nint((12.8-tsync)*12000.0/ND)
+    jsps=nint((12.8-tsync)*12000.0/NN)
     iloc=minloc(abs(isps-jsps))
     NSPS=isps(iloc(1))  !Samples per symbol
     NS=nint(tsync*12000.0/NSPS)
     if(mod(NS,2).eq.1) NS=NS+1
     NSYNC=NS*NSPS       !Samples in sync waveform
-    NZ=NSPS*(ND+NS)     !Samples in full Tx waveform
+    NZ=NSPS*(NN+NS)     !Samples in full Tx waveform
     NFFT=32768          !Length of FFT for sync waveform
     NFFT1=2*NSPS        !Length of FFTs for symbol spectra
 

lib/superfox/sfoxtest.f90

@@ -6,6 +6,7 @@ program sfoxtest
   use sfox_mod
   type(hdr) h                            !Header for .wav file
   integer*2 iwave(NMAX)                  !Generated i*2 waveform
+  integer nparam(0:7)
   real*4 xnoise(NMAX)                    !Random noise
   real*4 dat(NMAX)                       !Generated real data
   complex cdat(NMAX)                     !Generated complex waveform
@@ -13,7 +14,7 @@ program sfoxtest
   complex cnoise(NMAX)                   !Complex noise
   complex crcvd(NMAX)                    !Signal as received
   real a(3)
-  real, allocatable :: s3(:,:)           !Symbol spectra: will be s3(NQ,ND)
+  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
@@ -23,6 +24,8 @@ program sfoxtest
   integer, allocatable :: rxprob(:)
   integer, allocatable :: rxdat2(:)
   integer, allocatable :: rxprob2(:)
+  integer, allocatable :: correct(:)
+  
   character fname*17,arg*12,itu*2
 
   nargs=iargc()
@@ -75,16 +78,17 @@ program sfoxtest
           '   MaxErr:',i3,'  tsync:',f4.1,'   TxT:',f5.1/)
 
 ! Allocate storage for arrays that depend on code parameters.
-  allocate(s3(0:NQ-1,0:ND-1))
+  allocate(s3(0:NQ-1,0:NN-1))
   allocate(msg0(1:KK))
   allocate(parsym(1:NN-KK))
   allocate(chansym0(1:NN))
   allocate(chansym(1:NN))
   allocate(iera(1:NN))
-  allocate(rxdat(0:ND-1))
+  allocate(rxdat(0:NN-1))
-  allocate(rxprob(0:ND-1))
+  allocate(rxprob(0:NN-1))
-  allocate(rxdat2(0:ND-1))
+  allocate(rxdat2(0:NN-1))
-  allocate(rxprob2(0:ND-1))
+  allocate(rxprob2(0:NN-1))
+  allocate(correct(0:NN-1))
 
   rms=100.
   fsample=12000.0                   !Sample rate (Hz)
@@ -138,7 +142,8 @@ program sfoxtest
         f1=f0
         if(f0.eq.0.0) then
            f1=1500.0 + 200.0*(ran1(idummy)-0.5)
-           xdt=0.6*(ran1(idummy)-0.5)
+!           xdt=0.6*(ran1(idummy)-0.5)
+           xdt=0.3*ran1(idummy)
            call sfox_gen(chansym0,f1,fsample,syncwidth,cdat,clo)
         endif
         
@@ -157,10 +162,15 @@ program sfoxtest
         call sfox_sync(crcvd,clo,nv,f,t)
         ferr=f-f1
         terr=t-xdt
+        igoodsync=0
         if(abs(ferr).lt.baud/2.0 .and. abs(terr).lt.tsym/8.0) then
+           igoodsync=1
            ngoodsync=ngoodsync+1
            sqt=sqt + terr*terr
            sqf=sqf + ferr*ferr
+!        else
+!           write(*,3003) ferr,terr
+!3003       format('Sync failed:',f8.1,f8.3)
         endif
 
         a=0.
@@ -169,21 +179,24 @@ program sfoxtest
         f=1500.0
         call sfox_demod(crcvd,f,t,s3,chansym)    !Get s3 and hard symbol values
         call sym_prob(s3,rxdat,rxprob,rxdat2,rxprob2)
-        chansym(1:ND)=rxdat  !### TEMPORARY ? ###
+        if(igoodsync.eq.1) then
-!        do j=0,ND-1
+           do j=0,NN-1
-!        do j=0,5
+              if(chansym(1+j).ne.rxdat(j)) write(*,3001) xdt,j,chansym(1+j),  &
-!           write(*,3001) j,chansym(1+j),rxdat(j),rxprob(j),rxdat2(j),rxprob2(j)
+                   rxdat(j),rxprob(j),rxdat2(j),rxprob2(j)
-!3001       format('prob'i5,5i8)
+3001          format(f7.3,i5,5i8)
-!        enddo
+           enddo
+        endif
 
         nera=0
         chansym=mod(chansym,nq)                        !Enforce 0 to nq-1
         nharderr=count(chansym.ne.chansym0)            !Count hard errors
-!        nhard2=count(rxdat.ne.chansym0(1:ND))            !Count hard errors
+!        nhard2=count(rxdat.ne.chansym0(1:NN))            !Count hard errors
 !        print*,'A',nharderr,nhard2
         ntot=ntot+nharderr
         nworst=max(nworst,nharderr)
-        call rs_decode_sf(chansym,iera,nera,nfixed)    !Call the decoder
+        call rs_decode_sf(rxdat,iera,nera,nfixed)    !Call the BM decoder
+        ntrials=1000
+!        call ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials,nparam,correct,ntry)
         
         if(iand(nv,1).ne.0) then
            fname='000000_000001.wav'