1. Improvements to the OSD to allow deeper wideband decoding. 2. Add a third decoding pass. 3. Change symbol metric from max-log to max-amplitude.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8031 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2017-08-22 00:14:51 +00:00
parent 55b74f5a0e
commit d5d8abc29a
4 changed files with 238 additions and 85 deletions

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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)