Add routines that will implement 77-bit messages for MSK144. Initial tests will use a (128,90) code.

This commit is contained in:
Steve Franke 2018-06-13 09:40:04 -05:00
parent a5837bd5c7
commit 41920af89a
7 changed files with 610 additions and 0 deletions

113
lib/bpdecode128_90.f90 Normal file
View File

@ -0,0 +1,113 @@
subroutine bpdecode128_90(llr,apmask,maxiterations,decoded,cw,nharderror,iter)
!
! A log-domain belief propagation decoder for the (128,90) code.
!
integer, parameter:: N=128, K=90, M=N-K
integer*1 codeword(N),cw(N),apmask(N)
integer*1 decoded(K)
integer Nm(11,M)
integer Mn(3,N) ! 3 checks per bit
integer synd(M)
real tov(3,N)
real toc(11,M)
real tanhtoc(11,M)
real zn(N)
real llr(N)
real Tmn
integer nrw(M),ncw
include "ldpc_128_90_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))
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
codeword=cw
decoded=codeword(1:K)
nerr=0
do i=1,N
if( (2*cw(i)-1)*llr(i) .lt. 0.0 ) nerr=nerr+1
enddo
nharderror=nerr
return
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. 5 .and. iter .ge. 10 .and. ncheck .gt. 15) 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,ncw ! 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:7,i)=tanh(-toc(1:7,i)/2)
enddo
do j=1,N
do i=1,ncw
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

31
lib/crc13.cpp Normal file
View File

@ -0,0 +1,31 @@
#include <boost/crc.hpp>
#include <boost/config.hpp>
extern "C"
{
short crc13 (unsigned char const * data, int length);
bool crc13_check (unsigned char const * data, int length);
}
#define POLY 0x15D7
#ifdef BOOST_NO_CXX11_CONSTEXPR
#define TRUNCATED_POLYNOMIAL POLY
#else
namespace
{
unsigned long constexpr TRUNCATED_POLYNOMIAL = POLY;
}
#endif
// assumes CRC is last 13 bits of the data and is set to zero
// caller should assign the returned CRC into the message in big endian byte order
short crc13 (unsigned char const * data, int length)
{
return boost::augmented_crc<13, TRUNCATED_POLYNOMIAL> (data, length);
}
bool crc13_check (unsigned char const * data, int length)
{
return !boost::augmented_crc<13, TRUNCATED_POLYNOMIAL> (data, length);
}

46
lib/encode128_90.f90 Normal file
View File

@ -0,0 +1,46 @@
subroutine encode128_90(message,codeword)
! Encode an 90-bit message and return a 128-bit codeword.
! The generator matrix has dimensions (38,90).
! The code is a (128,90) regular ldpc code with column weight 3.
!
integer, parameter:: N=128, K=90, M=N-K
integer*1 codeword(N)
integer*1 gen(M,K)
integer*1 message(K)
integer*1 pchecks(M)
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
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 encode128_90

View File

@ -0,0 +1,34 @@
subroutine extractmessage128_90(decoded,msgreceived,ncrcflag)
use iso_c_binding, only: c_loc,c_size_t
use crc
use packjt
character*22 msgreceived
character*90 cbits
integer*1 decoded(90)
integer*1, target:: i1Dec8BitBytes(12)
integer*4 i4Dec6BitWords(12)
! Write decoded bits into cbits: 77-bit message plus 13-bit CRC
write(cbits,1000) decoded
1000 format(90i1)
read(cbits,1001) i1Dec8BitBytes
1001 format(12b8)
read(cbits,1002) ncrc13 !Received CRC12
1002 format(77x,b13)
i1Dec8BitBytes(10)=iand(i1Dec8BitBytes(10),128+64+32+16+8)
i1Dec8BitBytes(11:12)=0
icrc13=crc13(c_loc(i1Dec8BitBytes),12) !CRC13 computed from 77 msg bits
if(ncrc13.eq.icrc13 .or. sum(decoded(57:87)).eq.0) then !### Kludge ### ???
! CRC13 checks out --- unpack 72-bit message
read(cbits,'(12b6)') i4Dec6BitWords
call unpackmsg(i4Dec6BitWords,msgreceived,.false.,' ')
ncrcflag=1
else
msgreceived=' '
ncrcflag=-1
endif
return
end subroutine extractmessage128_90

View File

