Soft-decision RS decoding now basically working. Needs better tuning, no doubt.

CMakeLists.txt

@@ -593,6 +593,7 @@ set (wsjt_FSRCS
   lib/superfox/sfox_clo.f90
   lib/superfox/sym_prob.f90
   lib/superfox/getpp3.f90
+  lib/superfox/ftrsd3.f90
   lib/superfox/ran1.f90
   )
 
@@ -638,7 +639,6 @@ 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

@@ -1,4 +1,5 @@
-subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
+subroutine ftrsd3(s3,chansym0,rxdat,rxprob,rxdat2,rxprob2,ntrials0,  &
+     correct,param,ntry)
 
 ! Soft-decision decoder for Reed-Solomon codes.
  
@@ -14,8 +15,8 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
 
   use sfox_mod
 
-  integer, dimension(0:NN-1) :: rxdat,rxprob,rxdat2,rxprob2,workdat,   &
+  real s3(0:NQ-1,0:NN-1)          !Symbol spectra
-       correct,indexes
+  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
@@ -26,8 +27,9 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
   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
+  integer pass,tmp,thresh
   
   integer perr(0:7,0:7)
   data perr/ 4, 9,11,13,14,14,15,15, &
@@ -57,12 +59,13 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
   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
-        enddo
+        endif
      enddo
   enddo
 
@@ -83,7 +86,8 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
      param(5)=0
      param(7)=1000*1000               !???
      ntry=0
-     return
+!     print*,'AA1',nerr
+     go to 900
   endif
 
 ! Hard-decision decoding failed.  Try the FT soft-decision method.
@@ -101,10 +105,10 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
   do i=0,NN-1
      nsum=nsum+rxprob(i)
      j=indexes(NN-1-i)
-     ratio=(float)rxprob2(j)/((float)rxprob(j)+0.01)
+     ratio=float(rxprob2(j))/(float(rxprob(j))+0.01)
      ii=7.999*ratio
-     jj=(NN-1-i)/8
+     jj=int((7.999/NN)*(NN-1-i))
-     thresh0(i)=1.3*perr(jj)(ii)
+     thresh0(i)=1.3*perr(jj,ii)
   enddo
   if(nsum.le.0) return
 
@@ -118,30 +122,35 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
 ! 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(126-i)
+        j=indexes(NN-1-i)
         thresh=thresh0(i)
-
 ! Generate a random number ir, 0 <= ir <= 100 (see POSIX.1-2001 example).
-        nseed=nseed*1103515245 + 12345
+!        nseed=nseed*1103515245 + 12345
-        ir=mod(nseed/65536),32768)
+!        ir=mod(nseed/65536,32768)
-        ir=(100*ir)/32768
+!        ir=(100*ir)/32768
-        nseed=iand(ir,4294967295)
+!        nseed=iand(ir,2147483647)
 
-        if((ir.lt.thresh ) .and. numera.lt.(NN-KK)) then
+        ir=100.0*ran1(nseed)
+        if((ir.lt.thresh) .and. numera.lt.(NN-KK)) then
            era_pos(numera)=j
            numera=numera+1
+           if(rxdat(j).ne.chansym0(j)) then
+              ncaught=ncaught+1
+           endif
         endif
      enddo
-
-!     nerr=decode_rs_int(rs,workdat,era_pos,numera,0);
      call rs_decode_sf(workdat,era_pos,numera,nerr)    !Call the decoder
-
+     do i=0,NN-1
+        write(60,3101) i,chansym0(i),workdat(i),workdat(i)-chansym0(i)
+3101    format(4i8)
+     enddo
      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.
+      ! s3(NQ,NN) of synchronized symbol spectra.
         ncandidates=ncandidates+1
         nhard=0
         nsoft=0
@@ -151,7 +160,7 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
               if(workdat(i) .ne. rxdat2(i)) nsoft=nsoft+rxprob(i)
            endif
         enddo
-        nsoft=(NN-1)*nsoft/nsum
+        nsoft=NN*nsoft/nsum
         ntotal=nsoft+nhard
 
         pp=0.
@@ -168,8 +177,8 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
         else
            if(pp.gt.pp2 .and. pp.ne.pp1) pp2=pp
         endif
-        if(nhard_min <= 41 && ntotal_min <= 71) exit   !### New values ###
+        if(nhard_min.le.60 .and. ntotal_min.le.90) exit   !### Needs tuning
-     enddo
+     endif
      if(k.eq.ntrials) ntry=k
   enddo
 
