Working toward whole-hog switchover to new LDPC code with 77 bit messages in MSK144.

This commit is contained in:
Steve Franke 2018-07-07 09:44:52 -05:00
parent 42f75bf404
commit 622ed4a3ab
6 changed files with 405 additions and 115 deletions

View File

@ -8,18 +8,18 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter)
integer*1 cw(N),apmask(N) integer*1 cw(N),apmask(N)
integer*1 decoded(K) integer*1 decoded(K)
integer*1 message77(77) integer*1 message77(77)
integer Nm(12,M) integer Nm(11,M)
integer Mn(4,N) integer Mn(3,N)
integer nrw(M),ncw(N) integer nrw(M)
integer synd(M) integer synd(M)
real tov(4,N) real tov(4,N)
real toc(12,M) real toc(11,M)
real tanhtoc(12,M) real tanhtoc(11,M)
real zn(N) real zn(N)
real llr(N) real llr(N)
real Tmn real Tmn
include "ldpc_128_90_b_reordered_parity.f90" include "ldpc_128_90_reordered_parity.f90"
decoded=0 decoded=0
toc=0 toc=0
@ -39,7 +39,7 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter)
! Update bit log likelihood ratios (tov=0 in iteration 0). ! Update bit log likelihood ratios (tov=0 in iteration 0).
do i=1,N do i=1,N
if( apmask(i) .ne. 1 ) then if( apmask(i) .ne. 1 ) then
zn(i)=llr(i)+sum(tov(1:ncw(i),i)) zn(i)=llr(i)+sum(tov(1:ncw,i))
else else
zn(i)=llr(i) zn(i)=llr(i)
endif endif
@ -74,7 +74,7 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter)
ncnt=ncnt+1 ncnt=ncnt+1
endif endif
! write(*,*) iter,ncheck,nd,ncnt ! write(*,*) iter,ncheck,nd,ncnt
if( ncnt .ge. 5 .and. iter .ge. 10 .and. ncheck .gt. 15) then if( ncnt .ge. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) then
nharderror=-1 nharderror=-1
return return
endif endif
@ -86,7 +86,7 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter)
do i=1,nrw(j) do i=1,nrw(j)
ibj=Nm(i,j) ibj=Nm(i,j)
toc(i,j)=zn(ibj) toc(i,j)=zn(ibj)
do kk=1,4 ! subtract off what the bit had received from the check do kk=1,ncw ! subtract off what the bit had received from the check
if( Mn(kk,ibj) .eq. j ) then if( Mn(kk,ibj) .eq. j ) then
toc(i,j)=toc(i,j)-tov(kk,ibj) toc(i,j)=toc(i,j)-tov(kk,ibj)
endif endif
@ -96,11 +96,11 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter)
! send messages from check nodes to variable nodes ! send messages from check nodes to variable nodes
do i=1,M do i=1,M
tanhtoc(1:12,i)=tanh(-toc(1:12,i)/2) tanhtoc(1:11,i)=tanh(-toc(1:11,i)/2)
enddo enddo
do j=1,N do j=1,N
do i=1,ncw(j) do i=1,ncw
ichk=Mn(i,j) ! Mn(:,j) are the checks that include bit j ichk=Mn(i,j) ! Mn(:,j) are the checks that include bit j
Tmn=product(tanhtoc(1:nrw(ichk),ichk),mask=Nm(1:nrw(ichk),ichk).ne.j) Tmn=product(tanhtoc(1:nrw(ichk),ichk),mask=Nm(1:nrw(ichk),ichk).ne.j)
call platanh(-Tmn,y) call platanh(-Tmn,y)

116
lib/bpdecode128_90.f90.save Normal file
View File