@ -0,0 +1,41 @@
character*23 g(38)
data g/ &
"a08ea80879050a5e94da994", &
"59f3b48040ca089c81ee880", &
"e4070262802e31b7b17d3dc", &
"95cbcbaf032dc3d960bacc8", &
"c4d79b5dcc21161a254ffbc", &
"93fde9cdbf2622a70868424", &
"e73b888bb1b01167379ba28", &
"45a0d0a0f39a7ad2439949c", &
"759acef19444bcad79c4964", &
"71eb4dddf4f5ed9e2ea17e0", &
"80f0ad76fb247d6b4ca8d38", &
"184fff3aa1b82dc66640104", &
"ca4e320bb382ed14cbb1094", &
"52514447b90e25b9e459e28", &
"dd10c1666e071956bd0df38", &
"99c332a0b792a2da8ef1ba8", &
"7bd9f688e7ed402e231aaac", &
"00fcad76eb647d6a0ca8c38", &
"6ac8d0499c43b02eed78d70", &
"2c2c764baf795b4788db010", &
"0e907bf9e280d2624823dd0", &
"b857a6e315afd8c1c925e64", &
"8deb58e22d73a141cae3778", &
"22d3cb80d92d6ac132dfe08", &
"754763877b28c187746855c", &
"1d1bb7cf6953732e04ebca4", &
"2c65e0ea4466ab9f5e1deec", &
"6dc530ca37fc916d1f84870", &
"49bccbbee152355be7ac984", &
"e8387f3f4367cf45a150448", &
"8ce25e03d67d51091c81884", &
"b798012ffa40a93852752c8", &
"2e43307933adfca37adc3c8", &
"ca06e0a42ca1ec782d6c06c", &
"c02b762927556a7039e638c", &
"4a3e9b7d08b6807f8619fac", &
"45e8030f68997bb68544424", &
"7e79362c16773efc6482e30"/

View File

