diff --git a/CMakeLists.txt b/CMakeLists.txt index 12b25c7ef..415be33b4 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/lib/bpdecode144.f90 b/lib/bpdecode144.f90 index e1f5b5d2e..b67e5f944 100644 --- a/lib/bpdecode144.f90 +++ b/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*8 toc(8,M) -real*8 tanhtoc(8,M) -real*8 zn(N) -real*8 llr(N) -real*8 Tmn -real*8 xth +real tov(3,N) ! single precision seems to be adequate in log-domain +real toc(8,M) +real tanhtoc(8,M) +real zn(N) +real llr(N) +real Tmn 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)=mod(synd(i),2) - if( synd(i) .ne. 0 ) ncheck=ncheck+1 + synd(i)=sum(cw(Nm(:,i))) + if( mod(synd(i),2) .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) - toc(i,j)=toc(i,j)-tov(kk,Nm(i,j)) + if( Mn(kk,ibj) .eq. j ) then ! Mn(3,128) + toc(i,j)=toc(i,j)-tov(kk,ibj) endif enddo enddo @@ -265,26 +312,18 @@ do iter=0,maxiterations ! send messages from check nodes to variable nodes do i=1,M - tanhtoc(1:nrw,i)=tanh(-toc(1:nrw,i)/2) + tanhtoc(1:nrw,i)=tanh(-toc(1:nrw,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 - do kk=1,nrw - if( Nm(kk,ichk) .ne. j ) then - Tmn=Tmn*tanhtoc(kk,ichk) - endif - enddo - tov(i,j)=2*atanh(-Tmn) + Tmn=product(tanhtoc(:,ichk),mask=Nm(:,ichk).ne.j) + call platanh(-Tmn,y) + tov(i,j)=2*y enddo enddo - xth=35.0 - where(tov .gt. xth) tov=xth - where(tov .lt. -xth) tov=-xth - enddo niterations=-1 end subroutine bpdecode144 diff --git a/lib/bpdecode40.f90 b/lib/bpdecode40.f90 index 73eaf0a0c..7839980b3 100644 --- a/lib/bpdecode40.f90 +++ b/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*8 toc(7,M) -real*8 tanhtoc(7,M) -real*8 zn(N) -real*8 llr(N) -real*8 Tmn -real*8 xth +real tov(3,N) +real toc(7,M) +real tanhtoc(7,M) +real zn(N) +real llr(N) +real Tmn 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( synd(i) .ne. 0 ) ncheck=ncheck+1 + if( mod(synd(i),2) .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 - toc(i,j)=toc(i,j)-tov(kk,Nm(i,j)) + if( Mn(kk,ibj) .eq. j ) then + 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 - do kk=1,nrw(ichk) - if( Nm(kk,ichk) .ne. j ) then - Tmn=Tmn*tanhtoc(kk,ichk) - endif - enddo - tov(i,j)=2.*atanh(-Tmn) + Tmn=product(tanhtoc(1:nrw(i),ichk),mask=Nm(1:nrw(i),ichk).ne.j) + call platanh(-Tmn,y) + tov(i,j)=y enddo enddo - xth=35.0 - where(tov .gt. xth) tov=xth - where(tov .lt. -xth) tov=-xth - enddo niterations=-1 return diff --git a/lib/detectmsk144.f90 b/lib/detectmsk144.f90 index 8e004fcae..12c0ab2c3 100644 --- a/lib/detectmsk144.f90 +++ b/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 @@ -200,7 +202,7 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00) ! remove coarse freq error - should now be within a few Hz call tweak1(cdat,NPTS,-(1500+ferr),cdat) - + ! attempt frame synchronization ! correlate with sync word waveforms cc=0 @@ -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) - ccr1(i)=sum(cdat(i:i+41)*conjg(cbr)) - ccr2(i)=sum(cdat(i+40*6:i+40*6+41)*conjg(cbr)) - enddo - ccr=ccr1+ccr2 - ddr=abs(ccr1)*abs(ccr2) +! do i=1,NPTS-(40*6+41) +! ccr1(i)=sum(cdat(i:i+41)*conjg(cbr)) +! ccr2(i)=sum(cdat(i+40*6:i+40*6+41)*conjg(cbr)) +! enddo +! ccr=ccr1+ccr2 +! ddr=abs(ccr1)*abs(ccr2) cmax=maxval(abs(cc)) - crmax=maxval(abs(ccr)) - ishort=0 - if( crmax .gt. cmax ) ishort=1 +! crmax=maxval(abs(ccr)) +! if( crmax .gt. cmax ) then +! 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, & - max_iterations, niterations, max_dither, ndither) +! call ldpc_decode(lratio, decoded, & +! max_iterations, niterations, max_dither, ndither) + call bpdecode144(llr,max_iterations,decoded,niterations) call timer('ldpcdecod',1) if( niterations .ge. 0.0 ) then diff --git a/lib/detectmsk40.f90 b/lib/detectmsk40.f90 index 599d08922..3920c3a0a 100644 --- a/lib/detectmsk40.f90 +++ b/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 diff --git a/lib/msk144d.f90 b/lib/msk144d.f90 index 9a886ae2c..846b81b47 100644 --- a/lib/msk144d.f90 +++ b/lib/msk144d.f90 @@ -18,34 +18,38 @@ program msk144d character*12 mycall,hiscall character(len=500) optarg - type (option) :: long_options(5) = [ & - option ('help',.false.,'h','Display this help message',''), & - option ('mycall',.true.,'c','mycall',''), & + type (option) :: long_options(6) = [ & + option ('dxcall',.true.,'d','hiscall',''), & 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') - display_help = .true. - case ('c') - read (optarg(:narglen), *) mycall + case ('d') + read (optarg(:narglen), *) hiscall case ('e') t0=1e-4 + case ('h') + display_help = .true. + case ('m') + read (optarg(:narglen), *) mycall case ('n') read (optarg(:narglen), *) ntol - case ('x') - read (optarg(:narglen), *) hiscall + case ('s') + 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