@ -0,0 +1,116 @@
subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter)
!
! A log-domain belief propagation decoder for the (128,90) code.
!
use iso_c_binding, only: c_loc,c_size_t
use crc
integer, parameter:: N=128, K=90, M=N-K
integer*1 cw(N),apmask(N)
integer*1 decoded(K)
integer*1 message77(77)
integer Nm(12,M)
integer Mn(4,N)
integer nrw(M),ncw(N)
integer synd(M)
real tov(4,N)
real toc(12,M)
real tanhtoc(12,M)
real zn(N)
real llr(N)
real Tmn
include "ldpc_128_90_b_reordered_parity.f90"
decoded=0
toc=0
tov=0
tanhtoc=0
! initialize messages to checks
do j=1,M
do i=1,nrw(j)
toc(i,j)=llr((Nm(i,j)))
enddo
enddo
ncnt=0
do iter=0,maxiterations
! Update bit log likelihood ratios (tov=0 in iteration 0).
do i=1,N
if( apmask(i) .ne. 1 ) then
zn(i)=llr(i)+sum(tov(1:ncw(i),i))
else
zn(i)=llr(i)
endif
enddo
! Check to see if we have a codeword (check before we do any iteration).
cw=0
where( zn .gt. 0. ) cw=1
ncheck=0
do i=1,M
synd(i)=sum(cw(Nm(1:nrw(i),i)))
if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
! if( mod(synd(i),2) .ne. 0 ) write(*,*) 'check ',i,' unsatisfied'
enddo
! write(*,*) 'number of unsatisfied parity checks ',ncheck
if( ncheck .eq. 0 ) then ! we have a codeword - reorder the columns and return it
decoded=cw(1:K)
call chkcrc13a(decoded,nbadcrc)
if(nbadcrc.eq.0) then
message77=decoded(1:77)
nharderror=count( (2*cw-1)*llr .lt. 0.0 )
return
endif
endif
if( iter.gt.0 ) then ! this code block implements an early stopping criterion
! if( iter.gt.10000 ) then ! this code block implements an early stopping criterion
nd=ncheck-nclast
if( nd .lt. 0 ) then ! # of unsatisfied parity checks decreased
ncnt=0 ! reset counter
else
ncnt=ncnt+1
endif
! write(*,*) iter,ncheck,nd,ncnt
if( ncnt .ge. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) then
nharderror=-1
return
endif
endif
nclast=ncheck
! Send messages from bits to check nodes
do j=1,M
do i=1,nrw(j)
ibj=Nm(i,j)
toc(i,j)=zn(ibj)
do kk=1,4 ! subtract off what the bit had received from the check
if( Mn(kk,ibj) .eq. j ) then
toc(i,j)=toc(i,j)-tov(kk,ibj)
endif
enddo
enddo
enddo
! send messages from check nodes to variable nodes
do i=1,M
tanhtoc(1:12,i)=tanh(-toc(1:12,i)/2)
enddo
do j=1,N
do i=1,ncw(j)
ichk=Mn(i,j) ! Mn(:,j) are the checks that include bit j
Tmn=product(tanhtoc(1:nrw(ichk),ichk),mask=Nm(1:nrw(ichk),ichk).ne.j)
call platanh(-Tmn,y)
! y=atanh(-Tmn)
tov(i,j)=2*y
enddo
enddo
enddo
nharderror=-1
return
end subroutine bpdecode128_90

58
lib/encode_128_90.f90 Normal file
View File

@ -0,0 +1,58 @@
subroutine encode_128_90(message77,codeword)
!
! Add a 13-bit CRC to a 77-bit message and return a 128-bit codeword
!
use, intrinsic :: iso_c_binding
use iso_c_binding, only: c_loc,c_size_t
use crc
integer, parameter:: N=128, K=90, M=N-K
character*90 tmpchar
integer*1 codeword(N)
integer*1 gen(M,K)
integer*1 message77(77),message(K)
integer*1 pchecks(M)
integer*1, target :: i1MsgBytes(12)
include "ldpc_128_90_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=2
do jj=1, ibmax
icol=(j-1)*4+jj
if( btest(istr,4-jj) ) gen(i,icol)=1
enddo
enddo
enddo
first=.false.
endif
! Add 13 bit CRC to form 90-bit message+CRC13
write(tmpchar,'(77i1)') message77
tmpchar(78:80)='000'
i1MsgBytes=0
read(tmpchar,'(10b8)') i1MsgBytes(1:10)
ncrc13 = crc13 (c_loc (i1MsgBytes), 12)
write(tmpchar(78:90),'(b13)') ncrc13
read(tmpchar,'(90i1)') message
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 encode_128_90