@ -0,0 +1,176 @@
data Mn/ &
21, 34, 36, &
1, 8, 28, &
2, 9, 37, &
3, 7, 19, &
4, 16, 32, &
2, 5, 22, &
6, 13, 25, &
10, 31, 33, &
11, 24, 27, &
12, 15, 23, &
14, 18, 26, &
17, 20, 29, &
17, 30, 34, &
6, 34, 35, &
1, 10, 30, &
3, 18, 23, &
4, 12, 25, &
5, 28, 36, &
7, 14, 21, &
8, 15, 31, &
9, 27, 32, &
11, 19, 35, &
13, 16, 37, &
20, 24, 38, &
21, 22, 26, &
12, 29, 33, &
1, 17, 35, &
2, 28, 30, &
3, 10, 32, &
4, 8, 36, &
5, 19, 29, &
6, 20, 27, &
7, 22, 37, &
9, 11, 33, &
13, 24, 26, &
14, 31, 34, &
15, 16, 25, &
13, 18, 38, &
8, 20, 23, &
1, 32, 33, &
2, 17, 19, &
3, 24, 34, &
4, 7, 38, &
5, 11, 31, &
6, 18, 21, &
9, 15, 36, &
10, 16, 28, &
12, 26, 30, &
14, 27, 29, &
22, 25, 35, &
23, 30, 32, &
4, 11, 37, &
1, 14, 23, &
2, 8, 25, &
3, 13, 27, &
5, 10, 37, &
6, 16, 31, &
7, 15, 18, &
9, 22, 24, &
12, 19, 36, &
17, 26, 38, &
20, 21, 33, &
20, 28, 35, &
4, 29, 34, &
1, 26, 36, &
2, 23, 34, &
3, 9, 38, &
5, 6, 17, &
7, 27, 35, &
8, 14, 32, &
10, 15, 22, &
11, 18, 29, &
12, 13, 28, &
16, 19, 33, &
21, 25, 31, &
24, 30, 37, &
1, 3, 21, &
2, 18, 31, &
4, 6, 9, &
5, 8, 33, &
7, 29, 32, &
10, 13, 19, &
11, 22, 23, &
12, 27, 34, &
14, 15, 30, &
16, 27, 38, &
17, 28, 37, &
20, 25, 26, &
5, 24, 35, &
3, 6, 36, &
1, 12, 31, &
2, 4, 33, &
3, 16, 30, &
1, 2, 24, &
5, 23, 27, &
6, 28, 32, &
7, 17, 36, &
8, 22, 38, &
9, 18, 20, &
10, 21, 29, &
11, 13, 34, &
4, 14, 20, &
11, 30, 38, &
14, 35, 37, &
15, 19, 26, &
3, 28, 29, &
7, 8, 9, &
5, 18, 34, &
13, 15, 17, &
12, 16, 35, &
10, 23, 25, &
19, 21, 37, &
17, 27, 31, &
24, 25, 36, &
1, 18, 19, &
6, 26, 33, &
22, 31, 32, &
3, 20, 22, &
4, 21, 27, &
2, 13, 29, &
6, 7, 12, &
15, 24, 32, &
9, 25, 30, &
23, 37, 38, &
5, 16, 26, &
11, 14, 28, &
33, 36, 38, &
8, 10, 35/
data Nm/ &
2, 15, 27, 40, 53, 65, 77, 91, 94, 115, 0, &
3, 6, 28, 41, 54, 66, 78, 92, 94, 120, 0, &
4, 16, 29, 42, 55, 67, 77, 90, 93, 106, 118, &
5, 17, 30, 43, 52, 64, 79, 92, 102, 119, 0, &
6, 18, 31, 44, 56, 68, 80, 89, 95, 108, 125, &
7, 14, 32, 45, 57, 68, 79, 90, 96, 116, 121, &
4, 19, 33, 43, 58, 69, 81, 97, 107, 121, 0, &
2, 20, 30, 39, 54, 70, 80, 98, 107, 128, 0, &
3, 21, 34, 46, 59, 67, 79, 99, 107, 123, 0, &
8, 15, 29, 47, 56, 71, 82, 100, 111, 128, 0, &
9, 22, 34, 44, 52, 72, 83, 101, 103, 126, 0, &
10, 17, 26, 48, 60, 73, 84, 91, 110, 121, 0, &
7, 23, 35, 38, 55, 73, 82, 101, 109, 120, 0, &
11, 19, 36, 49, 53, 70, 85, 102, 104, 126, 0, &
10, 20, 37, 46, 58, 71, 85, 105, 109, 122, 0, &
5, 23, 37, 47, 57, 74, 86, 93, 110, 125, 0, &
12, 13, 27, 41, 61, 68, 87, 97, 109, 113, 0, &
11, 16, 38, 45, 58, 72, 78, 99, 108, 115, 0, &
4, 22, 31, 41, 60, 74, 82, 105, 112, 115, 0, &
12, 24, 32, 39, 62, 63, 88, 99, 102, 118, 0, &
1, 19, 25, 45, 62, 75, 77, 100, 112, 119, 0, &
6, 25, 33, 50, 59, 71, 83, 98, 117, 118, 0, &
10, 16, 39, 51, 53, 66, 83, 95, 111, 124, 0, &
9, 24, 35, 42, 59, 76, 89, 94, 114, 122, 0, &
7, 17, 37, 50, 54, 75, 88, 111, 114, 123, 0, &
11, 25, 35, 48, 61, 65, 88, 105, 116, 125, 0, &
9, 21, 32, 49, 55, 69, 84, 86, 95, 113, 119, &
2, 18, 28, 47, 63, 73, 87, 96, 106, 126, 0, &
12, 26, 31, 49, 64, 72, 81, 100, 106, 120, 0, &
13, 15, 28, 48, 51, 76, 85, 93, 103, 123, 0, &
8, 20, 36, 44, 57, 75, 78, 91, 113, 117, 0, &
5, 21, 29, 40, 51, 70, 81, 96, 117, 122, 0, &
8, 26, 34, 40, 62, 74, 80, 92, 116, 127, 0, &
1, 13, 14, 36, 42, 64, 66, 84, 101, 108, 0, &
14, 22, 27, 50, 63, 69, 89, 104, 110, 128, 0, &
1, 18, 30, 46, 60, 65, 90, 97, 114, 127, 0, &
3, 23, 33, 52, 56, 76, 87, 104, 112, 124, 0, &
24, 38, 43, 61, 67, 86, 98, 103, 124, 127, 0/
data nrw/ &
10,10,11,10,11,11,10,10,10,10,10,10,10,10,10,10,10,10, &
10,10,10,10,10,10,10,10,11,10,10,10,10,10,10,10,10,10, &
10,10/
ncw=3

169
lib/ldpcsim128_90.f90 Normal file
View File

