w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-10 06:22:15 -05:00
Code Issues Packages Projects Releases Wiki Activity

Implements decoding of FST4W messages as (240,50) crc-less codewords.

Browse Source 
By cascading the full 24-bit crc generator matrix with the (240,74) LDPC code generator, create a (240,50) generator that is used to decode with approximately 1 dB better sensitivity than the (240,64) with 14-bit CRC approach that is normally used. This approach treats the CRC bits as additional parity bits and provides no means for identifying incorrect codewords. All codewords on the list generated by the OSD algorithm have CRCs that match the CRC of the message payload. Codewords are validated by unpacking the message and comparing the unpacked message with the list of stored callsign/grid pairs stored in the fst4w_calls.txt file.
This commit is contained in:
Steven Franke 2020-12-13 15:33:12 -06:00
parent 5b60ddada1
commit 535c02d900
5 changed files with 452 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -570,6 +570,7 @@ set (wsjt_FSRCS
lib/fst4/ldpcsim240_74.f90
lib/fst4/osd240_101.f90
lib/fst4/osd240_74.f90
lib/fst4/fastosd240_74.f90
lib/fst4/get_crc24.f90
lib/fst4/fst4_baseline.f90
)

3
lib/fst4/decode240_74.f90
View File

@ -135,7 +135,8 @@ subroutine decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw,ntype,nharder
do i=1,nosd
zn=zsave(:,i)
call osd240_74(zn,Keff,apmask,norder,message74,cw,nharderror,dminosd)
! call osd240_74(zn,Keff,apmask,norder,message74,cw,nharderror,dminosd)
call fastosd240_74(zn,Keff,apmask,norder,message74,cw,nharderror,dminosd)
if(nharderror.gt.0) then
hdec=0
where(llr .ge. 0) hdec=1

291
lib/fst4/fastosd240_74.f90 Normal file
View File

