1. add error-pattern rejection to osd174 (factor of 3 speedup for norder=3)

2. osd codeword selected based on correlation discrepancy instead of Hamming dist (improves sensitivity)
3. trap zero-valued spectrum values caused by large dt.
4. trap all-zero codeword
4. both decoders return number of hard errors
5. osd174 returns correlation discrepancy (rejection threshold still needs to be tuned).
6. both decoders return the codeword with bits in as-transmitted order, for subtraction.


git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7735 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/fsk4hf/bpdecode174.f90

@@ -1,4 +1,4 @@
-subroutine bpdecode174(llr,apmask,maxiterations,decoded,niterations)
+subroutine bpdecode174(llr,apmask,maxiterations,decoded,cw,nharderror)
 !
 ! A log-domain belief propagation decoder for the (174,87) code.
 !
@@ -306,6 +306,7 @@ data nrw/   &
 
 ncw=3
 
+decoded=0
 toc=0
 tov=0
 tanhtoc=0
@@ -340,14 +341,13 @@ do iter=0,maxiterations
   enddo
 ! write(*,*) 'number of unsatisfied parity checks ',ncheck
   if( ncheck .eq. 0 ) then ! we have a codeword - reorder the columns and return it
-!    niterations=iter
     codeword=cw(colorder+1)
     decoded=codeword(M+1:N)
     nerr=0
     do i=1,N
       if( (2*cw(i)-1)*llr(i) .lt. 0.0 ) nerr=nerr+1
     enddo
-    niterations=nerr
+    nharderror=nerr
     return
   endif
 
@@ -361,7 +361,7 @@ do iter=0,maxiterations
     endif
 !    write(*,*) iter,ncheck,nd,ncnt
     if( ncnt .ge. 5 .and. iter .ge. 10 .and. ncheck .gt. 15) then
-      niterations=-1
+      nharderror=-1
       return
     endif
   endif
@@ -396,6 +396,6 @@ do iter=0,maxiterations
   enddo
 
 enddo
-niterations=-1
+nharderror=-1
 return
 end subroutine bpdecode174

lib/fsk4hf/ft8b.f90

@@ -35,10 +35,9 @@ subroutine ft8b(datetime,s,candidate,ncand)
         j=j+1
         s1(0:7,j)=s(ia:ib:2,n)
      enddo
-     
      do j=1,ND
         ps=s1(0:7,j)
-        ps=log(ps)
+        where (ps.gt.0.0) ps=log(ps)
         r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6))
         r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5))
         r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3))
@@ -46,28 +45,44 @@ subroutine ft8b(datetime,s,candidate,ncand)
         rxdata(3*j-1)=r2
         rxdata(3*j)=r1
      enddo
-     rxav=sum(rxdata)/ND
+
-     rx2av=sum(rxdata*rxdata)/ND
+     rxav=sum(rxdata)/(3.0*ND)
-     rxsig=sqrt(rx2av-rxav*rxav)
+     rx2av=sum(rxdata*rxdata)/(3.0*ND)
+     var=rx2av-rxav*rxav
+     if( var .gt. 0.0 ) then
+       rxsig=sqrt(var)
+     else
+       rxsig=sqrt(rx2av)
+     endif
      rxdata=rxdata/rxsig
      ss=0.84
      llr=2.0*rxdata/(ss*ss)
      apmask=0
-     call bpdecode174(llr,apmask,max_iterations,decoded,niterations)
+     cw=0
-     if(niterations.lt.0) call osd174(llr,norder,decoded,nharderrors,cw)
+! cw will be needed for subtraction.
-     nbadcrc=0
+! dmin is the correlation discrepancy of a returned codeword - it is 
-     call chkcrc12a(decoded,nbadcrc)
+!      used to select the best codeword within osd174.
-
+     call bpdecode174(llr,apmask,max_iterations,decoded,cw,nharderrors)
+     dmin=0.0
+     if(nharderrors.lt.0) then
+       call osd174(llr,norder,decoded,cw,nharderrors,dmin)
+! This threshold needs to be tuned. 99.0 should pass everything.
+       if( dmin .gt. 99.0 ) nharderrors=-1
+     endif
+! Reject the all-zero codeword
+     if( count(cw.eq.0) .eq. 174 ) cycle
+     nbadcrc=1
+     if( nharderrors .ge. 0 ) call chkcrc12a(decoded,nbadcrc)
      message='                      '
      if(nbadcrc.eq.0) then
         call extractmessage174(decoded,message,ncrcflag,recent_calls,nrecent)
         nsnr=nint(10.0*log10(sync) - 25.5)    !### empirical ###