@ -0,0 +1,169 @@
program ldpcsim
use, intrinsic :: iso_c_binding
use iso_c_binding, only: c_loc,c_size_t
use crc
use packjt
integer, parameter:: NRECENT=10, N=128, K=90, M=N-K
character*12 recent_calls(NRECENT)
character*22 msg,msgsent,msgreceived
character*96 tmpchar
character*8 arg
integer*1, allocatable :: codeword(:), decoded(:), message(:)
integer*1, target:: i1Msg8BitBytes(12)
integer*1 i1hash(4)
integer*1 msgbits(80)
integer*4 i4Msg6BitWords(13)
integer ihash
integer nerrtot(N),nerrdec(N),nmpcbad(K)
real*8, allocatable :: lratio(:), rxdata(:), rxavgd(:)
real, allocatable :: yy(:), llr(:)
equivalence(ihash,i1hash)
do i=1,NRECENT
recent_calls(i)=' '
enddo
nerrtot=0
nerrdec=0
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
! don't count hash bits as data bits
rate=real(K)/real(N)
write(*,*) "rate: ",rate
write(*,*) "niter= ",max_iterations," navg= ",navg," s= ",s
allocate ( codeword(N), decoded(K), message(K) )
allocate ( lratio(N), rxdata(N), rxavgd(N), yy(N), llr(N) )
msg="K9AN K1JT EN50"
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
tmpchar=' '
write(tmpchar,'(12b6.6)') i4Msg6BitWords
tmpchar(73:77)="00000" !i5bit
write(*,*) tmpchar
read(tmpchar,'(10b8)') i1Msg8BitBytes(1:10)
write(*,*) i1Msg8BitBytes
i1Msg8BitBytes(10:12)=0
checksum = crc13 (c_loc (i1Msg8ZBitZBytes), 12)
i1Msg8BitBytes(11)=checksum/256
i1Msg8BitBytes(12)=iand (checksum,255)
checksumok = crc13_check(c_loc (i1Msg8ZBitBytes), 12)
if( checksumok ) write(*,*) 'Good checksum'
write(tmpchar,'(12b8.8)') i1Msg8BitBytes(1:9)
read(tmpchar,'(77b)') msgbits(1:77)
read(tmpchar(84:96),'(6b)') msgbits(78:90)
call encode_128_90(msgbits,codeword)
call init_random_seed()
write(*,*) "Eb/N0 SNR2500 ngood nundetected nbadhash sigma"
do idb = -6, 14
db=idb/2.0-1.0
sigma=1/sqrt( 2*rate*(10**(db/10.0)) )
ngood=0
nue=0
nbadhash=0
do itrial=1, ntrials
rxavgd=0d0
do iav=1,navg
call sgran()
! Create a realization of a noisy received word
do i=1,N
rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran()
enddo
rxavgd=rxavgd+rxdata
enddo
rxdata=rxavgd
nerr=0
do i=1,N
if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1
enddo
nerrtot(nerr)=nerrtot(nerr)+1
! Correct signal normalization is important for this decoder.
rxav=sum(rxdata)/N
rx2av=sum(rxdata*rxdata)/N
rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig
! To match the metric to the channel, s should be set to the noise standard deviation.
! For now, set s to the value that optimizes decode probability near threshold.
! The s parameter can be tuned to trade a few tenth's dB of threshold for an order of
! magnitude in UER
if( s .lt. 0 ) then
ss=sigma
else
ss=s
endif
llr=2.0*rxdata/(ss*ss)
lratio=exp(llr)
yy=rxdata
! max_iterations is max number of belief propagation iterations
call bpdecode128_90(llr, apmask, max_iterations, decoded, cw, nharderrors, niterations)
! If the decoder finds a valid codeword, nharderrors will be .ge. 0.
if( nharderrors .ge. 0 ) then
call extractmessage1128_90(decoded,msgreceived,ncrcflag)
if( nncrcflag .ne. 1 ) then
nbadcrc=nbadcrc+1
endif
nueflag=0
nerrmpc=0
do i=1,K
if( msgbits(i) .ne. decoded(i) ) then
nueflag=1
nerrmpc=nerrmpc+1
endif
enddo
if(nerrmpc.ge.1) nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1
if( ncrcflag .eq. 1) then
ngood=ngood+1
if(nerr.ge.1) nerrdec(nerr)=nerrdec(nerr)+1
else if(nueflag .eq. 1 ) then
nue=nue+1;
endif
endif
enddo
snr2500=db-3.5
pberr=real(nberr)/real(ntrials*N)
write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2,8x,e10.3)") db,snr2500,ngood,nue,nbadcrc,ss,pberr
enddo
open(unit=23,file='nerrhisto.dat',status='unknown')
do i=1,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)
open(unit=25,file='nmpcbad.dat',status='unknown')
do i=1,K
write(25,'(i4,2x,i10)') i,nmpcbad(i)
enddo
close(25)
end program ldpcsim