@ -0,0 +1,291 @@
subroutine fastosd240_74(llr,k,apmask,ndeep,message74,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (240,74) code.
! Message payload is 50 bits. Any or all of a 24-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.24) is the number of CRC24 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 [50,74].
!
character*24 c24
integer, parameter:: N=240
integer*1 apmask(N),apmaskr(N)
integer*1, allocatable, save :: gen(:,:)
integer*1, allocatable :: genmrb(:,:),g2(:,:)
integer*1, allocatable :: temp(:),temprow(:),m0(:),me(:),mi(:)
integer indices(N),indices2(N),nxor(N)
integer*1 cw(N),ce(N),c0(N),hdec(N)
integer*1, allocatable :: decoded(:)
integer*1 message74(74)
integer*1, allocatable :: sp(:)
integer indx(N),ksave
real llr(N),rx(N),absrx(N)
logical first
data first/.true./,ksave/64/
save first,ksave
allocate( genmrb(k,N), g2(N,k) )
allocate( temp(k), temprow(n), m0(k), me(k), mi(k) )
allocate( decoded(k) )
if( first .or. k.ne.ksave) then ! fill the generator matrix
!
! Create generator matrix for partial CRC cascaded with LDPC code.
!
! Let p2=74-k and p1+p2=24.
!
! The last p2 bits of the CRC24 are cascaded with the LDPC code.
!
! The first p1=k-50 CRC24 bits will be used for error detection.
!
if( allocated(gen) ) deallocate(gen)
allocate( gen(k,N) )
gen=0
do i=1,k
message74=0
message74(i)=1
if(i.le.50) then
call get_crc24(message74,74,ncrc24)
write(c24,'(b24.24)') ncrc24
read(c24,'(24i1)') message74(51:74)
message74(51:k)=0
endif
call encode240_74(message74,cw)
gen(i,:)=cw
enddo
first=.false.
ksave=k
endif
! Use best k elements from the sorted list for the first basis. For the 2nd basis replace
! the nswap lowest quality symbols with the best nswap elements from the parity symbols.
nswap=20
do ibasis=1,2
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(llr)
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
if(ibasis.eq.2) then
do i=k-nswap+1,k
temp(1:k)=genmrb(1:k,i)
genmrb(1:k,i)=genmrb(1:k,i+nswap)
genmrb(1:k,i+nswap)=temp(1:k)
itmp=indices(i)
indices(i)=indices(i+nswap)
indices(i+nswap)=itmp
enddo
endif
! 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).
icol=1
indices2=0
nskipped=0
do id=1,k
iflag=0
do while(iflag.eq.0)
if(genmrb(id,icol).ne.1) then
do j=id+1,k
if(genmrb(j,icol).eq.1) then
temprow=genmrb(id,:)
genmrb(id,:)=genmrb(j,:)
genmrb(j,:)=temprow
iflag=1
endif
enddo
if(iflag.eq.0) then ! skip this column
nskipped=nskipped+1
indices2(k+nskipped)=icol ! put icol where skipped columns go
icol=icol+1 ! look at the next column
endif
else
iflag=1
endif
enddo
indices2(id)=icol
do j=1,k
if(id.ne.j .and. genmrb(j,icol).eq.1) then
genmrb(j,:)=ieor(genmrb(id,:),genmrb(j,:))
endif
enddo
icol=icol+1
enddo
do i=k+nskipped+1,240
indices2(i)=i
enddo
genmrb(1:k,:)=genmrb(1:k,indices2)
indices=indices(indices2)
!************************************
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=abs(llr)
absrx=absrx(indices)
rx=rx(indices)
apmaskr=apmaskr(indices)
call mrbencode74(m0,c0,g2,N,k)
nxor=ieor(c0,hdec)
nhardmin=sum(nxor)
dmin=sum(nxor*absrx)
np=32
if(ibasis.eq.1) allocate(sp(np))
cw=c0
ntotal=0
nrejected=0
if(ndeep.eq.0) goto 998 ! norder=0
if(ndeep.gt.4) ndeep=4
if( ndeep.eq. 1) then
nord=1
xlambda=0.0
nsyncmax=np
elseif(ndeep.eq.2) then
nord=2
xlambda=0.0
nsyncmax=np
elseif(ndeep.eq.3) then
nord=3
xlambda=4.0
nsyncmax=11
elseif(ndeep.eq.4) then
nord=4
xlambda=3.5
nsyndmax=11
endif
s1=sum(absrx(1:k))
s2=sum(absrx(k+1:N))
rho=s1/(s1+xlambda*s2)
rhodmin=rho*dmin
nerr64=-1
do iorder=1,nord
!beta=0.0
!if(iorder.ge.3) beta=0.4
!spnc_order=sum(absrx(k-iorder+1:k))+beta*(N-k)
!if(dmin.lt.spnc_order) cycle
mi(1:k-iorder)=0
mi(k-iorder+1:k)=1
iflag=k-iorder+1
do while(iflag .ge.0)
ntotal=ntotal+1
me=ieor(m0,mi)
d1=sum(mi(1:k)*absrx(1:k))
if(d1.gt.rhodmin) exit
call partial_syndrome(me,sp,np,g2,N,K)
nwhsp=sum(ieor(sp(1:np),hdec(k:k+np-1)))
if(nwhsp.le.nsyndmax) then
call mrbencode74(me,ce,g2,N,k)
nxor=ieor(ce,hdec)
dd=sum(nxor*absrx(1:N))
if( dd .lt. dmin ) then
dmin=dd
rhodmin=rho*dmin
cw=ce
nhardmin=sum(nxor)
nwhspmin=nwhsp
nerr64=sum(nxor(1:K))
endif
endif
! Get the next test error pattern, iflag will go negative
! when the last pattern with weight iorder has been generated.
call nextpat74(mi,k,iorder,iflag)
enddo
enddo
998 continue
! Re-order the codeword to [message bits][parity bits] format.
cw(indices)=cw
hdec(indices)=hdec
message74=cw(1:74)
call get_crc24(message74,74,nbadcrc)
if(nbadcrc.eq.0) exit
nhardmin=-nhardmin
enddo ! basis loop
return
end subroutine fastosd240_74
subroutine mrbencode74(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 mrbencode74
subroutine partial_syndrome(me,sp,np,g2,N,K)
integer*1 me(K),sp(np),g2(N,K)
! compute partial syndrome
sp=0
do i=1,K
if( me(i) .eq. 1 ) then
sp=ieor(sp,g2(K:K+np-1,i))
endif
enddo
return
end subroutine partial_syndrome
subroutine nextpat74(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 nextpat74

2
lib/fst4/ldpcsim240_74.f90
View File

@ -101,7 +101,7 @@ write(*,'(24i1)') msgbits(51:74)
llr=2.0*rxdata/(ss*ss)
apmask=0
dmin=0.0
maxosd=0
maxosd=2
call decode240_74(llr, Keff, maxosd, norder, apmask, message74, cw, ntype, nharderror, dmin)
if(nharderror.ge.0) then
n2err=0

232
lib/fst4_decode.f90
View File

@ -33,15 +33,17 @@ contains
ndepth,ntrperiod,nexp_decode,ntol,emedelay,lagain,lapcqonly,mycall, &
hiscall,iwspr)
use prog_args
use timer_module, only: timer
use packjt77
use, intrinsic :: iso_c_binding
include 'fst4/fst4_params.f90'
parameter (MAXCAND=100)
parameter (MAXCAND=100,MAXWCALLS=100)
class(fst4_decoder), intent(inout) :: this
procedure(fst4_decode_callback) :: callback
character*37 decodes(100)
character*37 msg,msgsent
character*20 wcalls(MAXWCALLS), wpart
character*77 c77
character*12 mycall,hiscall
character*12 mycall0,hiscall0
@ -66,7 +68,8 @@ contains
integer mcq(29),mrrr(19),m73(19),mrr73(19)
logical badsync,unpk77_success,single_decode
logical first,nohiscall,lwspr,ex
logical first,nohiscall,lwspr,ex,wcalls_exists,donocrcdecode
logical new_callsign
integer*2 iwave(30*60*12000)
@ -80,6 +83,7 @@ contains
0,1,0,1,0,1,1,0,1,1,1,1,1,0,0,0,1,0,1/
data first/.true./,hmod/1/
save first,apbits,nappasses,naptypes,mycall0,hiscall0
save wcalls,nwcalls
this%callback => callback
dxcall13=hiscall ! initialize for use in packjt77
@ -88,6 +92,24 @@ contains
if(iwspr.ne.0.and.iwspr.ne.1) return
if(first) then
! read the fst4_calls.txt file
write(*,*) 'data_dir is:',trim(data_dir)
inquire(file=trim(data_dir)//'/fst4w_calls.txt',exist=wcalls_exists)
if( wcalls_exists ) then
open(42,file=trim(data_dir)//'/fst4w_calls.txt',status='unknown')
write(*,*) 'fst4w_calls.txt exists'
do i=1,MAXWCALLS
wcalls(i)=''
read(42,fmt='(a)',end=2867) wcalls(i)
wcalls(i)=adjustl(wcalls(i))
if(len(trim(wcalls(i))).eq.0) exit
write(*,*) 'record ',i,':',wcalls(i),':',len(trim(wcalls(i)))
enddo
2867 nwcalls=i-1
write(*,*) 'nwcalls ',nwcalls
close(42)
endif
mcq=2*mod(mcq+rvec(1:29),2)-1
mrrr=2*mod(mrrr+rvec(59:77),2)-1
m73=2*mod(m73+rvec(59:77),2)-1
@ -216,12 +238,15 @@ contains
if(ndepth.eq.3) then
nblock=4
jittermax=2
donocrcdecode=.true.
elseif(ndepth.eq.2) then
nblock=3
jittermax=0
nblock=4
jittermax=2
donocrcdecode=.false.
elseif(ndepth.eq.1) then
nblock=1
nblock=4
jittermax=0
donocrcdecode=.false.
endif
ndropmax=1
@ -244,7 +269,7 @@ contains
if(iwspr.eq.1) then !FST4W
!300 Hz wide noise-fit window
nfa=max(100,nint(nfqso+1.5*baud-150))
nfa=max(100,nint(nfqso+1.5*baud-150))
nfb=min(4800,nint(nfqso+1.5*baud+150))
fa=max(100,nint(nfqso+1.5*baud-ntol)) ! signal search window
fb=min(4800,nint(nfqso+1.5*baud+ntol))
@ -252,18 +277,19 @@ contains
fa=max(100,nint(nfa+1.5*baud))
fb=min(4800,nint(nfb+1.5*baud))
! extend noise fit 100 Hz outside of search window
nfa=max(100,nfa-100)
nfa=max(100,nfa-100)
nfb=min(4800,nfb+100)
else
fa=max(100,nint(nfa+1.5*baud))
fb=min(4800,nint(nfb+1.5*baud))
! extend noise fit 100 Hz outside of search window
nfa=max(100,nfa-100)
nfa=max(100,nfa-100)
nfb=min(4800,nfb+100)
endif
ndecodes=0
decodes=' '
new_callsign=.false.
do inb=0,inb1,inb2
if(nb.lt.0) npct=inb
call blanker(iwave,nfft1,ndropmax,npct,c_bigfft)
@ -275,16 +301,16 @@ contains
nsyncoh=8
minsync=1.20
if(ntrperiod.eq.15) minsync=1.15
! Get first approximation of candidate frequencies
call get_candidates_fst4(c_bigfft,nfft1,nsps,hmod,fs,fa,fb,nfa,nfb, &
minsync,ncand,candidates0)
minsync,ncand,candidates0)
isbest=0
fc2=0.
do icand=1,ncand
fc0=candidates0(icand,1)
if(iwspr.eq.0 .and. nb.lt.0 .and. npct.ne.0 .and. &
abs(fc0-(nfqso+1.5*baud)).gt.ntol) cycle
abs(fc0-(nfqso+1.5*baud)).gt.ntol) cycle
detmet=candidates0(icand,2)
! Downconvert and downsample a slice of the spectrum centered on the
@ -330,7 +356,7 @@ contains
endif
enddo
ncand=ic
! If FST4 and Single Decode is not checked, then find candidates within
! 20 Hz of nfqso and put them at the top of the list
if(iwspr.eq.0 .and. .not.single_decode) then
@ -373,7 +399,7 @@ contains
bitmetrics=0
call timer('bitmetrc',0)
call get_fst4_bitmetrics(cframe,nss,nblock,nhicoh,bitmetrics, &
s4,nsync_qual,badsync)
s4,nsync_qual,badsync)
call timer('bitmetrc',1)
if(badsync) cycle
@ -405,7 +431,7 @@ contains
iaptype=0
endif
if(itry.gt.nblock) then ! do ap passes
if(itry.gt.nblock .and. iwspr.eq.0) then ! do ap passes
llr=llrs(:,nblock) ! Use largest blocksize as the basis for AP passes
iaptype=naptypes(nQSOProgress,itry-nblock)
if(lapcqonly) iaptype=1
@ -428,7 +454,7 @@ contains
apmask(1:58)=1
llr(1:58)=apmag*apbits(1:58)
endif
if(iaptype.eq.4 .or. iaptype.eq.5 .or. iaptype .eq.6) then
apmask=0
apmask(1:77)=1
@ -438,7 +464,7 @@ contains
if(iaptype.eq.6) llr(59:77)=apmag*mrr73(1:19)
endif
endif
dmin=0.0
nharderrors=-1
unpk77_success=.false.
@ -448,83 +474,133 @@ contains
norder=3
call timer('d240_101',0)
call decode240_101(llr,Keff,maxosd,norder,apmask,message101, &
cw,ntype,nharderrors,dmin)
cw,ntype,nharderrors,dmin)
call timer('d240_101',1)
elseif(iwspr.eq.1) then
maxosd=2
call timer('d240_74 ',0)
Keff=64
norder=4
call decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw, &
ntype,nharderrors,dmin)
call timer('d240_74 ',1)
endif
if(nharderrors .ge.0) then
if(count(cw.eq.1).eq.0) then
nharderrors=-nharderrors
cycle
endif
if(iwspr.eq.0) then
write(c77,'(77i1)') mod(message101(1:77)+rvec,2)
write(c77,'(77i1)') mod(message101(1:77)+rvec,2)
call unpack77(c77,1,msg,unpk77_success)
elseif(iwspr.eq.1) then
iaptype=0
if( donocrcdecode ) then
maxosd=1
call timer('d240_74 ',0)
Keff=50
norder=4
call decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw, &
ntype,nharderrors,dmin)
call timer('d240_74 ',1)
if(count(cw.eq.1).eq.0) then
nharderrors=-nharderrors
cycle
endif
write(c77,'(50i1)') message74(1:50)
c77(51:77)='000000000000000000000110000'
call unpack77(c77,1,msg,unpk77_success)
else
if(unpk77_success) then
unpk77_success=.false.
do i=1,nwcalls
if(index(msg,trim(wcalls(i))).gt.0) then
unpk77_success=.true.
iaptype=8
endif
enddo
endif
endif
if(.not. unpk77_success) then
maxosd=2
call timer('d240_74 ',0)
Keff=64
norder=3
call decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw, &
ntype,nharderrors,dmin)
call timer('d240_74 ',1)
if(nharderrors.lt.0) cycle
if(count(cw.eq.1).eq.0) then
nharderrors=-nharderrors
cycle
endif
write(c77,'(50i1)') message74(1:50)
c77(51:77)='000000000000000000000110000'
call unpack77(c77,1,msg,unpk77_success)
endif
if(unpk77_success) then
idupe=0
do i=1,ndecodes
if(decodes(i).eq.msg) idupe=1
if(unpk77_success.and.Keff.eq.64) then
i1=index(msg,' ')
i2=i1+index(msg(i1+1:),' ')
wpart=trim(msg(1:i2))
ifound=0
do i=1,nwcalls
if(index(wcalls(i),wpart).ne.0) ifound=1
enddo
if(idupe.eq.1) goto 800
ndecodes=ndecodes+1
decodes(ndecodes)=msg
if(iwspr.eq.0) then
call get_fst4_tones_from_bits(message101,itone,0)
else
call get_fst4_tones_from_bits(message74,itone,1)
endif
inquire(file='plotspec',exist=ex)
fmid=-999.0
call timer('dopsprd ',0)
if(ex) then
call dopspread(itone,iwave,nsps,nmax,ndown,hmod, &
isbest,fc_synced,fmid,w50)
endif
call timer('dopsprd ',1)
xsig=0
do i=1,NN
xsig=xsig+s4(itone(i),i)
enddo
base=candidates(icand,5)
arg=600.0*(xsig/base)-1.0
if(arg.gt.0.0) then
xsnr=10*log10(arg)-35.5-12.5*log10(nsps/8200.0)
if(ntrperiod.eq. 15) xsnr=xsnr+2
if(ntrperiod.eq. 30) xsnr=xsnr+1
if(ntrperiod.eq. 900) xsnr=xsnr+1
if(ntrperiod.eq.1800) xsnr=xsnr+2
else
xsnr=-99.9
if(ifound.eq.0) then ! This is a new callsign
new_callsign=.true.
if(nwcalls.lt.MAXWCALLS) then
nwcalls=nwcalls+1
wcalls(nwcalls)=wpart
else
wcalls(1:nwcalls-1)=wcalls(2:nwcalls)
wcalls(nwcalls)=wpart
endif
endif
endif
endif
if(nharderrors .ge.0 .and. unpk77_success) then
idupe=0
do i=1,ndecodes
if(decodes(i).eq.msg) idupe=1
enddo
if(idupe.eq.1) goto 800
ndecodes=ndecodes+1
decodes(ndecodes)=msg
if(iwspr.eq.0) then
call get_fst4_tones_from_bits(message101,itone,0)
else
cycle
call get_fst4_tones_from_bits(message74,itone,1)
endif
inquire(file='plotspec',exist=ex)
fmid=-999.0
call timer('dopsprd ',0)
if(ex) then
call dopspread(itone,iwave,nsps,nmax,ndown,hmod, &
isbest,fc_synced,fmid,w50)
endif
call timer('dopsprd ',1)
xsig=0
do i=1,NN
xsig=xsig+s4(itone(i),i)
enddo
base=candidates(icand,5)
arg=600.0*(xsig/base)-1.0
if(arg.gt.0.0) then
xsnr=10*log10(arg)-35.5-12.5*log10(nsps/8200.0)
if(ntrperiod.eq. 15) xsnr=xsnr+2
if(ntrperiod.eq. 30) xsnr=xsnr+1
if(ntrperiod.eq. 900) xsnr=xsnr+1
if(ntrperiod.eq.1800) xsnr=xsnr+2
else
xsnr=-99.9
endif
nsnr=nint(xsnr)
qual=0.
qual=0.0
if(iaptype.eq.8) qual=1.
fsig=fc_synced - 1.5*baud
if(ex) then
write(21,3021) nutc,icand,itry,nsyncoh,iaptype, &
ijitter,ntype,nsync_qual,nharderrors,dmin, &
sync,xsnr,xdt,fsig,w50,trim(msg)
ijitter,ntype,nsync_qual,nharderrors,dmin, &
sync,xsnr,xdt,fsig,w50,trim(msg)
3021 format(i6.6,6i3,2i4,f6.1,f7.2,f6.1,f6.2,f7.1,f7.3,1x,a)
flush(21)
endif
call this%callback(nutc,smax1,nsnr,xdt,fsig,msg, &
iaptype,qual,ntrperiod,lwspr,fmid,w50)
iaptype,qual,ntrperiod,lwspr,fmid,w50)
if(iwspr.eq.0 .and. nb.lt.0) go to 900
goto 800
endif
@ -532,9 +608,15 @@ contains
enddo ! istart jitter
800 enddo !candidate list
enddo ! noise blanker loop
if(new_callsign) then ! re-write the fst4w_calls.txt file
open(42,file=trim(data_dir)//'/fst4w_calls.txt',status='unknown')
do i=1,nwcalls
write(42,'(a20)') trim(wcalls(i))
enddo
close(42)
endif
900 return
end subroutine decode
end subroutine decode
subroutine sync_fst4(cd0,i0,f0,hmod,ncoh,np,nss,ntr,fs,sync)
@ -902,6 +984,6 @@ contains
enddo
return
end subroutine dopspread
end subroutine dopspread
end module fst4_decode