Many tweaks to msk144 decoders. WSJT-X now calls the standalone fortran decoders.

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

CMakeLists.txt

@@ -367,8 +367,6 @@ set (wsjt_FSRCS
   lib/geniscat.f90
   lib/genmsk.f90
   lib/genmsk144.f90
-  lib/genmsk32.f90
-  lib/genmsk40.f90
   lib/genmsk40.f90
   lib/genmsk_short.f90
   lib/genqra64.f90
@@ -1038,12 +1036,6 @@ target_link_libraries (jt65sim wsjt_fort wsjt_cxx)
 add_executable (qra64sim lib/qra/qra64/qra64sim.f90 wsjtx.rc)
 target_link_libraries (qra64sim wsjt_fort wsjt_cxx)
 
-add_executable (msk32d lib/msk32d.f90 wsjtx.rc)
-target_link_libraries (msk32d wsjt_fort wsjt_cxx)
-
-add_executable (msk32d_ldpc lib/msk32d_ldpc.f90 wsjtx.rc)
-target_link_libraries (msk32d_ldpc wsjt_fort wsjt_cxx)
-
 add_executable (jt9sim lib/jt9sim.f90 wsjtx.rc)
 target_link_libraries (jt9sim wsjt_fort wsjt_cxx)
 
@@ -1220,7 +1212,7 @@ install (TARGETS udp_daemon message_aggregator
   BUNDLE DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime
   )
 
-install (TARGETS jt9 jt65code qra64code qra64sim jt9code jt4code wsprd msk32d msk32d_ldpc 
+install (TARGETS jt9 jt65code qra64code qra64sim jt9code jt4code wsprd 
   RUNTIME DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime
   BUNDLE DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime
   )

lib/bpdecode144.f90

@@ -1,3 +1,51 @@
+subroutine pltanh(x,y)
+  isign=+1
+  if( x.lt.0 ) then
+    isign=-1
+    z=abs(x)
+  endif
+  if( z.le. 0.8 ) then
+    y=0.83*x
+    return
+  elseif( z.le. 1.6 ) then
+    y=isign*(0.322*z+0.4064)
+    return  
+  elseif( z.le. 3.0 ) then
+    y=isign*(0.0524*z+0.8378)
+    return
+  elseif( z.lt. 7.0 ) then
+    y=isign*(0.0012*z+0.9914)
+    return
+  else
+    y=isign*0.9998
+    return
+  endif
+end subroutine pltanh
+
+subroutine platanh(x,y)
+  isign=+1
+  if( x.lt.0 ) then
+    isign=-1
+    z=abs(x)
+  endif
+  if( z.le. 0.664 ) then
+    y=x/0.83
+    return
+  elseif( z.le. 0.9217 ) then
+    y=isign*(z-0.4064)/0.322
+    return
+  elseif( z.le. 0.9951 ) then
+    y=isign*(z-0.8378)/0.0524
+    return
+  elseif( z.le. 0.9998 ) then
+    y=isign*(z-0.9914)/0.0012
+    return
+  else
+    y=isign*7.0
+    return
+  endif
+end subroutine platanh
+
 subroutine bpdecode144(llr,maxiterations,decoded,niterations)
 !
 ! A log-domain belief propagation decoder for the msk144 code.
@@ -11,13 +59,12 @@ integer*1 decoded(K)
 integer Nm(8,M)  ! 8 bits per check 
 integer Mn(3,N)  ! 3 checks per bit
 integer synd(M)
-real*8 tov(3,N)
+real tov(3,N)    ! single precision seems to be adequate in log-domain
-real*8 toc(8,M)
+real toc(8,M)
-real*8 tanhtoc(8,M)
+real tanhtoc(8,M)
-real*8 zn(N)
+real zn(N)
-real*8 llr(N)
+real llr(N)
-real*8 Tmn
+real Tmn
-real*8 xth
 
 data colorder/0,1,2,3,4,5,6,7,8,9, &
               10,11,12,13,14,15,24,26,29,30, &
@@ -239,12 +286,11 @@ do iter=0,maxiterations
   where( zn .gt. 0. ) cw=1
   ncheck=0
   do i=1,M
-    synd(i)=sum(cw(Nm(1:nrw,i)))
+    synd(i)=sum(cw(Nm(:,i)))
-    synd(i)=mod(synd(i),2)
+    if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
-    if( synd(i) .ne. 0 ) ncheck=ncheck+1
   enddo
 
-  if( ncheck .eq. 0 ) then ! we have a codeword - reorder the columns and return it.
+  if( ncheck .eq. 0 ) then ! we have a codeword
     niterations=iter
     codeword=cw(colorder+1)
     decoded=codeword(M+1:N)
@@ -254,10 +300,11 @@ do iter=0,maxiterations
 ! Send messages from bits to check nodes 
   do j=1,M
     do i=1,nrw
-      toc(i,j)=zn(Nm(i,j))  
+      ibj=Nm(i,j)
+      toc(i,j)=zn(ibj)  
       do kk=1,ncw ! subtract off what the bit had received from the check
-        if( Mn(kk,Nm(i,j)) .eq. j ) then  ! Mn(3,128)
+        if( Mn(kk,ibj) .eq. j ) then  ! Mn(3,128)
-          toc(i,j)=toc(i,j)-tov(kk,Nm(i,j))
+          toc(i,j)=toc(i,j)-tov(kk,ibj)
         endif
       enddo
     enddo
@@ -271,20 +318,12 @@ do iter=0,maxiterations
   do j=1,N
     do i=1,ncw
       ichk=Mn(i,j)  ! Mn(:,j) are the checks that include bit j
-      Tmn=1.0
+      Tmn=product(tanhtoc(:,ichk),mask=Nm(:,ichk).ne.j)
-      do kk=1,nrw
+      call platanh(-Tmn,y)
-        if( Nm(kk,ichk) .ne. j ) then
+      tov(i,j)=2*y
-          Tmn=Tmn*tanhtoc(kk,ichk)
-        endif
-      enddo 
-      tov(i,j)=2*atanh(-Tmn)
     enddo
   enddo
 
-  xth=35.0
-  where(tov .gt. xth) tov=xth
-  where(tov .lt. -xth) tov=-xth
-
 enddo
 niterations=-1
 end subroutine bpdecode144

lib/bpdecode40.f90

@@ -1,23 +1,22 @@
 subroutine bpdecode40(llr,maxiterations,decoded,niterations)
 !
 ! A log-domain belief propagation decoder for the msk40 code.
-! The code is a regular (32,16) code with column weight 3 and row weights 5,6,7.
+! The code is a regular (32,16) code with column weight 3, row weights 5,6,7.
 ! k9an August, 2016
 !
 integer, parameter:: N=32, K=16, M=N-K
 integer*1 codeword(N),cw(N)
 integer*1 colorder(N)
 integer*1 decoded(K)
-integer Nm(7,M)   
+integer Nm(7,M)  ! 5,6 or 7 bits per check 
 integer Mn(3,N)  ! 3 checks per bit
 integer synd(M)
-real*8 tov(3,N)
+real tov(3,N)
-real*8 toc(7,M)
+real toc(7,M)
-real*8 tanhtoc(7,M)
+real tanhtoc(7,M)
-real*8 zn(N)
+real zn(N)
-real*8 llr(N)
+real llr(N)
-real*8 Tmn
+real Tmn
-real*8 xth
 integer nrw(M)
 
 data colorder/ &
@@ -106,8 +105,7 @@ do iter=0,maxiterations
   ncheck=0
   do i=1,M
     synd(i)=sum(cw(Nm(1:nrw(i),i)))
-    synd(i)=mod(synd(i),2)
+    if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
-    if( synd(i) .ne. 0 ) ncheck=ncheck+1
   enddo
 
   if( ncheck .eq. 0 ) then ! we have a codeword - reorder the columns and return it
@@ -120,10 +118,11 @@ do iter=0,maxiterations
 ! Send messages from bits to check nodes 
   do j=1,M
     do i=1,nrw(j)
-      toc(i,j)=zn(Nm(i,j))  
+      ibj=Nm(i,j)
+      toc(i,j)=zn(ibj)  
       do kk=1,ncw ! subtract off what the bit had received from the check
-        if( Mn(kk,Nm(i,j)) .eq. j ) then  
+        if( Mn(kk,ibj) .eq. j ) then  
-          toc(i,j)=toc(i,j)-tov(kk,Nm(i,j))
+          toc(i,j)=toc(i,j)-tov(kk,ibj)
         endif
       enddo
     enddo
@@ -131,26 +130,18 @@ do iter=0,maxiterations
 
 ! send messages from check nodes to variable nodes
   do i=1,M
-      tanhtoc(1:nrw(i),i)=tanh(-toc(1:nrw(i),i)/2.)
+    tanhtoc(1:7,i)=tanh(-toc(1:7,i)/2)
   enddo
 
   do j=1,N
     do i=1,ncw
       ichk=Mn(i,j)  ! Mn(:,j) are the checks that include bit j
-      Tmn=1.0
+      Tmn=product(tanhtoc(1:nrw(i),ichk),mask=Nm(1:nrw(i),ichk).ne.j)
-      do kk=1,nrw(ichk)
+      call platanh(-Tmn,y)
-        if( Nm(kk,ichk) .ne. j ) then
+      tov(i,j)=y
-          Tmn=Tmn*tanhtoc(kk,ichk)
-        endif
-      enddo 
-      tov(i,j)=2.*atanh(-Tmn)
     enddo
   enddo
 
-  xth=35.0
-  where(tov .gt. xth) tov=xth
-  where(tov .lt. -xth) tov=-xth
-
 enddo
 niterations=-1
 return

lib/detectmsk144.f90

@@ -32,7 +32,7 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
   real detfer(MAXSTEPS)
   real rcw(12)
   real dd(NPTS)
-  real ddr(NPTS)
+!  real ddr(NPTS)
   real ferrs(MAXCAND)
   real pp(12)                          !Half-sine pulse shape
   real snrs(MAXCAND)
@@ -41,6 +41,7 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
   real*8 dt, df, fs, pi, twopi
   real softbits(144)
   real*8 lratio(128)
+  real llr(128)
   logical first
   data first/.true./
   data s8/0,1,1,1,0,0,1,0/
@@ -49,7 +50,7 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
 
   i=index(pchk_file,".pchk")
   gen_file=pchk_file(1:i-1)//".gen"
-  call init_ldpc(trim(pchk_file)//char(0),trim(gen_file)//char(0))
+!  call init_ldpc(trim(pchk_file)//char(0),trim(gen_file)//char(0))
 
   if(first) then
      nmatchedfilter=1
@@ -186,6 +187,7 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
   nmessages=0
   allmessages=char(0)
   lines=char(0)
+  nshort=0
 
   do ip=1,ndet  ! Try to sync/demod/decode each candidate.
     imid=times(ip)*fs
@@ -215,16 +217,17 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
     enddo
     cc=cc1+cc2
     dd=abs(cc1)*abs(cc2)
-    do i=1,NPTS-(40*6+41)
+!    do i=1,NPTS-(40*6+41)
-      ccr1(i)=sum(cdat(i:i+41)*conjg(cbr))
+!      ccr1(i)=sum(cdat(i:i+41)*conjg(cbr))
-      ccr2(i)=sum(cdat(i+40*6:i+40*6+41)*conjg(cbr))
+!      ccr2(i)=sum(cdat(i+40*6:i+40*6+41)*conjg(cbr))
-    enddo
+!    enddo
-    ccr=ccr1+ccr2
+!    ccr=ccr1+ccr2
-    ddr=abs(ccr1)*abs(ccr2)
+!    ddr=abs(ccr1)*abs(ccr2)
     cmax=maxval(abs(cc))
-    crmax=maxval(abs(ccr))
+!    crmax=maxval(abs(ccr))
-    ishort=0
+!    if( crmax .gt. cmax ) then
-    if( crmax .gt. cmax ) ishort=1
+!      nshort=nshort+1
+!    endif
  
 ! Find 6 largest peaks
     do ipk=1, 6
@@ -374,13 +377,15 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
               sigma=0.75
               lratio(1:48)=softbits(9:9+47)
               lratio(49:128)=softbits(65:65+80-1)
+              llr=2.0*lratio/(sigma*sigma)
               lratio=exp(2.0*lratio/(sigma*sigma))
   
               max_iterations=10
               max_dither=1
               call timer('ldpcdecod',0)
-              call ldpc_decode(lratio, decoded, &
+!              call ldpc_decode(lratio, decoded, &
-                           max_iterations, niterations, max_dither, ndither)
+!                           max_iterations, niterations, max_dither, ndither)
+              call bpdecode144(llr,max_iterations,decoded,niterations)
               call timer('ldpcdecod',1)
 
               if( niterations .ge. 0.0 ) then

lib/detectmsk40.f90

@@ -37,6 +37,7 @@ subroutine detectmsk40(cbig,n,pchk_file,mycall,hiscall,lines,nmessages,   &
   real rcw(12)
   real ddr(NPTS)
   real ferrs(MAXCAND)
+  real llr(32)
   real pp(12)                          !Half-sine pulse shape
   real snrs(MAXCAND)
   real softbits(40)
@@ -112,7 +113,7 @@ subroutine detectmsk40(cbig,n,pchk_file,mycall,hiscall,lines,nmessages,   &
   pchk_file40=pchk_file(1:i-1)//"32-16"//pchk_file(i+6:)
   i=index(pchk_file40,".pchk")
   gen_file40=pchk_file40(1:i-1)//".gen"
-  call init_ldpc(trim(pchk_file40)//char(0),trim(gen_file40)//char(0))
+!  call init_ldpc(trim(pchk_file40)//char(0),trim(gen_file40)//char(0))
  
 ! Fill the detmet, detferr arrays
   nstepsize=60  ! 5ms steps
@@ -378,16 +379,19 @@ subroutine detectmsk40(cbig,n,pchk_file,mycall,hiscall,lines,nmessages,   &
 
               sigma=0.75
               if(xsnr.lt.1.5) sigma=1.1 - 0.0875*(xsnr+4.0) 
-              lratio(1:32)=exp(2.0*softbits(9:40)/(sigma*sigma))
+              lratio(1:32)=exp(2.0*softbits(9:40)/(sigma*sigma)) ! Use this for Radford Neal's routines
+              llr(1:32)=2.0*softbits(9:40)/(sigma*sigma)  ! Use log likelihood for bpdecode40
 
               max_iterations=5
               max_dither=1
-              call ldpc_decode(lratio,decoded,max_iterations,niterations,max_dither,ndither)
+!              call ldpc_decode(lratio,decoded,max_iterations,niterations,max_dither,ndither)
+              call bpdecode40(llr,max_iterations, decoded, niterations)
               ncalls=ncalls+1
                
               nhashflag=0
               if( niterations .ge. 0 ) then
-                call ldpc_encode(decoded,cw)
+                call encode_msk40(decoded,cw)
+!                call ldpc_encode(decoded,cw)
                 nhammd=0
                 cord=0.0
                 do i=1,32

lib/msk144d.f90

@@ -18,34 +18,38 @@ program msk144d
   character*12 mycall,hiscall
   character(len=500) optarg
 
-  type (option) :: long_options(5) = [ &
+  type (option) :: long_options(6) = [ &
-       option ('help',.false.,'h','Display this help message',''), &
+       option ('dxcall',.true.,'d','hiscall',''), &  
-       option ('mycall',.true.,'c','mycall',''), &
        option ('evemode',.true.,'e','',''), &
+       option ('help',.false.,'h','Display this help message',''), &
+       option ('mycall',.true.,'m','mycall',''), &
        option ('nftol',.true.,'n','nftol',''), &
-       option ('hiscall',.true.,'x','hiscall','') &  
+       option ('short',.false.,'s','enable Sh','') &
        ]
   t0=0.0
   ntol=100
   mycall=''
   hiscall=''
+  bShMsgs=.false.
  
   do
-     call getopt('c:ehn:x:',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.)
+     call getopt('d:ehm:ns:',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.)
      if( nstat .ne. 0 ) then
         exit
      end if
      select case (c)
-     case ('h')
+     case ('d')
-        display_help = .true.
+        read (optarg(:narglen), *) hiscall
-     case ('c')
-        read (optarg(:narglen), *) mycall
      case ('e')
         t0=1e-4
+     case ('h')
+        display_help = .true.
+     case ('m')
+        read (optarg(:narglen), *) mycall
      case ('n')
         read (optarg(:narglen), *) ntol
-     case ('x')
+     case ('s')
-        read (optarg(:narglen), *) hiscall
+        bShMsgs=.true. 
      end select
   end do
 
@@ -56,7 +60,6 @@ program msk144d
      print *, '       msk144 decode pre-recorded .WAV file(s)'
      print *, ''
      print *, 'OPTIONS:'
-     print *, ''
      do i = 1, size (long_options)
         call long_options(i) % print (6)
      end do
@@ -65,7 +68,6 @@ program msk144d
 
   call init_timer ('timer.out')
   call timer('msk144  ',0)
-  bShMsgs=.true.
   pchk_file='./peg-128-80-reg3.pchk'
   ndecoded=0
   do ifile=noffset+1,noffset+nremain