114
lib/genmsk_128_90.f90 Normal file
View File

@ -0,0 +1,114 @@
subroutine genmsk_128_90(msg0,mygrid,ichk,bcontest,msgsent,i4tone,itype)
! s8 + 48bits + s8 + 80 bits = 144 bits (72ms message duration)
!
! Encode an MSK144 message
! Input:
! - msg0 requested message to be transmitted
! - ichk if ichk=1, return only msgsent
! if ichk.ge.10000, set imsg=ichk-10000 for short msg
! - msgsent message as it will be decoded
! - i4tone array of audio tone values, 0 or 1
! - itype message type
! 1 = standard message "Call_1 Call_2 Grid/Rpt"
! 2 = type 1 prefix
! 3 = type 1 suffix
! 4 = type 2 prefix
! 5 = type 2 suffix
! 6 = free text (up to 13 characters)
! 7 = short message "<Call_1 Call2> Rpt"
use iso_c_binding, only: c_loc,c_size_t
use packjt77
character*37 msg0
character*37 message !Message to be generated
character*37 msgsent !Message as it will be received
character*77 c77
character*6 mygrid
integer*4 i4tone(144)
integer*1 codeword(128)
integer*1 msgbits(77)
integer*1 bitseq(144) !Tone #s, data and sync (values 0-1)
integer*1 s8(8)
logical bcontest
real*8 pp(12)
real*8 xi(864),xq(864),pi,twopi
data s8/0,1,1,1,0,0,1,0/
equivalence (ihash,i1hash)
logical first
data first/.true./
save
if(first) then
first=.false.
nsym=128
pi=4.0*atan(1.0)
twopi=8.*atan(1.0)
do i=1,12
pp(i)=sin((i-1)*pi/12)
enddo
endif
if(msg0(1:1).eq.'@') then !Generate a fixed tone
read(msg0(2:5),*,end=1,err=1) nfreq !at specified frequency
go to 2
1 nfreq=1000
2 i4tone(1)=nfreq
else
message=msg0
do i=1,22
if(ichar(message(i:i)).eq.0) then
message(i:)=' '
exit
endif
enddo
do i=1,22 !Strip leading blanks
if(message(1:1).ne.' ') exit
message=message(i+1:)
enddo
if(message(1:1).eq.'<') then
call genmsk40(message,msgsent,ichk,i4tone,itype)
if(itype.lt.0) go to 999
i4tone(41)=-40
go to 999
endif
call pack77(message,i3,n3,c77)
call unpack77(c77,msgsent) !Unpack to get msgsent
if(ichk.eq.1) go to 999
read(c77,"77i1") msgbits
call encode_128_90(msgbits,codeword)
!Create 144-bit channel vector:
!8-bit sync word + 48 bits + 8-bit sync word + 80 bits
bitseq=0
bitseq(1:8)=s8
bitseq(9:56)=codeword(1:48)
bitseq(57:64)=s8
bitseq(65:144)=codeword(49:128)
bitseq=2*bitseq-1
xq(1:6)=bitseq(1)*pp(7:12) !first bit is mapped to 1st half-symbol on q
do i=1,71
is=(i-1)*12+7
xq(is:is+11)=bitseq(2*i+1)*pp
enddo
xq(864-5:864)=bitseq(1)*pp(1:6) !last half symbol
do i=1,72
is=(i-1)*12+1
xi(is:is+11)=bitseq(2*i)*pp
enddo
! Map I and Q to tones.
i4tone=0
do i=1,72
i4tone(2*i-1)=(bitseq(2*i)*bitseq(2*i-1)+1)/2;
i4tone(2*i)=-(bitseq(2*i)*bitseq(mod(2*i,144)+1)-1)/2;
enddo
endif
! Flip polarity
i4tone=-i4tone+1
999 return
end subroutine genmsk_128_90

View File

