diff --git a/CMakeLists.txt b/CMakeLists.txt index 640777d73..8671a03a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -534,6 +534,7 @@ set (wsjt_FSRCS lib/77bit/my_hash.f90 lib/wsprd/osdwspr.f90 lib/ft8/osd174_91.f90 + lib/osd128_90.f90 lib/pctile.f90 lib/peakdt9.f90 lib/peakup.f90 diff --git a/lib/ldpcsim128_90.f90 b/lib/ldpcsim128_90.f90 index 2b6baf407..28e7fef03 100644 --- a/lib/ldpcsim128_90.f90 +++ b/lib/ldpcsim128_90.f90 @@ -4,8 +4,8 @@ program ldpcsim128_90 ! the second incarnation of MSK144. use packjt77 - integer, parameter:: NRECENT=10, N=128, K=90, M=N-K - character*12 recent_calls(NRECENT) + integer, parameter:: N=128, K=90, M=N-K +! character*12 recent_calls(NRECENT) character*37 msg,msgsent,msgreceived character*77 c77 character*8 arg @@ -16,27 +16,29 @@ program ldpcsim128_90 integer nerrtot(0:N),nerrdec(0:N) logical unpk77_success real*8 rxdata(N), rxavgd(N) - real llr(N) + real llr(N),llra(N) - do i=1,NRECENT + do i=1,MAXNRECENT recent_calls(i)=' ' enddo nerrtot=0 nerrdec=0 nargs=iargc() - if(nargs.ne.4) then - print*,'Usage: ldpcsim niter navg #trials s ' - print*,'eg: ldpcsim 10 1 1000 0.75' + if(nargs.ne.5) then + print*,'Usage: ldpcsim niter ndeep navg #trials s ' + print*,'eg: ldpcsim 10 2 1 1000 0.75' return endif call getarg(1,arg) read(arg,*) max_iterations call getarg(2,arg) - read(arg,*) navg + read(arg,*) ndeep call getarg(3,arg) - read(arg,*) ntrials + read(arg,*) navg call getarg(4,arg) + read(arg,*) ntrials + call getarg(5,arg) read(arg,*) s rate=real(77)/real(N) @@ -62,7 +64,7 @@ program ldpcsim128_90 ! call init_random_seed() write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr" - do idb = 14,0,-1 + do idb = 6,6,-1 db=idb/2.0-1.0 sigma=1/sqrt( 2*rate*(10**(db/10.0)) ) ngood=0 @@ -103,8 +105,13 @@ program ldpcsim128_90 apmask=0 ! max_iterations is max number of belief propagation iterations + call bpdecode128_90(llr, apmask, max_iterations, message77, cw, nharderrors, niterations) + if(ndeep.ge.0 .and. nharderrors.lt.0) then + call osd128_90(llr, apmask, ndeep, message77, cw, nharderrors, dmin) + endif + ! If the decoder finds a valid codeword, nharderrors will be .ge. 0. if( nharderrors .ge. 0 ) then write(c77,'(77i1)') message77 diff --git a/lib/msk144decodeframe.f90 b/lib/msk144decodeframe.f90 index a2dbd0dff..aeb9bd891 100644 --- a/lib/msk144decodeframe.f90 +++ b/lib/msk144decodeframe.f90 @@ -1,4 +1,4 @@ -subroutine msk144decodeframe(c,softbits,msgreceived,nsuccess) +subroutine msk144decodeframe(c,softbits,msgreceived,nsuccess,doosd) ! use timer_module, only: timer use packjt77 parameter (NSPM=864) @@ -14,7 +14,7 @@ subroutine msk144decodeframe(c,softbits,msgreceived,nsuccess) real pp(12) real softbits(144) real llr(128) - logical first,unpk77_success + logical first,unpk77_success,doosd data first/.true./ data s8/0,1,1,1,0,0,1,0/ save first,cb,fs,pi,twopi,dt,s8,pp @@ -96,6 +96,7 @@ subroutine msk144decodeframe(c,softbits,msgreceived,nsuccess) ! call timer('bpdec128_90 ',0) apmask=0 call bpdecode128_90(llr,apmask,max_iterations,decoded77,cw,nharderror,niterations) + if(doosd .and. nharderror .lt. 0) call osd128_90(llr,apmask,3,decoded77,cw,nharderror,dmin) ! call timer('bpdec128_90 ',1) if( nharderror .ge. 0 .and. nharderror .lt. 18 ) then nsuccess=1 diff --git a/lib/msk144spd.f90 b/lib/msk144spd.f90 index 2bc6dcfef..85758a1dd 100644 --- a/lib/msk144spd.f90 +++ b/lib/msk144spd.f90 @@ -31,6 +31,7 @@ subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg,ct, & real tpat(NPATTERNS) real*8 dt, df, fs, pi, twopi logical first + logical doosd data first/.true./ data navpatterns/ & 0,1,0, & @@ -43,6 +44,7 @@ subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg,ct, & save df,first,fs,pi,twopi,dt,tframe,rcw + if( ntol .le. 50 ) doosd = .true. if(first) then nmatchedfilter=1 ! define half-sine pulse and raised-cosine edge window @@ -179,7 +181,7 @@ subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg,ct, & if( is.eq.2) ic0=max(1,ic0-1) if( is.eq.3) ic0=min(NSPM,ic0+1) ct=cshift(c,ic0-1) - call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess) + call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess,doosd) if( ndecodesuccess .gt. 0 ) then tret=(nstart(icand)+NSPM/2)/fs fret=fest diff --git a/lib/mskrtd.f90 b/lib/mskrtd.f90 index 7cc688071..54bcf29e7 100644 --- a/lib/mskrtd.f90 +++ b/lib/mskrtd.f90 @@ -45,6 +45,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, & logical*1 seenb4 logical*1 bflag logical*1 bvar + logical*1 doosd data first/.true./ data iavpatterns/ & @@ -143,6 +144,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, & npat=NPATTERNS if( ndepth .eq. 1 ) npat=0 if( ndepth .eq. 2 ) npat=2 + if( ntol .le. 50 ) doosd = .true. do iavg=1,npat iavmask=iavpatterns(1:8,iavg) navg=sum(iavmask) @@ -158,7 +160,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, & if(is.eq.2) ic0=max(1,ic0-1) if(is.eq.3) ic0=min(NSPM,ic0+1) ct=cshift(c,ic0-1) - call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess) + call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess,doosd) if(ndecodesuccess .gt. 0) then tdec=tsec+xmc(iavg)*tframe goto 900 diff --git a/lib/osd128_90.f90 b/lib/osd128_90.f90 new file mode 100644 index 000000000..d26aa4fea --- /dev/null +++ b/lib/osd128_90.f90 @@ -0,0 +1,371 @@ +subroutine osd128_90(llr,apmask,ndeep,message77,cw,nhardmin,dmin) +! +! An ordered-statistics decoder for the (128,90) code. +! +integer, parameter:: N=128, K=90, M=N-K +integer*1 apmask(N),apmaskr(N) +integer*1 gen(K,N) +integer*1 genmrb(K,N),g2(N,K) +integer*1 temp(K),m0(K),me(K),mi(K),misub(K),e2sub(N-K),e2(N-K),ui(N-K) +integer*1 r2pat(N-K) +integer indices(N),nxor(N) +integer*1 cw(N),ce(N),c0(N),hdec(N) +integer*1 decoded(K) +integer*1 message77(77) +integer indx(N) +real llr(N),rx(N),absrx(N) + +include "ldpc_128_90_generator.f90" + +logical first,reset +data first/.true./ +save first,gen + +if( first ) then ! fill the generator matrix + gen=0 + do i=1,M + do j=1,23 + read(g(i)(j:j),"(Z1)") istr + ibmax=4 + if(j.eq.23) ibmax=2 + do jj=1, ibmax + irow=(j-1)*4+jj + if( btest(istr,4-jj) ) gen(irow,K+i)=1 + enddo + enddo + enddo + do irow=1,K + gen(irow,irow)=1 + enddo +first=.false. +endif + +rx=llr +apmaskr=apmask + +! Hard decisions on the received word. +hdec=0 +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 decreasing reliability. +do i=1,N + 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 in positions 1:K in order of decreasing reliability (more or less). +do id=1,K ! diagonal element indices + do icol=id,K+20 ! The 20 is ad hoc - beware + iflag=0 + if( genmrb(id,icol) .eq. 1 ) then + iflag=1 + if( icol .ne. id ) then ! reorder column + temp(1:K)=genmrb(1:K,id) + genmrb(1:K,id)=genmrb(1:K,icol) + genmrb(1:K,icol)=temp(1:K) + itmp=indices(id) + indices(id)=indices(icol) + indices(icol)=itmp + endif + do ii=1,K + if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then + genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N)) + endif + enddo + exit + endif + enddo +enddo + +g2=transpose(genmrb) + +! The hard decisions for the K MRB bits define the order 0 message, m0. +! Encode m0 using the modified generator matrix to find the "order 0" codeword. +! Flip various combinations of bits in m0 and re-encode to generate a list of +! codewords. Return the member of the list that has the smallest Euclidean +! distance to the received word. + +hdec=hdec(indices) ! hard decisions from received symbols +m0=hdec(1:K) ! zero'th order message +absrx=absrx(indices) +rx=rx(indices) +apmaskr=apmaskr(indices) + +call mrbencode90(m0,c0,g2,N,K) +nxor=ieor(c0,hdec) +nhardmin=sum(nxor) +dmin=sum(nxor*absrx) + +cw=c0 +ntotal=0 +nrejected=0 + +if(ndeep.eq.0) goto 998 ! norder=0 +if(ndeep.gt.5) ndeep=5 +if( ndeep.eq. 1) then + nord=1 + npre1=0 + npre2=0 + nt=12 + ntheta=3 +elseif(ndeep.eq.2) then + nord=1 + npre1=1 + npre2=0 + nt=18 + ntheta=4 +elseif(ndeep.eq.3) then + nord=1 + npre1=1 + npre2=0 + nt=12 + ntheta=4 +elseif(ndeep.eq.4) then + nord=1 + npre1=1 + npre2=1 + nt=12 + ntheta=4 + ntau=15 +elseif(ndeep.eq.5) then + nord=1 + npre1=1 + npre2=1 + nt=12 + ntheta=4 + ntau=5 +endif + +do iorder=1,nord + misub(1:K-iorder)=0 + misub(K-iorder+1:K)=1 + iflag=K-iorder+1 + do while(iflag .ge.0) + if(iorder.eq.nord .and. npre1.eq.0) then + iend=iflag + else + iend=1 + endif + do n1=iflag,iend,-1 + mi=misub + mi(n1)=1 + if(any(iand(apmaskr(1:K),mi).eq.1)) cycle + ntotal=ntotal+1 + me=ieor(m0,mi) + if(n1.eq.iflag) then + call mrbencode90(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + e2=e2sub + nd1Kpt=sum(e2sub(1:nt))+1 + d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K)) + else + e2=ieor(e2sub,g2(K+1:N,n1)) + nd1Kpt=sum(e2(1:nt))+2 + endif + if(nd1Kpt .le. ntheta) then + call mrbencode90(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + if(n1.eq.iflag) then + dd=d1+sum(e2sub*absrx(K+1:N)) + else + dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N)) + endif + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + nd1Kptbest=nd1Kpt + endif + else + nrejected=nrejected+1 + endif + enddo +! Get the next test error pattern, iflag will go negative +! when the last pattern with weight iorder has been generated. + call nextpat90(misub,k,iorder,iflag) + enddo +enddo + +if(npre2.eq.1) then + reset=.true. + ntotal=0 + do i1=K,1,-1 + do i2=i1-1,1,-1 + ntotal=ntotal+1 + mi(1:ntau)=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2)) + call boxit90(reset,mi(1:ntau),ntau,ntotal,i1,i2) + enddo + enddo + + ncount2=0 + ntotal2=0 + reset=.true. +! Now run through again and do the second pre-processing rule + misub(1:K-nord)=0 + misub(K-nord+1:K)=1 + iflag=K-nord+1 + do while(iflag .ge.0) + me=ieor(m0,misub) + call mrbencode90(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + do i2=0,ntau + ntotal2=ntotal2+1 + ui=0 + if(i2.gt.0) ui(i2)=1 + r2pat=ieor(e2sub,ui) +778 continue + call fetchit90(reset,r2pat(1:ntau),ntau,in1,in2) + if(in1.gt.0.and.in2.gt.0) then + ncount2=ncount2+1 + mi=misub + mi(in1)=1 + mi(in2)=1 + if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:K),mi).eq.1)) cycle + me=ieor(m0,mi) + call mrbencode90(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + dd=sum(nxor*absrx) + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + endif + goto 778 + endif + enddo + call nextpat90(misub,K,nord,iflag) + enddo +endif + +998 continue +! Re-order the codeword to [message bits][parity bits] format. +cw(indices)=cw +hdec(indices)=hdec +decoded=cw(1:K) +call chkcrc13a(decoded,nbadcrc) +message77=decoded(1:77) +if(nbadcrc.eq.1) nhardmin=-nhardmin + +return +end subroutine osd128_90 + +subroutine mrbencode90(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 mrbencode90 + +subroutine nextpat90(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 + do i=1,k ! iflag will point to the lowest-index 1 in mi + if(mi(i).eq.1) then + iflag=i + exit + endif + enddo + return +end subroutine nextpat90 + +subroutine boxit90(reset,e2,ntau,npindex,i1,i2) + integer*1 e2(1:ntau) + integer indexes(5000,2),fp(0:525000),np(5000) + logical reset + common/boxes/indexes,fp,np + + if(reset) then + patterns=-1 + fp=-1 + np=-1 + sc=-1 + indexes=-1 + reset=.false. + endif + + indexes(npindex,1)=i1 + indexes(npindex,2)=i2 + ipat=0 + do i=1,ntau + if(e2(i).eq.1) then + ipat=ipat+ishft(1,ntau-i) + endif + enddo + + ip=fp(ipat) ! see what's currently stored in fp(ipat) + if(ip.eq.-1) then + fp(ipat)=npindex + else + do while (np(ip).ne.-1) + ip=np(ip) + enddo + np(ip)=npindex + endif + return +end subroutine boxit90 + +subroutine fetchit90(reset,e2,ntau,i1,i2) + integer indexes(5000,2),fp(0:525000),np(5000) + integer lastpat + integer*1 e2(ntau) + logical reset + common/boxes/indexes,fp,np + save lastpat,inext + + if(reset) then + lastpat=-1 + reset=.false. + endif + + ipat=0 + do i=1,ntau + if(e2(i).eq.1) then + ipat=ipat+ishft(1,ntau-i) + endif + enddo + index=fp(ipat) + + if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices + i1=indexes(index,1) + i2=indexes(index,2) + inext=np(index) + elseif(lastpat.eq.ipat .and. inext.gt.0) then + i1=indexes(inext,1) + i2=indexes(inext,2) + inext=np(inext) + else + i1=-1 + i2=-1 + inext=-1 + endif + lastpat=ipat + return +end subroutine fetchit90 +