From 41920af89a50acdc8e46207932a3b070af2d06aa Mon Sep 17 00:00:00 2001 From: Steve Franke Date: Wed, 13 Jun 2018 09:40:04 -0500 Subject: [PATCH] Add routines that will implement 77-bit messages for MSK144. Initial tests will use a (128,90) code. --- lib/bpdecode128_90.f90 | 113 +++++++++++++++++ lib/crc13.cpp | 31 +++++ lib/encode128_90.f90 | 46 +++++++ lib/extractmessage128_90.f90 | 34 ++++++ lib/ldpc_128_90_generator.f90 | 41 +++++++ lib/ldpc_128_90_reordered_parity.f90 | 176 +++++++++++++++++++++++++++ lib/ldpcsim128_90.f90 | 169 +++++++++++++++++++++++++ 7 files changed, 610 insertions(+) create mode 100644 lib/bpdecode128_90.f90 create mode 100644 lib/crc13.cpp create mode 100644 lib/encode128_90.f90 create mode 100644 lib/extractmessage128_90.f90 create mode 100644 lib/ldpc_128_90_generator.f90 create mode 100644 lib/ldpc_128_90_reordered_parity.f90 create mode 100644 lib/ldpcsim128_90.f90 diff --git a/lib/bpdecode128_90.f90 b/lib/bpdecode128_90.f90 new file mode 100644 index 000000000..89bff000e --- /dev/null +++ b/lib/bpdecode128_90.f90 @@ -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 diff --git a/lib/crc13.cpp b/lib/crc13.cpp new file mode 100644 index 000000000..03265088a --- /dev/null +++ b/lib/crc13.cpp @@ -0,0 +1,31 @@ +#include +#include + +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); +} diff --git a/lib/encode128_90.f90 b/lib/encode128_90.f90 new file mode 100644 index 000000000..48a13fbfb --- /dev/null +++ b/lib/encode128_90.f90 @@ -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 diff --git a/lib/extractmessage128_90.f90 b/lib/extractmessage128_90.f90 new file mode 100644 index 000000000..1282ceea1 --- /dev/null +++ b/lib/extractmessage128_90.f90 @@ -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 diff --git a/lib/ldpc_128_90_generator.f90 b/lib/ldpc_128_90_generator.f90 new file mode 100644 index 000000000..200d17036 --- /dev/null +++ b/lib/ldpc_128_90_generator.f90 @@ -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"/ diff --git a/lib/ldpc_128_90_reordered_parity.f90 b/lib/ldpc_128_90_reordered_parity.f90 new file mode 100644 index 000000000..b3a5e9b85 --- /dev/null +++ b/lib/ldpc_128_90_reordered_parity.f90 @@ -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 diff --git a/lib/ldpcsim128_90.f90 b/lib/ldpcsim128_90.f90 new file mode 100644 index 000000000..1bb0b4e59 --- /dev/null +++ b/lib/ldpcsim128_90.f90 @@ -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