-        write(13,1110) datetime,0,nsnr,xdt,f1,xdta,f1a,niterations,    &
-             nharderrors,message
-1110    format(a13,2i4,2(f6.2,f7.1),2i4,2x,a22)
         write(*,1112) datetime(8:13),nsnr,xdt,nint(f1),message
 1112    format(a6,i4,f5.1,i6,2x,a22)
      endif
+     write(13,1110) datetime,0,nsnr,xdt,f1,xdta,f1a,    &
+          nharderrors,dmin,message
+1110 format(a13,2i4,2(f6.2,f7.1),i4,2x,f6.2,2x,a22)
   enddo
 
   return

lib/fsk4hf/ldpcsim174.f90

@@ -188,12 +188,10 @@ do idb = 20,-10,-1
     apmask(colorder(174-87+1:174-87+nap)+1)=1
 
 ! max_iterations is max number of belief propagation iterations
-    call bpdecode174(llr, apmask, max_iterations, decoded, niterations)
+    call bpdecode174(llr, apmask, max_iterations, decoded, cw, nharderrors)
-ni1=niterations
+    if( norder .ge. 0 .and. nharderrors .lt. 0 ) call osd174(llr, norder, decoded, cw,  nharderrors)
-    if( norder .ge. 0 .and. niterations .lt. 0 ) call osd174(llr, norder, decoded, niterations, cw)
+! If the decoder finds a valid codeword, nharderrors will be .ge. 0.
-ni2=niterations
+    if( nharderrors .ge. 0 ) then
-! If the decoder finds a valid codeword, niterations will be .ge. 0.
-    if( niterations .ge. 0 ) then
       call extractmessage174(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
       if( ncrcflag .ne. 1 ) then
         nbadcrc=nbadcrc+1

lib/fsk4hf/osd174.f90

@@ -1,18 +1,14 @@
-subroutine osd174(llr,norder,decoded,niterations,cw)
+subroutine osd174(llr,norder,decoded,cw,nhardmin,dmin)
 !
-! An ordered-statistics decoder based on ideas from: 
+! An ordered-statistics decoder for the (174,87) code.
-! "Soft-decision decoding of linear block codes based on ordered statistics,"
-! by Marc P. C. Fossorier and Shu Lin, 
-! IEEE Trans Inf Theory, Vol 41, No 5, Sep 1995
 ! 
-
 include "ldpc_174_87_params.f90"
 
 integer*1 gen(K,N)
-integer*1 genmrb(K,N)
+integer*1 genmrb(K,N),g2(N,K)
-integer*1 temp(K),m0(K),me(K)
+integer*1 temp(K),m0(K),me(K),mi(K)
-integer indices(N)
+integer indices(N),nxor(N)
-integer*1 codeword(N),cw(N),hdec(N)
+integer*1 cw(N),ce(N),c0(N),hdec(N)
 integer*1 decoded(K)
 integer indx(N)
 real llr(N),rx(N),absrx(N)
@@ -28,9 +24,7 @@ if( first ) then ! fill the generator matrix
       read(g(i)(j:j),"(Z1)") istr
         do jj=1, 4 
           irow=(j-1)*4+jj
-          if( irow .le. K ) then
+          if( btest(istr,4-jj) ) gen(irow,i)=1
-            if( btest(istr,4-jj) ) gen(irow,i)=1
-          endif
         enddo
     enddo
   enddo
@@ -50,36 +44,33 @@ where(rx .ge. 0) hdec=1
 ! use magnitude of received symbols as a measure of reliability.
 absrx=abs(rx) 
 call indexx(absrx,N,indx)  
-! re-order the columns of the generator matrix in order of increasing reliability.
+
+! re-order the columns of the generator matrix in order of decreasing reliability.
 do i=1,N
-  genmrb(1:K,N+1-i)=gen(1:K,indx(N+1-i))
+  genmrb(1:K,i)=gen(1:K,indx(N+1-i))
+  indices(i)=indx(N+1-i)
 enddo
 
 ! do gaussian elimination to create a generator matrix with the most reliable
-! received bits as the systematic bits. if it happens that the K most reliable
+! received bits in positions 1:K in order of decreasing reliability (more or less). 
-! bits are not independent, then we dip into the bits just below the K best bits
+! reliability will not be strictly decreasing because column re-ordering is needed
-! to find K independent most reliable bits. the "indices" array tracks column 
+! to put the generator matrix in systematic form. the "indices" array tracks 
-! permutations caused by reliability sorting and gaussian elimination.
+! column permutations caused by reliability sorting and gaussian elimination.
-do i=1,N
-  indices(i)=indx(i)
-enddo
 do id=1,K ! diagonal element indices 
-  do ic=id,K+20  ! The 20 is ad hoc - beware
+  do icol=id,K+20  ! The 20 is ad hoc - beware
-    icol=N-K+ic
-    if( icol .gt. N ) icol=M+1-(icol-N)
     iflag=0
     if( genmrb(id,icol) .eq. 1 ) then
       iflag=1
-      if( icol-M .ne. id ) then ! reorder column
+      if( icol .ne. id ) then ! reorder column
-        temp(1:K)=genmrb(1:K,M+id)
+        temp(1:K)=genmrb(1:K,id)
-        genmrb(1:K,M+id)=genmrb(1:K,icol)
+        genmrb(1:K,id)=genmrb(1:K,icol)
         genmrb(1:K,icol)=temp(1:K) 
-        itmp=indices(M+id)
+        itmp=indices(id)
-        indices(M+id)=indices(icol)
+        indices(id)=indices(icol)
         indices(icol)=itmp
       endif
       do ii=1,K
-        if( ii .ne. id .and. genmrb(ii,N-K+id) .eq. 1 ) then
+        if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then
           genmrb(ii,1:N)=mod(genmrb(ii,1:N)+genmrb(id,1:N),2)
         endif
       enddo
@@ -88,64 +79,102 @@ do id=1,K ! diagonal element indices
   enddo
 enddo
 
-! use the hard decisions for the K MRB bits to define the order 0 
+g2=transpose(genmrb)
-! message, m0. Encode m0 using the modified generator matrix to 
-! find the order 0 codeword. Flip all combinations of N bits in m0
-! and re-encode to find the list of order N codewords. Test all such
-! codewords against the received word to decide which codeword is 
-! most likely to be correct.
-m0=0
-where (rx(indices(M+1:N)).ge.0.0) m0=1
 
-nhardmin=N
+! The hard decisions for the K MRB bits define the order 0 message, m0. 
-corrmax=-1.0e32
+! Encode m0 using the modified generator matrix to find the "order 0" codeword. 
-j0=0
+! Flip various combinations of bits in m0 and re-encode to generate a list of
-j1=0
+! codewords. Test all such codewords against the received word to decide which
-j2=0
+! codeword is most likely to be correct.
-j3=0
+
-if( norder.ge.4 ) j0=K
+hdec=hdec(indices)   ! hard decisions from received symbols
-if( norder.ge.3 ) j1=K
+m0=hdec(1:K)         ! zero'th order message
-if( norder.ge.2 ) j2=K
+absrx=absrx(indices) 
-if( norder.ge.1 ) j3=K
+rx=rx(indices)       
+
+s1=sum(absrx(1:K))
+s2=sum(absrx(K+1:N))
+xlam=7.0  ! larger values reject more error patterns 
+rho=s1/(s1+xlam*s2)
+call mrbencode(m0,c0,g2,N,K)
+nxor=ieor(c0,hdec)
+nhardmin=sum(nxor)
+dmin=sum(nxor*absrx)
+thresh=rho*dmin
+
+cw=c0
 nt=0
-do i1=0,j0
+nrejected=0
-  do i2=i1,j1
+do iorder=1,norder
-    do i3=i2,j2
+  mi(1:K-iorder)=0
-      do i4=i3,j3
+  mi(K-iorder+1:K)=1
-nt=nt+1
+  iflag=0
-        me=m0
+  do while(iflag .ge. 0 ) 
-        if( i1 .ne. 0 ) me(i1)=1-me(i1)
+    dpat=sum(mi*absrx(1:K))
-        if( i2 .ne. 0 ) me(i2)=1-me(i2)
+    nt=nt+1
-        if( i3 .ne. 0 ) me(i3)=1-me(i3)
+    if( dpat .lt. thresh ) then  ! reject unlikely error patterns
-        if( i4 .ne. 0 ) me(i4)=1-me(i4)
+      me=ieor(m0,mi)
-
+      call mrbencode(me,ce,g2,N,K)
-! me is the m0 + error pattern. encode this message using genmrb to
+      nxor=ieor(ce,hdec)
-! produce a codeword. test the codeword against the received vector
+      dd=sum(nxor*absrx)
-! and save it if it's the best that we've seen so far. 
+      if( dd .lt. dmin ) then
-        do i=1,N 
+        dmin=dd
-          nsum=sum(iand(me,genmrb(1:K,i)))
+        cw=ce
-          codeword(i)=mod(nsum,2)
+        nhardmin=sum(nxor)
-        enddo
+        thresh=rho*dmin
-! undo the bit re-ordering to put the "real" message bits at the end
+      endif
-        codeword(indices)=codeword
+    else
-        nhard=count(codeword .ne. hdec)
+      nrejected=nrejected+1
-!        corr=sum(codeword*rx)  ! to save time use nhard to pick best codeword
+    endif
-        if( nhard .lt. nhardmin ) then
+! get the next test error pattern, iflag will go negative
-!         if( corr .gt. corrmax ) then
+! when the last pattern with weight iorder has been generated
-          cw=codeword
+    call nextpat(mi,k,iorder,iflag)
-          nhardmin=nhard
+  enddo
-!          corrmax=corr
-          i1min=i1
-          i2min=i2
-          i3min=i3
-          i4min=i4
-          if( nhardmin .le. 15 ) goto 200 ! early exit - tune for each code 
-        endif
-      enddo
-    enddo
-  enddo 
 enddo
-200 decoded=cw(M+1:N)
+
-niterations=nhardmin
+!write(*,*) 'nhardmin ',nhardmin
+!write(*,*) 'total patterns ',nt,' number rejected ',nrejected
+
+! re-order the codeword to place message bits at the end
+cw(indices)=cw
+hdec(indices)=hdec
+decoded=cw(M+1:N)
+cw(colorder+1)=cw ! put the codeword back into received-word order
 return
 end subroutine osd174
+
+subroutine mrbencode(me,codeword,g2,N,K)
+integer*1 me(K),codeword(N),g2(N,K)
+! fast encoding for low-weight test patterns
+  codeword=0
+  do i=1,K
+    if( me(i) .eq. 1 ) then
+      codeword=ieor(codeword,g2(1:N,i))
+    endif
+  enddo
+return
+end subroutine mrbencode
+
+subroutine nextpat(mi,k,iorder,iflag)
+  integer*1 mi(k),ms(k)
+! generate the next test error pattern
+  ind=-1
+  do i=1,k-1
+     if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i 
+  enddo
+  if( ind .lt. 0 ) then ! no more patterns of this order
+    iflag=ind
+    return
+  endif
+  ms=0
+  ms(1:ind-1)=mi(1:ind-1)
+  ms(ind)=1
+  ms(ind+1)=0
+  if( ind+1 .lt. k ) then
+     nz=iorder-sum(ms)
+     ms(k-nz+1:k)=1
+  endif
+  mi=ms
+  iflag=ind
+  return
+end subroutine nextpat