diff --git a/lib/fsk4hf/ft8b.f90 b/lib/fsk4hf/ft8b.f90 index 630d8cd89..3d3a62299 100644 --- a/lib/fsk4hf/ft8b.f90 +++ b/lib/fsk4hf/ft8b.f90 @@ -12,7 +12,7 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & real a(5) real s1(0:7,ND),s2(0:7,NN) real ps(0:7) - real rxdata(3*ND),rxdatap(3*ND) + real bmeta(3*ND),bmetb(3*ND),bmetap(3*ND) real llr(3*ND),llra(3*ND),llr0(3*ND),llrap(3*ND) !Soft symbols real dd0(15*12000) integer*1 decoded(KK),apmask(3*ND),cw(3*ND) @@ -122,7 +122,7 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & csymb=cmplx(0.0,0.0) if( i1.ge.1 .and. i1+31 .le. NP2 ) csymb=cd0(i1:i1+31) call four2a(csymb,32,1,-1,1) - s2(0:7,k)=abs(csymb(1:8)) + s2(0:7,k)=abs(csymb(1:8))/1e3 enddo ! sync quality check @@ -154,31 +154,30 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & enddo do j=1,ND - ps=s1(0:7,j) - 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)) i4=3*j-2 i2=3*j-1 i1=3*j - rxdata(i4)=r4 - rxdata(i2)=r2 - rxdata(i1)=r1 - rxdatap(i4)=r4 - rxdatap(i2)=r2 - rxdatap(i1)=r1 + ps=s1(0:7,j) + 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)) + bmeta(i4)=r4 + bmeta(i2)=r2 + bmeta(i1)=r1 + bmetap(i4)=r4 + bmetap(i2)=r2 + bmetap(i1)=r1 if(nQSOProgress .eq. 0 .or. nQSOProgress .eq. 5) then ! When bits 88:115 are set as ap bits, bit 115 lives in symbol 39 along ! with no-ap bits 116 and 117. Take care of metrics for bits 116 and 117. if(j.eq.39) then ! take care of bits that live in symbol 39 if(apsym(28).lt.0) then - rxdatap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1)) - rxdatap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2)) + bmetap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1)) + bmetap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2)) else - rxdatap(i2)=max(ps(6),ps(7))-max(ps(4),ps(5)) - rxdatap(i1)=max(ps(5),ps(7))-max(ps(4),ps(6)) + bmetap(i2)=max(ps(6),ps(7))-max(ps(4),ps(5)) + bmetap(i1)=max(ps(5),ps(7))-max(ps(4),ps(6)) endif endif endif @@ -187,43 +186,34 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & ! with ap bits 116 and 117. Take care of metric for bit 115. ! if(j.eq.39) then ! take care of bit 115 ! iii=2*(apsym(29)+1)/2 + (apsym(30)+1)/2 ! known values of bits 116 & 117 -! if(iii.eq.0) rxdatap(i4)=ps(4)-ps(0) -! if(iii.eq.1) rxdatap(i4)=ps(5)-ps(1) -! if(iii.eq.2) rxdatap(i4)=ps(6)-ps(2) -! if(iii.eq.3) rxdatap(i4)=ps(7)-ps(3) +! if(iii.eq.0) bmetap(i4)=ps(4)-ps(0) +! if(iii.eq.1) bmetap(i4)=ps(5)-ps(1) +! if(iii.eq.2) bmetap(i4)=ps(6)-ps(2) +! if(iii.eq.3) bmetap(i4)=ps(7)-ps(3) ! endif ! bit 144 lives in symbol 48 and will be 1 if it is set as an ap bit. ! take care of metrics for bits 142 and 143 if(j.eq.48) then ! bit 144 is always 1 - rxdatap(i4)=max(ps(5),ps(7))-max(ps(1),ps(3)) - rxdatap(i2)=max(ps(3),ps(7))-max(ps(1),ps(5)) + bmetap(i4)=max(ps(5),ps(7))-max(ps(1),ps(3)) + bmetap(i2)=max(ps(3),ps(7))-max(ps(1),ps(5)) endif ! bit 154 lives in symbol 52 and will be 0 if it is set as an ap bit ! take care of metrics for bits 155 and 156 if(j.eq.52) then ! bit 154 will be 0 if it is set as an ap bit. - rxdatap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1)) - rxdatap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2)) + bmetap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1)) + bmetap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2)) endif enddo - rxav=sum(rxdata)/(3.0*ND) - 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 -! Let's just assume that rxsig is OK for rxdatap too... - rxdatap=rxdatap/rxsig + call normalizebmet(bmeta,3*ND) + call normalizebmet(bmetap,3*ND) ss=0.84 - llr0=2.0*rxdata/(ss*ss) - llra=2.0*rxdatap/(ss*ss) ! llr's for use with ap + llr0=2.0*bmeta/(ss*ss) + llra=2.0*bmetap/(ss*ss) ! llr's for use with ap apmag=4.0 ! pass # @@ -245,11 +235,8 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & do ipass=1,npasses llr=llr0 - if(ipass.ne.2 .and. ipass.ne.3) nblank=0 - if(ipass.eq.2) nblank=24 - if(ipass.eq.3) nblank=48 - if(nblank.gt.0) llr(1:nblank)=0. - + if(ipass.eq.2) llr(1:24)=0. + if(ipass.eq.3) llr(1:48)=0. if(ipass.le.3) then apmask=0 llrap=llr @@ -308,16 +295,17 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & call timer('bpd174 ',1) dmin=0.0 if(ndepth.eq.3 .and. nharderrors.lt.0) then - norder=2 + ndeep=3 if(abs(nfqso-f1).le.napwid .or. abs(nftx-f1).le.napwid) then if((ipass.eq.2 .or. ipass.eq.3) .and. .not.nagain) then - norder=2 + ndeep=3 else - norder=3 ! for nagain, use norder=3 for all passes + ndeep=4 endif endif + if(nagain) ndeep=5 call timer('osd174 ',0) - call osd174(llrap,apmask,norder,decoded,cw,nharderrors,dmin) + call osd174(llrap,apmask,ndeep,decoded,cw,nharderrors,dmin) call timer('osd174 ',1) endif nbadcrc=1 @@ -356,3 +344,18 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & return end subroutine ft8b + +subroutine normalizebmet(bmet,n) + real bmet(n) + + bmetav=sum(bmet)/real(n) + bmet2av=sum(bmet*bmet)/real(n) + var=bmet2av-bmetav*bmetav + if( var .gt. 0.0 ) then + bmetsig=sqrt(var) + else + bmetsig=sqrt(bmet2av) + endif + bmet=bmet/bmetsig + return +end subroutine normalizebmet diff --git a/lib/fsk4hf/ldpcsim174.f90 b/lib/fsk4hf/ldpcsim174.f90 index 1b3b23b53..9c4a72893 100644 --- a/lib/fsk4hf/ldpcsim174.f90 +++ b/lib/fsk4hf/ldpcsim174.f90 @@ -39,16 +39,16 @@ nmpcbad=0 ! Used to collect the number of errors in the message+crc part of the nargs=iargc() if(nargs.ne.4) then - print*,'Usage: ldpcsim niter norder #trials s ' + print*,'Usage: ldpcsim niter ndepth #trials s ' print*,'eg: ldpcsim 10 2 1000 0.84' - print*,'belief propagation iterations: niter, ordered-statistics order: norder' + print*,'belief propagation iterations: niter, ordered-statistics depth: ndepth' print*,'If s is negative, then value is ignored and sigma is calculated from SNR.' return endif call getarg(1,arg) read(arg,*) max_iterations call getarg(2,arg) -read(arg,*) norder +read(arg,*) ndepth call getarg(3,arg) read(arg,*) ntrials call getarg(4,arg) @@ -128,13 +128,14 @@ allocate ( rxdata(N), llr(N) ) call encode174(msgbits,codeword) call init_random_seed() - call sgran() +! call sgran() write(*,*) 'codeword' write(*,'(22(8i1,1x))') codeword write(*,*) "Es/N0 SNR2500 ngood nundetected nbadcrc sigma" do idb = 20,-10,-1 +!do idb = -3,-3,-1 db=idb/2.0-1.0 sigma=1/sqrt( 2*(10**(db/10.0)) ) ngood=0 @@ -189,7 +190,7 @@ do idb = 20,-10,-1 ! max_iterations is max number of belief propagation iterations call bpdecode174(llr, apmask, max_iterations, decoded, cw, nharderrors,niterations) - if( norder .ge. 0 .and. nharderrors .lt. 0 ) call osd174(llr, apmask, norder, decoded, cw, nharderrors, dmin) + if( ndepth .ge. 0 .and. nharderrors .lt. 0 ) call osd174(llr, apmask, ndepth, decoded, cw, nharderrors, dmin) ! If the decoder finds a valid codeword, nharderrors will be .ge. 0. if( nharderrors .ge. 0 ) then call extractmessage174(decoded,msgreceived,ncrcflag,recent_calls,nrecent) @@ -206,7 +207,6 @@ do idb = 20,-10,-1 endif enddo nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1 - if( ncrcflag .eq. 1 ) then if( nueflag .eq. 0 ) then ngood=ngood+1 diff --git a/lib/fsk4hf/osd174.f90 b/lib/fsk4hf/osd174.f90 index fc5ac6441..1ec1c4e5b 100644 --- a/lib/fsk4hf/osd174.f90 +++ b/lib/fsk4hf/osd174.f90 @@ -1,4 +1,4 @@ -subroutine osd174(llr,apmask,norder,decoded,cw,nhardmin,dmin) +subroutine osd174(llr,apmask,ndeep,decoded,cw,nhardmin,dmin) ! ! An ordered-statistics decoder for the (174,87) code. ! @@ -7,13 +7,14 @@ include "ldpc_174_87_params.f90" 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) +integer*1 temp(K),m0(K),me(K),mep(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 cw(N),ce(N),c0(N),hdec(N),qi(N) integer*1 decoded(K) integer indx(N) real llr(N),rx(N),absrx(N) -logical first +logical first,reset data first/.true./ save first,gen @@ -38,7 +39,6 @@ endif rx=llr(colorder+1) apmaskr=apmask(colorder+1) - ! Hard decisions on the received word. hdec=0 where(rx .ge. 0) hdec=1 @@ -83,9 +83,8 @@ 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. A pre-processing step selects a subset of these codewords. -! Return the member of the subset with the smallest Euclidean distance to the -! the received word. +! 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 @@ -101,26 +100,42 @@ dmin=sum(nxor*absrx) cw=c0 ntotal=0 nrejected=0 -nt=40 ! Count the errors in the nt best bits in the redundancy part of the cw -ntheta=12 ! Reject the codeword without computing distance if # errors exceeds ntheta -! norder should be 1, 2, or 3. -! if norder = 1, do one loop, no pre-processing -! if norder = 2, do norder=1, then norder=2 using first W&H pre-processing rule -! if norder = 3, do norder=2, then norder=3 using first W&H pre-processing rule - -if(norder.lt.1) goto 998 ! norder=0 -if(norder.gt.3) norder=3 - -if( norder.eq. 1) then +if(ndeep.eq.0) goto 998 ! norder=0 +if(ndeep.gt.5) ndeep=5 +if( ndeep.eq. 1) then nord=1 - npre=0 -elseif(norder.eq.2) then + npre1=0 + npre2=0 + nt=40 + ntheta=12 +elseif(ndeep.eq.2) then nord=1 - npre=1 -elseif(norder.eq.3) then + npre1=1 + npre2=0 + nt=40 + ntheta=12 +elseif(ndeep.eq.3) then + nord=1 + npre1=1 + npre2=1 + nt=40 + ntheta=12 + ntau=14 +elseif(ndeep.eq.4) then nord=2 - npre=1 + npre1=1 + npre2=0 + nt=40 + ntheta=12 + ntau=19 +elseif(ndeep.eq.5) then + nord=2 + npre1=1 + npre2=1 + nt=40 + ntheta=12 + ntau=19 endif do iorder=1,nord @@ -134,7 +149,7 @@ do iorder=1,nord iflag=K-1 endif do while(iflag .ge.0) - if(iorder.eq.nord .and. npre.eq.0) then + if(iorder.eq.nord .and. npre1.eq.0) then iend=iflag else iend=1 @@ -179,7 +194,62 @@ do iorder=1,nord enddo enddo -!write(*,*) 'rejected, total, nd1Kptbest: ',nrejected, ntotal, nd1Kptbest +if(npre2.eq.1) then + reset=.true. + ntotal=0 + do i1=K,1,-1 + do i2=i1-1,1,-1 + ntotal=ntotal+1 + mi=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2)) + call boxit(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 + if(nord.eq.1) then + misub(1:K-1)=0 + misub(K)=1 + iflag=K + elseif(nord.eq.2) then + misub(1:K-1)=0 + misub(K-1:K)=1 + iflag=K-1 + endif + do while(iflag .ge.0) + me=ieor(m0,misub) + call mrbencode(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 fetchit(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 mrbencode(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 nextpat(misub,K,nord,iflag) + enddo +endif 998 continue ! Re-order the codeword to place message bits at the end. @@ -230,3 +300,78 @@ subroutine nextpat(mi,k,iorder,iflag) enddo return end subroutine nextpat + +subroutine boxit(reset,e2,ntau,npindex,i1,i2) + integer*1 e2(1:ntau) + integer indexes(4000,2),fp(0:525000),np(4000) + 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 boxit + +subroutine fetchit(reset,e2,ntau,i1,i2) + integer indexes(4000,2),fp(0:525000),np(4000) + integer lastpat + integer*1 e2(ntau) + logical reset + common/boxes/indexes,fp,np + save lastpat,inext,first + + 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 fetchit + diff --git a/lib/ft8_decode.f90 b/lib/ft8_decode.f90 index caf235a3d..20cbaf863 100644 --- a/lib/ft8_decode.f90 +++ b/lib/ft8_decode.f90 @@ -67,19 +67,24 @@ contains ! For now: ! ndepth=1: no subtraction, 1 pass, belief propagation only ! ndepth=2: subtraction, 2 passes, belief propagation only -! ndepth=3: subtraction, 2 passes, bp+osd2 at and near nfqso +! ndepth=3: subtraction, 2 passes, bp+osd if(ndepth.eq.1) npass=1 - if(ndepth.ge.2) npass=2 + if(ndepth.ge.2) npass=3 do ipass=1,npass newdat=.true. ! Is this a problem? I hijacked newdat. + syncmin=1.5 if(ipass.eq.1) then lsubtract=.true. if(ndepth.eq.1) lsubtract=.false. - syncmin=1.5 - else - lsubtract=.false. - syncmin=1.5 + elseif(ipass.eq.2) then + n2=ndecodes + if(ndecodes.eq.0) cycle + lsubtract=.true. + elseif(ipass.eq.3) then + if((ndecodes-n2).eq.0) cycle + lsubtract=.false. endif + call timer('sync8 ',0) call sync8(dd,ifa,ifb,syncmin,nfqso,s,candidate,ncand) call timer('sync8 ',1)