@ -1,132 +1,134 @@
program ldpcsim program ldpcsim128_90
use packjt ! Simulate the performance of the (128,90) code that is used in
integer, parameter:: NRECENT=10, N=128, K=90, M=N-K ! the second incarnation of MSK144.
character*12 recent_calls(NRECENT)
character*22 msg,msgsent,msgreceived
character*96 tmpchar
character*8 arg
integer*1 codeword(N), message77(77)
integer*1 apmask(N),cw(N)
integer*1 msgbits(77)
integer*4 i4Msg6BitWords(13)
integer nerrtot(0:N),nerrdec(0:N)
real*8 rxdata(N), rxavgd(N)
real llr(N)
do i=1,NRECENT use packjt77
recent_calls(i)=' ' integer, parameter:: NRECENT=10, N=128, K=90, M=N-K
enddo character*12 recent_calls(NRECENT)
nerrtot=0 character*37 msg,msgsent,msgreceived
nerrdec=0 character*77 c77
character*8 arg
integer*1 codeword(N), message77(77)
integer*1 apmask(N),cw(N)
integer*1 msgbits(77)
integer*4 i4Msg6BitWords(13)
integer nerrtot(0:N),nerrdec(0:N)
real*8 rxdata(N), rxavgd(N)
real llr(N)
nargs=iargc() do i=1,NRECENT
if(nargs.ne.4) then recent_calls(i)=' '
print*,'Usage: ldpcsim niter navg #trials s ' enddo
print*,'eg: ldpcsim 10 1 1000 0.75' nerrtot=0
return nerrdec=0
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) navg
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s
rate=real(K)/real(N) nargs=iargc()
if(nargs.ne.4) then
print*,'Usage: ldpcsim niter navg #trials s '
print*,'eg: ldpcsim 10 1 1000 0.75'
return
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) navg
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s
write(*,*) "rate: ",rate rate=real(77)/real(N)
write(*,*) "niter= ",max_iterations," navg= ",navg," s= ",s
write(*,*) "rate: ",rate
write(*,*) "niter= ",max_iterations," navg= ",navg," s= ",s
msg="K1ABC RR73; W9XYZ <KH1/KH7Z> -12"
i3=0
n3=1
call pack77(msg,i3,n3,c77)
call unpack77(c77,msgsent)
read(c77,'(77i1)') msgbits
!msg="K9AN K1JT EN50"
msg="G4WJS K1JT FN20"
call packmsg(msg,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes
call unpackmsg(i4Msg6BitWords,msgsent,.false.,' ') !Unpack to get msgsent
write(*,*) "message sent ",msgsent write(*,*) "message sent ",msgsent
tmpchar=' '
write(tmpchar,'(12b6.6)') i4Msg6BitWords(1:12)
tmpchar(73:77)="00000" !i5bit
read(tmpchar,'(77i1)') msgbits(1:77)
write(*,*) 'msgbits' write(*,*) 'msgbits'
write(*,'(28i1,1x,28i1,1x,16i1,1x,5i1)') msgbits write(*,'(77i1)') msgbits
! msgbits is the 77-bit message, codeword is 128 bits ! msgbits is the 77-bit message, codeword is 128 bits
call encode128_90(msgbits,codeword) call encode_128_90(msgbits,codeword)
call init_random_seed() ! call init_random_seed()
write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr" write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr"
do idb = 14,-6,-1 do idb = 14,0,-1
db=idb/2.0-1.0 db=idb/2.0-1.0
sigma=1/sqrt( 2*rate*(10**(db/10.0)) ) sigma=1/sqrt( 2*rate*(10**(db/10.0)) )
ngood=0 ngood=0
nue=0 nue=0
nbadcrc=0 nbadcrc=0
nsumerr=0 nsumerr=0
do itrial=1, ntrials do itrial=1, ntrials
rxavgd=0d0 rxavgd=0d0
do iav=1,navg do iav=1,navg
call sgran() ! call sgran()
! Create a realization of a noisy received word ! Create a realization of a noisy received word
do i=1,N do i=1,N
rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran() rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran()
enddo
rxavgd=rxavgd+rxdata
enddo enddo
rxavgd=rxavgd+rxdata rxdata=rxavgd
enddo nerr=0
rxdata=rxavgd do i=1,N
nerr=0 if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1
do i=1,N enddo
if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1 nerrtot(nerr)=nerrtot(nerr)+1
enddo
nerrtot(nerr)=nerrtot(nerr)+1
rxav=sum(rxdata)/N rxav=sum(rxdata)/N
rx2av=sum(rxdata*rxdata)/N rx2av=sum(rxdata*rxdata)/N
rxsig=sqrt(rx2av-rxav*rxav) rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig rxdata=rxdata/rxsig
! The s parameter can be tuned to trade a few tenth's dB of threshold for an order of ! The s parameter can be tuned to trade a few tenth's dB of threshold for an order of
! magnitude in UER ! magnitude in UER
if( s .lt. 0 ) then if( s .lt. 0 ) then
ss=sigma ss=sigma
else else
ss=s ss=s
endif endif
llr=2.0*rxdata/(ss*ss) llr=2.0*rxdata/(ss*ss)
apmask=0 apmask=0
! max_iterations is max number of belief propagation iterations ! max_iterations is max number of belief propagation iterations
call bpdecode128_90(llr, apmask, max_iterations, message77, cw, nharderrors, niterations) call bpdecode128_90(llr, apmask, max_iterations, message77, cw, nharderrors, niterations)
! If the decoder finds a valid codeword, nharderrors will be .ge. 0. ! If the decoder finds a valid codeword, nharderrors will be .ge. 0.
if( nharderrors .ge. 0 ) then if( nharderrors .ge. 0 ) then
call extractmessage77(message77,msgreceived) write(c77,'(77i1)') message77
nhw=count(cw.ne.codeword) call unpack77(c77,msgreceived)
if(nhw.eq.0) then ! this is a good decode nhw=count(cw.ne.codeword)
if(nhw.eq.0) then ! this is a good decode
ngood=ngood+1 ngood=ngood+1
nerrdec(nerr)=nerrdec(nerr)+1 nerrdec(nerr)=nerrdec(nerr)+1
else ! this is an undetected error else ! this is an undetected error
nue=nue+1 nue=nue+1
endif endif
endif endif
nsumerr=nsumerr+nerr nsumerr=nsumerr+nerr
enddo
snr2500=db+10*log10(rate*2000.0/2500.0) ! symbol rate is 2000 s^-1 and ref BW is 2500 Hz.
pberr=real(nsumerr)/real(ntrials*N)
write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,7x,f5.2,3x,e10.3)") db,snr2500,ngood,nue,ss,pberr
enddo enddo
snr2500=db-2.5 open(unit=23,file='nerrhisto.dat',status='unknown')
pberr=real(nsumerr)/real(ntrials*N) do i=0,N
write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,7x,f5.2,3x,e10.3)") db,snr2500,ngood,nue,ss,pberr write(23,'(i4,2x,i10,i10,f10.2)') i,nerrdec(i),nerrtot(i),real(nerrdec(i))/real(nerrtot(i)+1e-10)
enddo
enddo close(23)
open(unit=23,file='nerrhisto.dat',status='unknown') end program ldpcsim128_90
do i=0,N
write(23,'(i4,2x,i10,i10,f10.2)') i,nerrdec(i),nerrtot(i),real(nerrdec(i))/real(nerrtot(i)+1e-10)
enddo
close(23)
end program ldpcsim

View File

@ -41,7 +41,7 @@ call getarg(4,arg)
read(arg,*) s read(arg,*) s
! don't count hash bits as data bits ! don't count hash bits as data bits
K=80 K=72
N=128 N=128
rate=real(K)/real(N) rate=real(K)/real(N)
@ -161,7 +161,7 @@ do idb = 14,-6,-1
endif endif
endif endif
enddo enddo
snr2500=db-3.0 snr2500=db-3.5
write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2)") db,snr2500,ngood,nue,nbadhash,ss write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2)") db,snr2500,ngood,nue,nbadhash,ss
enddo enddo