diff --git a/CMakeLists.txt b/CMakeLists.txt index 131663a75..7cd7ae83b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -428,6 +428,7 @@ set (wsjt_FSRCS lib/encode_128_90.f90 lib/ft8/encode174.f90 lib/ft8/encode174_91.f90 + lib/ft8/encode174_91_nocrc.f90 lib/entail.f90 lib/ephem.f90 lib/extract.f90 @@ -603,12 +604,14 @@ set (wsjt_FSRCS lib/fsk4hf/genwsprcpm.f90 lib/fsk4hf/encode204.f90 lib/fsk4hf/decode174_74.f90 + lib/ft8/decode174_91.f90 lib/fsk4hf/decode240_101.f90 lib/fsk4hf/decode280_101.f90 lib/fsk4hf/ldpcsim174_91.f90 lib/fsk4hf/ldpcsim174_74.f90 lib/fsk4hf/ldpcsim240_101.f90 lib/fsk4hf/ldpcsim280_101.f90 + lib/fsk4hf/get_crc14.f90 lib/fsk4hf/get_crc24.f90 lib/fsk4hf/encode174_74.f90 lib/fsk4hf/encode240_101.f90 diff --git a/lib/fsk4hf/get_crc14.f90 b/lib/fsk4hf/get_crc14.f90 index e1f8df8ff..635808972 100644 --- a/lib/fsk4hf/get_crc14.f90 +++ b/lib/fsk4hf/get_crc14.f90 @@ -1,15 +1,19 @@ -subroutine get_crc14(mc,ncrc) - +subroutine get_crc14(mc,len,ncrc) +! +! 1. To calculate 14-bit CRC, mc(1:len-14) is the message and mc(len-13:len) are zero. +! 2. To check a received CRC, mc(1:len is the received message plus CRC. +! ncrc will be zero if the received message/CRC are consistent +! character c14*14 - - integer*1 mc(68),r(15),p(15) + integer*1 mc(len) + integer*1 r(15),p(15) integer ncrc ! polynomial for 14-bit CRC 0x6757 data p/1,1,0,0,1,1,1,0,1,0,1,0,1,1,1/ ! divide by polynomial r=mc(1:15) - do i=0,53 + do i=0,len-15 r(15)=mc(i+15) r=mod(r+r(1)*p,2) r=cshift(r,1) @@ -17,6 +21,5 @@ subroutine get_crc14(mc,ncrc) write(c14,'(14b1)') r(1:14) read(c14,'(b14.14)') ncrc -! mc(55:68)=r(1:14) end subroutine get_crc14 diff --git a/lib/fsk4hf/ldpcsim174_91.f90 b/lib/fsk4hf/ldpcsim174_91.f90 index a753dd8ad..2531910e0 100644 --- a/lib/fsk4hf/ldpcsim174_91.f90 +++ b/lib/fsk4hf/ldpcsim174_91.f90 @@ -10,7 +10,7 @@ character*6 grid character*96 tmpchar integer*1, allocatable :: codeword(:), decoded(:), message(:) integer*1 msgbits(77) -integer*1 message77(77) +integer*1 message77(77),message91(91) integer*1 apmask(N), cw(N) integer nerrtot(0:N),nerrdec(0:N) logical unpk77_success @@ -21,11 +21,12 @@ nerrtot=0 nerrdec=0 nargs=iargc() -if(nargs.ne.4) then - print*,'Usage: ldpcsim niter ndepth #trials s ' - print*,'eg: ldpcsim 10 2 1000 0.84' +if(nargs.ne.6) then + print*,'Usage: ldpcsim niter ndepth #trials s Keff BPOSD' + print*,'eg: ldpcsim 10 2 1000 0.84 91 1' print*,'belief propagation iterations: niter, ordered-statistics depth: ndepth' print*,'If s is negative, then value is ignored and sigma is calculated from SNR.' + print*,'If BPOSD=0, no coupling. BPOSD=1, BP output to OSD input.' return endif call getarg(1,arg) @@ -36,9 +37,13 @@ call getarg(3,arg) read(arg,*) ntrials call getarg(4,arg) read(arg,*) s +call getarg(5,arg) +read(arg,*) Keff +call getarg(6,arg) +read(arg,*) nbposd ! scale Eb/No for a (174,91) code -rate=real(K)/real(N) +rate=real(Keff)/real(N) write(*,*) "rate: ",rate write(*,*) "niter= ",max_iterations," s= ",s @@ -46,8 +51,7 @@ write(*,*) "niter= ",max_iterations," s= ",s allocate ( codeword(N), decoded(K), message(K) ) allocate ( rxdata(N), llr(N) ) -! msg="K1JT K9AN EN50" - msg="G4WJS K9AN EN50" + msg="K9ABC K1ABC FN20" i3=0 n3=1 call pack77(msg,i3,n3,c77) !Pack into 12 6-bit bytes @@ -55,20 +59,19 @@ allocate ( rxdata(N), llr(N) ) write(*,*) "message sent ",msgsent read(c77,'(77i1)') msgbits(1:77) - write(*,*) 'message' write(*,'(a71,1x,a3,1x,a3)') c77(1:71),c77(72:74),c77(75:77) - call encode174_91(msgbits,codeword) - call init_random_seed() + call encode174_91(msgbits,codeword) + write(*,*) 'crc14' + write(*,'(14i1)') codeword(78:91) write(*,*) 'codeword' write(*,'(22(8i1,1x))') codeword -write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr" -do idb = 20,-10,-1 -!do idb = 0,0,-1 +write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr" +do idb = 10,-4,-1 db=idb/2.0-1.0 sigma=1/sqrt( 2*rate*(10**(db/10.0)) ) ngood=0 @@ -101,21 +104,26 @@ do idb = 20,-10,-1 llr(1:nap)=5*(2.0*msgbits(1:nap)-1.0) apmask=0 apmask(1:nap)=1 - ! max_iterations is max number of belief propagation iterations - call bpdecode174_91(llr, apmask, max_iterations, message77, cw, nharderrors,niterations,ncheck) - if( ndepth .ge. 0 .and. nharderrors .lt. 0 ) call osd174_91(llr, apmask, ndepth, message77, cw, nharderrors, dmin) -!call decode174_91(llr, apmask, max_iterations, message77, cw, nharderrors,niterations,ncheck,dmin,nsuper) + call bpdecode174_91(llr, apmask, max_iterations, message77, cw, nhardbp,niterations,ncheck) + if( ndepth .ge. 0 .and. nhardbp .lt. 0 ) then + dmin=0.0 + if(nbposd.eq.0) then + call osd174_91(llr,Keff,apmask,ndepth,message91,cw,nhardosd,dmin) + else + maxsuper=2 + call decode174_91(llr,Keff,ndepth,apmask,maxsuper,message91,cw,nhardosd,niterations,ncheck,dmin) + endif ! If the decoder finds a valid codeword, nharderrors will be .ge. 0. - if( nharderrors .ge. 0 ) then - call extractmessage77(message77,msgreceived) - nhw=count(cw.ne.codeword) - if(nhw.eq.0) then ! this is a good decode - ngood=ngood+1 - nerrdec(nerr)=nerrdec(nerr)+1 - else - nue=nue+1 - endif + endif + if( nhardbp .ge. 0 .or. nhardosd.ge.0 ) then + nhw=count(cw.ne.codeword) + if(nhw.eq.0) then ! this is a good decode + ngood=ngood+1 + nerrdec(nerr)=nerrdec(nerr)+1 + else + nue=nue+1 + endif endif nsumerr=nsumerr+nerr enddo diff --git a/lib/ft4_decode.f90 b/lib/ft4_decode.f90 index 3de9f89f5..6207653d4 100644 --- a/lib/ft4_decode.f90 +++ b/lib/ft4_decode.f90 @@ -56,6 +56,7 @@ contains integer apmy_ru(28),aphis_fd(28) integer*2 iwave(NMAX) !Raw received data integer*1 message77(77),rvec(77),apmask(2*ND),cw(2*ND) + integer*1 message91(91) integer*1 hbits(2*NN) integer i4tone(103) integer nappasses(0:5) ! # of decoding passes for QSO States 0-5 @@ -412,12 +413,19 @@ contains call timer('bpdec174',1) if(doosd .and. nharderror.lt.0) then - ndeep=3 +! ndeep=3 + ndeep=2 if(abs(nfqso-f1).le.napwid) then - ndeep=4 +! ndeep=4 + ndeep=3 endif call timer('osd174_91 ',0) - call osd174_91(llr,apmask,ndeep,message77,cw,nharderror,dmin) + Keff=91 + maxsuper=2 +! call osd174_91(llr,Keff,apmask,ndeep,message91,cw,nharderror,dmin) + call decode174_91(llr,Keff,ndeep,apmask,maxsuper,message91,cw,nharderror, & + niterations,ncheck,dmin) + message77=message91(1:77) call timer('osd174_91 ',1) endif diff --git a/lib/ft8/decode174_91.f90 b/lib/ft8/decode174_91.f90 index 9538141d2..baa9b8946 100644 --- a/lib/ft8/decode174_91.f90 +++ b/lib/ft8/decode174_91.f90 @@ -1,13 +1,12 @@ -subroutine decode174_91(llr,apmask,maxiterations,message77,cw,nharderror,iter,ncheck,dmin,isuper) +subroutine decode174_91(llr,Keff,ndeep,apmask,maxsuper,message91,cw,nharderror,iter,ncheck,dmin) ! ! A hybrid bp/osd decoder for the (174,91) code. ! -use iso_c_binding, only: c_loc,c_size_t -use crc integer, parameter:: N=174, K=91, M=N-K integer*1 cw(N),apmask(N) integer*1 decoded(K) -integer*1 message77(77) +integer*1 nxor(N),hdec(N) +integer*1 message91(91) integer nrw(M),ncw integer Nm(7,M) integer Mn(3,N) ! 3 checks per bit @@ -35,7 +34,6 @@ enddo ncnt=0 nclast=0 -maxsuper=5 maxiterations=1 zsum=0.0 @@ -63,10 +61,12 @@ zsum=zsum+zn ! write(*,*) 'number of unsatisfied parity checks ',ncheck if( ncheck .eq. 0 ) then ! we have a codeword - if crc is good, return it decoded=cw(1:K) - call chkcrc14a(decoded,nbadcrc) + call get_crc14(decoded,91,nbadcrc) +write(*,*) nbadcrc +write(*,'(91i1)') decoded nharderror=count( (2*cw-1)*llr .lt. 0.0 ) if(nbadcrc.eq.0) then - message77=decoded(1:77) + message91=decoded(1:91) dmin=0.0 return endif @@ -118,13 +118,18 @@ dmin=0.0 enddo enddo ! bp iterations - ndeep=1 -llr=zsum - call osd174_91(llr,apmask,ndeep,message77,cw,nharderror,dmin) + + call osd174_91(zsum,Keff,apmask,ndeep,message91,cw,nharderror,dminosd) if(nharderror.gt.0) then + hdec=0 + where(llr .ge. 0) hdec=1 + nxor=ieor(hdec,cw) + dmin=sum(nxor*abs(llr)) return endif enddo ! super iterations + nharderror=-1 + return end subroutine decode174_91 diff --git a/lib/ft8/encode174_91_nocrc.f90 b/lib/ft8/encode174_91_nocrc.f90 new file mode 100644 index 000000000..6d583d907 --- /dev/null +++ b/lib/ft8/encode174_91_nocrc.f90 @@ -0,0 +1,49 @@ +subroutine encode174_91_nocrc(message,codeword) +! +! Add a 14-bit CRC to a 77-bit message and return a 174-bit codeword +! +use, intrinsic :: iso_c_binding +use iso_c_binding, only: c_loc,c_size_t +use crc + +integer, parameter:: N=174, K=91, M=N-K +character*91 tmpchar +integer*1 codeword(N) +integer*1 gen(M,K) +integer*1 message(K) +integer*1 pchecks(M) +integer*1, target :: i1MsgBytes(12) +include "ldpc_174_91_c_generator.f90" +logical first +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=3 + do jj=1, ibmax + icol=(j-1)*4+jj + if( btest(istr,4-jj) ) gen(i,icol)=1 + enddo + enddo + enddo +first=.false. +endif + +do i=1,M + nsum=0 + do j=1,K + nsum=nsum+message(j)*gen(i,j) + enddo + pchecks(i)=mod(nsum,2) +enddo + +codeword(1:K)=message +codeword(K+1:N)=pchecks + +return +end subroutine encode174_91_nocrc diff --git a/lib/ft8/ft8b.f90 b/lib/ft8/ft8b.f90 index 1f77878ae..638b75557 100644 --- a/lib/ft8/ft8b.f90 +++ b/lib/ft8/ft8b.f90 @@ -18,7 +18,7 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,nzhsym,lapon, & real llra(174),llrb(174),llrc(174),llrd(174),llrz(174) !Soft symbols real dd0(15*12000) real ss(9) - integer*1 message77(77),apmask(174),cw(174) + integer*1 message77(77),message91(91),apmask(174),cw(174) integer apsym(58),aph10(10) integer mcq(29),mcqru(29),mcqfd(29),mcqtest(29),mcqww(29) integer mrrr(19),m73(19),mrr73(19) @@ -407,12 +407,19 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,nzhsym,lapon, & dmin=0.0 if(nharderrors.lt.0 .and. ncheck.le.30 .and. ndepth.ge.2) then ndeep=ndepth +!ndeep=ndepth-1 if(abs(nfqso-f1).le.napwid .or. abs(nftx-f1).le.napwid .or. ncontest.eq.7) then ndeep=4 +!ndeep=3 endif if(nagain) ndeep=5 +!if(nagain) ndeep=4 call timer('osd174_91 ',0) - call osd174_91(llrz,apmask,ndeep,message77,cw,nharderrors,dmin) + Keff=91 +! call osd174_91(llrz,Keff,apmask,ndeep,message91,cw,nharderrors,dmin) + maxsuper=2 + call decode174_91(llrz,Keff,ndeep,apmask,maxsuper,message91,cw,nharderrors,niterations,ncheck,dmin) + if(nharderrors.ge.0) message77=message91(1:77) call timer('osd174_91 ',1) endif diff --git a/lib/ft8/osd174_91.f90 b/lib/ft8/osd174_91.f90 index bea351101..d9584f56e 100644 --- a/lib/ft8/osd174_91.f90 +++ b/lib/ft8/osd174_91.f90 @@ -1,237 +1,269 @@ -subroutine osd174_91(llr,apmask,ndeep,message77,cw,nhardmin,dmin) +subroutine osd174_91(llr,k,apmask,ndeep,message91,cw,nhardmin,dmin) ! ! An ordered-statistics decoder for the (174,91) code. +! Message payload is 77 bits. Any or all of a 14-bit CRC can be +! used for detecting incorrect codewords. The remaining CRC bits are +! cascaded with the LDPC code for the purpose of improving the +! distance spectrum of the code. +! +! If p1 (0.le.p1.le.14) is the number of CRC14 bits that are +! to be used for bad codeword detection, then the argument k should +! be set to 77+p1. +! +! Valid values for k are in the range [77,91]. +! + character*14 c14 + integer, parameter:: N=174 + integer*1 apmask(N),apmaskr(N) + integer*1, allocatable, save :: gen(:,:) + integer*1, allocatable :: genmrb(:,:),g2(:,:) + integer*1, allocatable :: temp(:),m0(:),me(:),mi(:),misub(:),e2sub(:),e2(:),ui(:) + integer*1, allocatable :: r2pat(:) + integer indices(N),nxor(N) + integer*1 cw(N),ce(N),c0(N),hdec(N) + integer*1, allocatable :: decoded(:) + integer*1 message91(91),m96(96) + integer indx(N) + real llr(N),rx(N),absrx(N) + + logical first,reset + data first/.true./ + save first + + allocate( genmrb(k,N), g2(N,k) ) + allocate( temp(k), m0(k), me(k), mi(k), misub(k), e2sub(N-k), e2(N-k), ui(N-k) ) + allocate( r2pat(N-k), decoded(k) ) + + if( first ) then ! fill the generator matrix +! +! Create generator matrix for partial CRC cascaded with LDPC code. ! -integer, parameter:: N=174, K=91, 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_174_91_c_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=3 - do jj=1, ibmax - irow=(j-1)*4+jj - if( btest(istr,4-jj) ) gen(irow,K+i)=1 +! Let p2=91-k and p1+p2=14. +! +! The last p2 bits of the CRC14 are cascaded with the LDPC code. +! +! The first p1=k-77 CRC14 bits will be used for error detection. +! + allocate( gen(k,N) ) + gen=0 + do i=1,k + message91=0 + message91(i)=1 + if(i.le.77) then + m96=0 + m96(1:91)=message91 + call get_crc14(m96,96,ncrc14) + write(c14,'(b14.14)') ncrc14 + read(c14,'(14i1)') message91(78:91) + message91(78:k)=0 + endif + call encode174_91_nocrc(message91,cw) + gen(i,:)=cw enddo - enddo - enddo - do irow=1,K - gen(irow,irow)=1 - enddo -first=.false. -endif -rx=llr -apmaskr=apmask + first=.false. + endif + + rx=llr + apmaskr=apmask ! Hard decisions on the received word. -hdec=0 -where(rx .ge. 0) hdec=1 + 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) + 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 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 +! 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 - exit - endif - enddo -enddo + enddo -g2=transpose(genmrb) + 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. +! 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. +! 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) + 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 mrbencode91(m0,c0,g2,N,K) -nxor=ieor(c0,hdec) -nhardmin=sum(nxor) -dmin=sum(nxor*absrx) + call mrbencode91(m0,c0,g2,N,k) + nxor=ieor(c0,hdec) + nhardmin=sum(nxor) + dmin=sum(nxor*absrx) -cw=c0 -ntotal=0 -nrejected=0 -npre1=0 -npre2=0 - -if(ndeep.eq.0) goto 998 ! norder=0 -if(ndeep.gt.5) ndeep=5 -if( ndeep.eq. 1) then - nord=1 + cw=c0 + ntotal=0 + nrejected=0 npre1=0 npre2=0 - nt=40 - ntheta=12 -elseif(ndeep.eq.2) then - nord=1 - 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 - 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 - 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 - d1=0. - 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 mrbencode91(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)) + if(ndeep.eq.0) goto 998 ! norder=0 + if(ndeep.gt.6) ndeep=6 + if( ndeep.eq. 1) then + nord=1 + npre1=0 + npre2=0 + nt=40 + ntheta=12 + elseif(ndeep.eq.2) then + nord=1 + npre1=1 + npre2=0 + nt=40 +! ntheta=12 + ntheta=10 + elseif(ndeep.eq.3) then + nord=1 + npre1=1 + npre2=1 + nt=40 + ntheta=12 + ntau=14 + elseif(ndeep.eq.4) then + nord=2 + npre1=1 + npre2=1 + nt=40 + ntheta=12 + ntau=17 + elseif(ndeep.eq.5) then + nord=3 + npre1=1 + npre2=1 + nt=40 + ntheta=12 + ntau=15 + elseif(ndeep.eq.6) then + nord=4 + npre1=1 + npre2=1 + nt=95 + ntheta=12 + ntau=15 + 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 - e2=ieor(e2sub,g2(K+1:N,n1)) - nd1Kpt=sum(e2(1:nt))+2 + iend=1 endif - if(nd1Kpt .le. ntheta) then - call mrbencode91(me,ce,g2,N,K) - nxor=ieor(ce,hdec) + d1=0. + 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 - dd=d1+sum(e2sub*absrx(K+1:N)) + call mrbencode91(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 - dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N)) + e2=ieor(e2sub,g2(k+1:N,n1)) + nd1kpt=sum(e2(1:nt))+2 endif - if( dd .lt. dmin ) then - dmin=dd - cw=ce - nhardmin=sum(nxor) - nd1Kptbest=nd1Kpt + if(nd1kpt .le. ntheta) then + call mrbencode91(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 - else - nrejected=nrejected+1 - endif - enddo + enddo ! Get the next test error pattern, iflag will go negative ! when the last pattern with weight iorder has been generated. - call nextpat91(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 boxit91(reset,mi(1:ntau),ntau,ntotal,i1,i2) + call nextpat91(misub,k,iorder,iflag) enddo enddo - ncount2=0 - ntotal2=0 - reset=.true. + 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 boxit91(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 mrbencode91(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 + 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 mrbencode91(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 fetchit91(reset,r2pat(1:ntau),ntau,in1,in2) if(in1.gt.0.and.in2.gt.0) then ncount2=ncount2+1 - mi=misub + 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 + if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:k),mi).eq.1)) cycle me=ieor(m0,mi) - call mrbencode91(me,ce,g2,N,K) + call mrbencode91(me,ce,g2,N,k) nxor=ieor(ce,hdec) dd=sum(nxor*absrx) if( dd .lt. dmin ) then @@ -240,136 +272,138 @@ if(npre2.eq.1) then nhardmin=sum(nxor) endif goto 778 - endif + endif + enddo + call nextpat91(misub,k,nord,iflag) enddo - call nextpat91(misub,K,nord,iflag) - enddo -endif + endif 998 continue -! Re-order the codeword to [message bits][parity bits] format. -cw(indices)=cw -hdec(indices)=hdec -decoded=cw(1:K) -call chkcrc14a(decoded,nbadcrc) -message77=decoded(1:77) -if(nbadcrc.eq.1) nhardmin=-nhardmin +! Re-order the codeword to [message bits][parity bits] format. + cw(indices)=cw + hdec(indices)=hdec + message91=cw(1:91) + m96=0 + m96(1:77)=cw(1:77) + m96(83:96)=cw(78:91) + call get_crc14(m96,96,nbadcrc) + if(nbadcrc.ne.0) nhardmin=-nhardmin -return + return end subroutine osd174_91 subroutine mrbencode91(me,codeword,g2,N,K) -integer*1 me(K),codeword(N),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 + codeword=0 + do i=1,K + if( me(i) .eq. 1 ) then + codeword=ieor(codeword,g2(1:N,i)) + endif + enddo + return end subroutine mrbencode91 subroutine nextpat91(mi,k,iorder,iflag) - integer*1 mi(k),ms(k) + 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 + 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 nextpat91 subroutine boxit91(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 + 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 + if(reset) then + patterns=-1 + fp=-1 + np=-1 + sc=-1 + indexes=-1 + reset=.false. + endif - 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 + 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 boxit91 subroutine fetchit91(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 + 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 + 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) + 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 + 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 fetchit91 - +