@@ -177,12 +186,13 @@ subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry)
   param(1)=nhard_min
   param(2)=nsoft_min
   param(3)=nera_best
-!  param(4)= pp1 > 0 ? 1000.0*pp2/pp1 : 1000.0
+  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
 
-  return
+900 return
 end subroutine ftrsd3

lib/superfox/sfox_demod.f90

@@ -33,5 +33,8 @@ subroutine sfox_demod(crcvd,f,t,s3,chansym)
      chansym(n)=ipk(1) - 1
   enddo
 
+  call pctile(s3,NQ*NN,50,base)
+  s3=s3/base
+
   return
 end subroutine sfox_demod

lib/superfox/sfoxtest.f90

@@ -6,7 +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)
+  integer param(0:8)
   real*4 xnoise(NMAX)                    !Random noise
   real*4 dat(NMAX)                       !Generated real data
   complex cdat(NMAX)                     !Generated complex waveform
@@ -81,9 +81,9 @@ program sfoxtest
   allocate(s3(0:NQ-1,0:NN-1))
   allocate(msg0(1:KK))
   allocate(parsym(1:NN-KK))
-  allocate(chansym0(1:NN))
+  allocate(chansym0(0:NN-1))
-  allocate(chansym(1:NN))
+  allocate(chansym(0:NN-1))
-  allocate(iera(1:NN))
+  allocate(iera(0:NN-1))
   allocate(rxdat(0:NN-1))
   allocate(rxprob(0:NN-1))
   allocate(rxdat2(0:NN-1))
@@ -106,8 +106,8 @@ program sfoxtest
 
   call rs_init_sf(MM,NQ,NN,KK,NFZ)          !Initialize the Karn codec
   call rs_encode_sf(msg0,parsym)            !Compute parity symbols
-  chansym0(1:kk)=msg0(1:kk)
+  chansym0(0:kk-1)=msg0(1:kk)
-  chansym0(kk+1:nn)=parsym(1:nn-kk)
+  chansym0(kk:nn-1)=parsym(1:nn-kk)
 
 ! Generate clo, the LO for sync detection
   call sfox_clo(fsample,syncwidth,clo)
@@ -143,7 +143,6 @@ 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.3*ran1(idummy)
            call sfox_gen(chansym0,f1,fsample,syncwidth,cdat)
         endif
@@ -169,10 +168,7 @@ program sfoxtest
            ngoodsync=ngoodsync+1
            sqt=sqt + terr*terr
            sqf=sqf + ferr*ferr
-!        else
+       endif
-!           write(*,3003) ferr,terr
-!3003       format('Sync failed:',f8.1,f8.3)
-        endif
 
         a=0.
         a(1)=1500.0-f
@@ -180,42 +176,28 @@ 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)
-        if(igoodsync.eq.1) then
-           do j=0,NN-1
-              if(chansym(1+j).ne.rxdat(j)) write(*,3001) xdt,j,chansym(1+j),  &
-                   rxdat(j),rxprob(j),rxdat2(j),rxprob2(j)
-3001          format(f7.3,i5,5i8)
-           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:NN))            !Count hard errors
-!        print*,'A',nharderr,nhard2
         ntot=ntot+nharderr
         nworst=max(nworst,nharderr)
-        call rs_decode_sf(rxdat,iera,nera,nfixed)    !Call the BM 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)
+        call ftrsd3(s3,chansym0,rxdat,rxprob,rxdat2,rxprob2,ntrials,  &
-        
+             correct,param,ntry)
+
         if(iand(nv,1).ne.0) then
            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)
-!           write(*,1100) f1,xdt
+       endif
-!1100       format(/'f0:',f7.1,'  xdt:',f6.2)
-!           write(*,1112) f,t
-!1112       format('f: ',f7.1,'   DT:',f6.2)
-!           write(*,1110) ferr,terr
-!1110       format('err:',f6.1,f12.2)
-!           write(*,1120) nharderr
-!1120       format('Hard errors:',i4)
-        endif
 
-        if(nharderr.le.maxerr) ngood=ngood+1
+!        if(nharderr.le.maxerr) ngood=ngood+1
+       if(count(correct.ne.chansym0).eq.0) ngood=ngood+1
      enddo  ! ifile
      fgoodsync=float(ngoodsync)/nfiles
      fgood=float(ngood)/nfiles