From 2786c20ba2ddd5536dec4232b933807bb3a1dff5 Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Thu, 22 Mar 2018 14:20:07 +0000 Subject: [PATCH] Add support for a rate 1/3 (204,68) code. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8585 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/fsk4hf/bpdecode204.f90 | 482 ++++++++++++++++++++++++++++++ lib/fsk4hf/encode204.f90 | 48 +++ lib/fsk4hf/ldpc_204_68_params.f90 | 154 ++++++++++ lib/fsk4hf/ldpcsim204.f90 | 249 +++++++++++++++ lib/fsk4hf/osd204.f90 | 365 ++++++++++++++++++++++ 5 files changed, 1298 insertions(+) create mode 100644 lib/fsk4hf/bpdecode204.f90 create mode 100644 lib/fsk4hf/encode204.f90 create mode 100644 lib/fsk4hf/ldpc_204_68_params.f90 create mode 100644 lib/fsk4hf/ldpcsim204.f90 create mode 100644 lib/fsk4hf/osd204.f90 diff --git a/lib/fsk4hf/bpdecode204.f90 b/lib/fsk4hf/bpdecode204.f90 new file mode 100644 index 000000000..160b51a91 --- /dev/null +++ b/lib/fsk4hf/bpdecode204.f90 @@ -0,0 +1,482 @@ +subroutine bpdecode204(llr,apmask,maxiterations,decoded,cw,nharderror,iter) +! +! A log-domain belief propagation decoder for the (204,68) code. +! +integer, parameter:: N=204, K=68, M=N-K +integer*1 codeword(N),cw(N),apmask(N) +integer colorder(N) +integer*1 decoded(K) +integer Nm(6,M) ! 4, 5, or 6 bits per check +integer Mn(3,N) ! 3 checks per bit +integer synd(M) +real tov(3,N) +real toc(6,M) +real tanhtoc(6,M) +real zn(N) +real llr(N) +real Tmn +integer nrw(M) + +data colorder/ & + 0, 1, 2, 3, 4, 5, 47, 6, 7, 8, 9, 10, 11, 12, 58, 55, 13, & + 14, 15, 46, 17, 18, 60, 19, 20, 21, 22, 23, 24, 25, 57, 26, 27, 49, & + 28, 52, 65, 16, 50, 73, 59, 68, 63, 29, 30, 31, 32, 51, 62, 56, 66, & + 45, 33, 34, 53, 67, 35, 36, 37, 61, 69, 54, 38, 71, 82, 39, 77, 80, & + 83, 78, 84, 48, 41, 85, 40, 64, 75, 96, 74, 72, 76, 86, 87, 89, 90, & + 79, 70, 92, 99, 93,101, 95,100, 97, 94, 42, 98,103,105,102, 43,104, & + 88, 44,106, 81,107,110,108,111,112,109,113,114,117,118,116,121,115, & + 119,122,120,125,129,124,127,126,128, 91,123,133,131,130,134,135,137, & + 136,132,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152, & + 153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169, & + 170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186, & + 187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203/ + +data Mn/ & + 1, 38, 107, & + 2, 7, 114, & + 3, 48, 106, & + 4, 79, 94, & + 5, 97, 108, & + 6, 50, 122, & + 8, 78, 134, & + 9, 55, 65, & + 10, 62, 100, & + 11, 16, 99, & + 12, 113, 119, & + 13, 31, 125, & + 14, 15, 127, & + 17, 87, 103, & + 18, 81, 98, & + 19, 43, 77, & + 20, 102, 130, & + 21, 36, 111, & + 22, 23, 60, & + 24, 39, 112, & + 25, 37, 42, & + 26, 41, 51, & + 27, 67, 70, & + 28, 64, 136, & + 29, 61, 68, & + 30, 91, 124, & + 32, 80, 121, & + 33, 40, 117, & + 34, 35, 90, & + 44, 88, 93, & + 45, 128, 133, & + 46, 56, 69, & + 47, 49, 52, & + 53, 76, 131, & + 54, 104, 116, & + 57, 84, 86, & + 58, 120, 135, & + 59, 75, 92, & + 63, 71, 109, & + 66, 74, 126, & + 72, 85, 105, & + 73, 82, 95, & + 83, 89, 123, & + 96, 115, 118, & + 101, 110, 129, & + 52, 99, 132, & + 1, 3, 20, & + 2, 77, 89, & + 4, 72, 75, & + 5, 34, 79, & + 6, 24, 130, & + 7, 48, 88, & + 8, 36, 116, & + 9, 71, 114, & + 10, 87, 101, & + 11, 22, 121, & + 12, 50, 64, & + 13, 39, 53, & + 14, 41, 78, & + 15, 68, 96, & + 16, 83, 90, & + 17, 23, 45, & + 18, 47, 126, & + 19, 70, 91, & + 21, 57, 76, & + 25, 110, 117, & + 26, 82, 135, & + 27, 46, 58, & + 28, 37, 56, & + 29, 66, 102, & + 30, 62, 125, & + 31, 85, 93, & + 32, 104, 113, & + 33, 81, 92, & + 35, 100, 118, & + 38, 95, 133, & + 40, 86, 109, & + 42, 61, 124, & + 43, 59, 119, & + 44, 49, 134, & + 51, 97, 122, & + 54, 105, 107, & + 55, 128, 136, & + 60, 67, 84, & + 63, 112, 115, & + 65, 74, 131, & + 69, 80, 94, & + 73, 98, 123, & + 103, 130, 134, & + 46, 106, 111, & + 1, 84, 108, & + 120, 129, 132, & + 65, 75, 127, & + 2, 80, 101, & + 3, 118, 119, & + 4, 52, 124, & + 5, 13, 68, & + 6, 27, 81, & + 7, 51, 76, & + 8, 77, 108, & + 9, 31, 58, & + 10, 18, 57, & + 11, 63, 105, & + 12, 14, 132, & + 15, 56, 123, & + 16, 21, 128, & + 17, 37, 59, & + 19, 85, 126, & + 20, 71, 91, & + 22, 26, 117, & + 23, 79, 98, & + 24, 32, 95, & + 25, 90, 93, & + 28, 49, 109, & + 29, 116, 120, & + 30, 54, 136, & + 33, 53, 107, & + 34, 64, 103, & + 35, 39, 67, & + 36, 71, 73, & + 38, 47, 125, & + 40, 66, 94, & + 41, 70, 104, & + 42, 55, 112, & + 43, 44, 82, & + 29, 45, 88, & + 48, 86, 127, & + 50, 72, 135, & + 60, 74, 96, & + 61, 121, 131, & + 62, 78, 92, & + 69, 100, 133, & + 83, 122, 129, & + 87, 97, 106, & + 89, 102, 113, & + 24, 99, 108, & + 20, 72, 110, & + 111, 115, 117, & + 35, 52, 114, & + 1, 44, 94, & + 2, 23, 107, & + 3, 81, 136, & + 4, 8, 96, & + 5, 37, 70, & + 6, 43, 131, & + 7, 103, 115, & + 9, 94, 122, & + 10, 68, 82, & + 11, 56, 88, & + 12, 46, 126, & + 13, 16, 75, & + 14, 79, 112, & + 15, 47, 110, & + 17, 36, 39, & + 18, 63, 120, & + 19, 22, 55, & + 21, 49, 113, & + 25, 54, 57, & + 26, 89, 125, & + 27, 101, 109, & + 28, 31, 60, & + 30, 74, 97, & + 32, 92, 93, & + 33, 83, 91, & + 34, 58, 121, & + 38, 65, 111, & + 40, 99, 118, & + 3, 41, 61, & + 42, 50, 100, & + 45, 78, 106, & + 48, 95, 129, & + 51, 85, 133, & + 53, 59, 69, & + 11, 62, 66, & + 64, 73, 124, & + 67, 123, 134, & + 76, 104, 132, & + 77, 100, 127, & + 36, 80, 119, & + 84, 102, 135, & + 86, 105, 124, & + 4, 87, 128, & + 90, 106, 116, & + 65, 98, 130, & + 92, 108, 114, & + 1, 52, 121, & + 2, 84, 117, & + 5, 83, 105, & + 6, 15, 63, & + 7, 28, 82, & + 8, 32, 135, & + 9, 104, 134, & + 9, 10, 89, & + 12, 62, 107, & + 13, 40, 103, & + 14, 31, 95, & + 16, 27, 74, & + 17, 90, 132, & + 18, 34, 69, & + 19, 103, 129, & + 20, 76, 87, & + 21, 22, 130, & + 23, 25, 99, & + 24, 101, 126/ + +data Nm/ & + 1, 47, 91, 140, 186, 0, & + 2, 48, 94, 141, 187, 0, & + 3, 47, 95, 142, 168, 0, & + 4, 49, 96, 143, 182, 0, & + 5, 50, 97, 144, 188, 0, & + 6, 51, 98, 145, 189, 0, & + 2, 52, 99, 146, 190, 0, & + 7, 53, 100, 143, 191, 0, & + 8, 54, 101, 147, 192, 193, & + 9, 55, 102, 148, 193, 0, & + 10, 56, 103, 149, 174, 0, & + 11, 57, 104, 150, 194, 0, & + 12, 58, 97, 151, 195, 0, & + 13, 59, 104, 152, 196, 0, & + 13, 60, 105, 153, 189, 0, & + 10, 61, 106, 151, 197, 0, & + 14, 62, 107, 154, 198, 0, & + 15, 63, 102, 155, 199, 0, & + 16, 64, 108, 156, 200, 0, & + 17, 47, 109, 137, 201, 0, & + 18, 65, 106, 157, 202, 0, & + 19, 56, 110, 156, 202, 0, & + 19, 62, 111, 141, 203, 0, & + 20, 51, 112, 136, 204, 0, & + 21, 66, 113, 158, 203, 0, & + 22, 67, 110, 159, 0, 0, & + 23, 68, 98, 160, 197, 0, & + 24, 69, 114, 161, 190, 0, & + 25, 70, 115, 126, 0, 0, & + 26, 71, 116, 162, 0, 0, & + 12, 72, 101, 161, 196, 0, & + 27, 73, 112, 163, 191, 0, & + 28, 74, 117, 164, 0, 0, & + 29, 50, 118, 165, 199, 0, & + 29, 75, 119, 139, 0, 0, & + 18, 53, 120, 154, 179, 0, & + 21, 69, 107, 144, 0, 0, & + 1, 76, 121, 166, 0, 0, & + 20, 58, 119, 154, 0, 0, & + 28, 77, 122, 167, 195, 0, & + 22, 59, 123, 168, 0, 0, & + 21, 78, 124, 169, 0, 0, & + 16, 79, 125, 145, 0, 0, & + 30, 80, 125, 140, 0, 0, & + 31, 62, 126, 170, 0, 0, & + 32, 68, 90, 150, 0, 0, & + 33, 63, 121, 153, 0, 0, & + 3, 52, 127, 171, 0, 0, & + 33, 80, 114, 157, 0, 0, & + 6, 57, 128, 169, 0, 0, & + 22, 81, 99, 172, 0, 0, & + 33, 46, 96, 139, 186, 0, & + 34, 58, 117, 173, 0, 0, & + 35, 82, 116, 158, 0, 0, & + 8, 83, 124, 156, 0, 0, & + 32, 69, 105, 149, 0, 0, & + 36, 65, 102, 158, 0, 0, & + 37, 68, 101, 165, 0, 0, & + 38, 79, 107, 173, 0, 0, & + 19, 84, 129, 161, 0, 0, & + 25, 78, 130, 168, 0, 0, & + 9, 71, 131, 174, 194, 0, & + 39, 85, 103, 155, 189, 0, & + 24, 57, 118, 175, 0, 0, & + 8, 86, 93, 166, 184, 0, & + 40, 70, 122, 174, 0, 0, & + 23, 84, 119, 176, 0, 0, & + 25, 60, 97, 148, 0, 0, & + 32, 87, 132, 173, 199, 0, & + 23, 64, 123, 144, 0, 0, & + 39, 54, 109, 120, 0, 0, & + 41, 49, 128, 137, 0, 0, & + 42, 88, 120, 175, 0, 0, & + 40, 86, 129, 162, 197, 0, & + 38, 49, 93, 151, 0, 0, & + 34, 65, 99, 177, 201, 0, & + 16, 48, 100, 178, 0, 0, & + 7, 59, 131, 170, 0, 0, & + 4, 50, 111, 152, 0, 0, & + 27, 87, 94, 179, 0, 0, & + 15, 74, 98, 142, 0, 0, & + 42, 67, 125, 148, 190, 0, & + 43, 61, 133, 164, 188, 0, & + 36, 84, 91, 180, 187, 0, & + 41, 72, 108, 172, 0, 0, & + 36, 77, 127, 181, 0, 0, & + 14, 55, 134, 182, 201, 0, & + 30, 52, 126, 149, 0, 0, & + 43, 48, 135, 159, 193, 0, & + 29, 61, 113, 183, 198, 0, & + 26, 64, 109, 164, 0, 0, & + 38, 74, 131, 163, 185, 0, & + 30, 72, 113, 163, 0, 0, & + 4, 87, 122, 140, 147, 0, & + 42, 76, 112, 171, 196, 0, & + 44, 60, 129, 143, 0, 0, & + 5, 81, 134, 162, 0, 0, & + 15, 88, 111, 184, 0, 0, & + 10, 46, 136, 167, 203, 0, & + 9, 75, 132, 169, 178, 0, & + 45, 55, 94, 160, 204, 0, & + 17, 70, 135, 180, 0, 0, & + 14, 89, 118, 146, 195, 200, & + 35, 73, 123, 177, 192, 0, & + 41, 82, 103, 181, 188, 0, & + 3, 90, 134, 170, 183, 0, & + 1, 82, 117, 141, 194, 0, & + 5, 91, 100, 136, 185, 0, & + 39, 77, 114, 160, 0, 0, & + 45, 66, 137, 153, 0, 0, & + 18, 90, 138, 166, 0, 0, & + 20, 85, 124, 152, 0, 0, & + 11, 73, 135, 157, 0, 0, & + 2, 54, 139, 185, 0, 0, & + 44, 85, 138, 146, 0, 0, & + 35, 53, 115, 183, 0, 0, & + 28, 66, 110, 138, 187, 0, & + 44, 75, 95, 167, 0, 0, & + 11, 79, 95, 179, 0, 0, & + 37, 92, 115, 155, 0, 0, & + 27, 56, 130, 165, 186, 0, & + 6, 81, 133, 147, 0, 0, & + 43, 88, 105, 176, 0, 0, & + 26, 78, 96, 175, 181, 0, & + 12, 71, 121, 159, 0, 0, & + 40, 63, 108, 150, 204, 0, & + 13, 93, 127, 178, 0, 0, & + 31, 83, 106, 182, 0, 0, & + 45, 92, 133, 171, 200, 0, & + 17, 51, 89, 184, 202, 0, & + 34, 86, 130, 145, 0, 0, & + 46, 92, 104, 177, 198, 0, & + 31, 76, 132, 172, 0, 0, & + 7, 80, 89, 176, 192, 0, & + 37, 67, 128, 180, 191, 0, & + 24, 83, 116, 142, 0, 0/ + +data nrw/ & + 5,5,5,5,5,5,5,5,6,5,5,5,5,5,5,5,5, & + 5,5,5,5,5,5,5,5,4,5,5,4,4,5,5,4,5, & + 4,5,4,4,4,5,4,4,4,4,4,4,4,4,4,4,4, & + 5,4,4,4,4,4,4,4,4,4,5,5,4,5,4,4,4, & + 5,4,4,4,4,5,4,5,4,4,4,4,4,5,5,5,4, & + 4,5,4,5,5,4,5,4,5,5,4,4,4,5,5,5,4, & + 6,5,5,5,5,5,4,4,4,4,4,4,4,4,5,4,4, & + 4,5,4,4,5,4,5,4,4,5,5,4,5,4,5,5,4/ + +ncw=3 + +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(colorder+1) + decoded=codeword(M+1:N) + 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:6,i)=tanh(-toc(1:6,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 bpdecode204 diff --git a/lib/fsk4hf/encode204.f90 b/lib/fsk4hf/encode204.f90 new file mode 100644 index 000000000..b824196a7 --- /dev/null +++ b/lib/fsk4hf/encode204.f90 @@ -0,0 +1,48 @@ +subroutine encode204(message,codeword) +! Encode an 68-bit message and return a 204-bit codeword. +! The generator matrix has dimensions (136,68). +! The code is a (204,68) regular ldpc code with column weight 3. +! The code was generated using the PEG algorithm. +! After creating the codeword, the columns are re-ordered according to +! "colorder" to make the codeword compatible with the parity-check matrix +! + +include "ldpc_204_68_params.f90" + +integer*1 codeword(N) +integer*1 gen(M,K) +integer*1 itmp(N) +integer*1 message(K) +integer*1 pchecks(M) +logical first +data first/.true./ + +save first,gen + +if( first ) then ! fill the generator matrix + gen=0 + do i=1,M + do j=1,17 + read(g(i)(j:j),"(Z1)") istr + do jj=1, 4 + 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 +itmp(1:M)=pchecks +itmp(M+1:N)=message(1:K) +codeword(colorder+1)=itmp(1:N) + +return +end subroutine encode204 diff --git a/lib/fsk4hf/ldpc_204_68_params.f90 b/lib/fsk4hf/ldpc_204_68_params.f90 new file mode 100644 index 000000000..6e75ef114 --- /dev/null +++ b/lib/fsk4hf/ldpc_204_68_params.f90 @@ -0,0 +1,154 @@ +integer, parameter:: N=204, K=68, M=N-K +character*17 g(136) +integer colorder(N) +data g/ & !parity generator matrix for (204,68) code + "2de7435fd27c0031d", & + "f331b40671e20ea80", & + "48bd3f8cb9a24392f", & + "d4ed71c935162aa2a", & + "c437a3284ec58bce7", & + "35a806dd5be35627c", & + "396e797c33a4739a6", & + "768f331a59c15487b", & + "c214eac24ae5e1732", & + "0b5c53ff3a6da1192", & + "99624981d2703fb97", & + "e9f5447ef7f1ff6af", & + "bd8c730f0cfdf0727", & + "26f61e63e1e098f7f", & + "ef826566137b6526f", & + "af0e4fa251e9b4926", & + "75974a8b2a24292c5", & + "71caf0f2cd10f6d4f", & + "b1103f1f26e6898b7", & + "67ceb7d6f490da64f", & + "ee0e8fbefec23008a", & + "11cc2227e8bd676ca", & + "6e71626ba1e278046", & + "005d28da267e50e13", & + "a9ae4a130aaba8219", & + "d8ab72e0158d0da70", & + "56009d42b37bd66ff", & + "c39a75eca99b0e996", & + "6886de0bf7c0bf4bb", & + "1046cd8f64162f7b5", & + "da0f15843ac21e3a5", & + "e9bf9cd19f3db3913", & + "2fb9cb42d650f47a7", & + "a2b6c5a378fa75a65", & + "41a88f3cd60b79d6c", & + "fcf175794cc3ac96a", & + "8677a3447d40a9f71", & + "97a1f08c250b4bf12", & + "0168f090a1df6e8ea", & + "418a06bf372cc67d9", & + "0f17b880c1ff51239", & + "b2afd6d585deb961b", & + "60298ac5b58dbeee0", & + "8350c03c40119feff", & + "b29c964a8accf6af4", & + "9b46f036a5c178b5d", & + "917398bff051c300a", & + "5e52c03b2f8c5128c", & + "beae6c33c87ba38ab", & + "20843f7b056a02ebf", & + "66690d65acd9de598", & + "8f025841af5b54331", & + "b43cd869d3be2c3db", & + "c9c342fe63c18df50", & + "d331b40671e28ea80", & + "62406a0f4947e6ce9", & + "d67b1495883b22e1b", & + "734534c372408895b", & + "d88750e33d9677dcd", & + "6f96964da55138687", & + "80bee98bb75d50ef2", & + "c428ef3e3f06f4c56", & + "b1a1499b125883a35", & + "ac892d4b37fa9e395", & + "458dbda0f95ab11a5", & + "6f93c9e95b1094eed", & + "2e370d713914f848e", & + "758806dd5be35627c", & + "8c52e01caec798b49", & + "c286cc25bae3669cf", & + "87c56fb895c100884", & + "e89cb1376a18fd911", & + "156ffe5f30dc354e0", & + "f20d0b121d6a6b3ee", & + "7db08891b491a95d2", & + "191fac548d5077bdf", & + "023a37d7ea5660bbc", & + "6781668b363fee682", & + "bbfaf262cab7370da", & + "feea557965b7e474f", & + "c094eb223e1d305b8", & + "2be051abdd5beea35", & + "0790449880fda9d00", & + "f9029a39ec869e7b4", & + "5a29f48926ec9a552", & + "e0463306dc1470f87", & + "9251058334d790f86", & + "3019e1d4578e8a4dc", & + "887e46631502fa111", & + "c25fcd7a42465d326", & + "cf64bcc1056b555c4", & + "3e71c0fe5f0ad733b", & + "11055ec43b076e5b2", & + "3440f64dfa3c30a96", & + "2b73885b4d3299f60", & + "2e71627ba1e268046", & + "ad23743d5e6e5b80c", & + "c9757b05f29bfdc10", & + "f7112bea739247b51", & + "3664062387998b2b1", & + "90897a3b8785aefba", & + "29e126e3201fc1d46", & + "96c9001c84d5257fc", & + "067723447d40a9f71", & + "1a019cc68f7511402", & + "4bd48eb2330032763", & + "d139a5da936b37647", & + "765ab46a4dec5f04f", & + "706f475ad19b91955", & + "1755c988fa8a55e5c", & + "2fd9ed5777eb01d6a", & + "bec27d85b954d3fe8", & + "7135a3b92c45b3f8d", & + "353237872f002163a", & + "e31e4a97aef10c729", & + "da527d5e1cbc4edb6", & + "6e33cdede17c3207e", & + "ef2d2062e84dc401f", & + "8217c84c50c1bf833", & + "12ffbac7b2219c9e0", & + "3729178706f66881f", & + "2fdd748c382a608a1", & + "dd0a00076f9dcec73", & + "46b1d37bced447035", & + "7316f33a9c05ef178", & + "152c39a6de8954cc3", & + "16efffb7b62e12ba3", & + "9d9ec2bb467affd83", & + "467723445d40a9f61", & + "87994762b3bf50697", & + "b1bfa5b51526dde9b", & + "b0a6a19d709a96148", & + "990d567c0aba31a14", & + "171f190792461b1e0", & + "166011c27d2b6b8a4", & + "170c15831244ae73e"/ + +data colorder/ & + 0, 1, 2, 3, 4, 5, 47, 6, 7, 8, 9, 10, 11, 12, 58, 55, 13, & + 14, 15, 46, 17, 18, 60, 19, 20, 21, 22, 23, 24, 25, 57, 26, 27, 49, & + 28, 52, 65, 16, 50, 73, 59, 68, 63, 29, 30, 31, 32, 51, 62, 56, 66, & + 45, 33, 34, 53, 67, 35, 36, 37, 61, 69, 54, 38, 71, 82, 39, 77, 80, & + 83, 78, 84, 48, 41, 85, 40, 64, 75, 96, 74, 72, 76, 86, 87, 89, 90, & + 79, 70, 92, 99, 93,101, 95,100, 97, 94, 42, 98,103,105,102, 43,104, & + 88, 44,106, 81,107,110,108,111,112,109,113,114,117,118,116,121,115, & + 119,122,120,125,129,124,127,126,128, 91,123,133,131,130,134,135,137, & + 136,132,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152, & + 153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169, & + 170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186, & + 187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203/ diff --git a/lib/fsk4hf/ldpcsim204.f90 b/lib/fsk4hf/ldpcsim204.f90 new file mode 100644 index 000000000..224b2624f --- /dev/null +++ b/lib/fsk4hf/ldpcsim204.f90 @@ -0,0 +1,249 @@ +program ldpcsim204 + +! End-to-end test of the (300,60)/crc10 encoder and decoders. + +use crc +use packjt + +parameter(NRECENT=10) +character*12 recent_calls(NRECENT) +character*8 arg +integer*1, allocatable :: codeword(:), decoded(:), message(:) +integer*1, target:: i1Msg8BitBytes(9) +integer*1, target:: i1Dec8BitBytes(9) +integer*1 msgbits(68) +integer*1 apmask(204) +integer*1 cw(204) +integer*2 checksum +integer colorder(204) +integer nerrtot(204),nerrdec(204),nmpcbad(68) +logical checksumok,fsk,bpsk +real*8, allocatable :: rxdata(:) +real, allocatable :: llr(:) +real dllr(204),llrd(204) + +data colorder/ & + 0, 1, 2, 3, 4, 5, 47, 6, 7, 8, 9, 10, 11, 12, 58, 55, 13, & + 14, 15, 46, 17, 18, 60, 19, 20, 21, 22, 23, 24, 25, 57, 26, 27, 49, & + 28, 52, 65, 16, 50, 73, 59, 68, 63, 29, 30, 31, 32, 51, 62, 56, 66, & + 45, 33, 34, 53, 67, 35, 36, 37, 61, 69, 54, 38, 71, 82, 39, 77, 80, & + 83, 78, 84, 48, 41, 85, 40, 64, 75, 96, 74, 72, 76, 86, 87, 89, 90, & + 79, 70, 92, 99, 93,101, 95,100, 97, 94, 42, 98,103,105,102, 43,104, & + 88, 44,106, 81,107,110,108,111,112,109,113,114,117,118,116,121,115, & + 119,122,120,125,129,124,127,126,128, 91,123,133,131,130,134,135,137, & + 136,132,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152, & + 153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169, & + 170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186, & + 187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203/ + +do i=1,NRECENT + recent_calls(i)=' ' +enddo +nerrtot=0 +nerrdec=0 +nmpcbad=0 ! Used to collect the number of errors in the message+crc part of the codeword + +nargs=iargc() +if(nargs.ne.4) then + print*,'Usage: ldpcsim niter ndeep #trials s ' + print*,'eg: ldpcsim 100 4 1000 0.84' + 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,*) ndeep +call getarg(3,arg) +read(arg,*) ntrials +call getarg(4,arg) +read(arg,*) s + +fsk=.false. +bpsk=.true. + +N=204 +K=68 +rate=real(K)/real(N) + +write(*,*) "rate: ",rate +write(*,*) "niter= ",max_iterations," s= ",s + +allocate ( codeword(N), decoded(K), message(K) ) +allocate ( rxdata(N), llr(N) ) + +! The message should be packed into the first 7 bytes + i1Msg8BitBytes(1:6)=85 + i1Msg8BitBytes(7)=64 +! The CRC will be put into the last 2 bytes + i1Msg8BitBytes(8:9)=0 + checksum = crc10 (c_loc (i1Msg8BitBytes), 9) +! For reference, the next 3 lines show how to check the CRC + i1Msg8BitBytes(8)=checksum/256 + i1Msg8BitBytes(9)=iand (checksum,255) + checksumok = crc10_check(c_loc (i1Msg8BitBytes), 9) + if( checksumok ) write(*,*) 'Good checksum' +write(*,*) i1Msg8BitBytes(1:9) + + mbit=0 + do i=1, 7 + i1=i1Msg8BitBytes(i) + do ibit=1,8 + mbit=mbit+1 + msgbits(mbit)=iand(1,ishft(i1,ibit-8)) + enddo + enddo + i1=i1Msg8BitBytes(8) ! First 2 bits of crc10 are LSB of this byte + do ibit=1,2 + msgbits(50+ibit)=iand(1,ishft(i1,ibit-2)) + enddo + i1=i1Msg8BitBytes(9) ! Now shift in last 8 bits of the CRC + do ibit=1,8 + msgbits(52+ibit)=iand(1,ishft(i1,ibit-8)) + enddo + + write(*,*) 'message' + write(*,'(9(8i1,1x))') msgbits + + call encode204(msgbits,codeword) + call init_random_seed() + call sgran() + + write(*,*) 'codeword' + write(*,'(204i1)') codeword + +write(*,*) "Es/N0 SNR2500 ngood nundetected nbadcrc sigma" +do idb = 20,-18,-1 +!do idb = -16, -16, -1 + db=idb/2.0-1.0 +! sigma=1/sqrt( 2*rate*(10**(db/10.0)) ) ! to make db represent Eb/No + sigma=1/sqrt( 2*(10**(db/10.0)) ) ! db represents Es/No + ngood=0 + nue=0 + nbadcrc=0 + nberr=0 + do itrial=1, ntrials +! Create a realization of a noisy received word + do i=1,N + if( bpsk ) then + rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran() + elseif( fsk ) then + if( codeword(i) .eq. 1 ) then + r1=(1.0 + sigma*gran())**2 + (sigma*gran())**2 + r2=(sigma*gran())**2 + (sigma*gran())**2 + elseif( codeword(i) .eq. 0 ) then + r2=(1.0 + sigma*gran())**2 + (sigma*gran())**2 + r1=(sigma*gran())**2 + (sigma*gran())**2 + endif + rxdata(i)=0.35*(sqrt(r1)-sqrt(r2)) +! rxdata(i)=0.35*(exp(r1)-exp(r2)) +! rxdata(i)=0.12*(log(r1)-log(r2)) + endif + enddo + nerr=0 + do i=1,N + if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1 + enddo + if(nerr.ge.1) nerrtot(nerr)=nerrtot(nerr)+1 + nberr=nberr+nerr + +! 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) + apmask=0 +! max_iterations is max number of belief propagation iterations + call bpdecode204(llr,apmask,max_iterations,decoded,cw,nharderror,niterations) + if( (nharderror .lt. 0) .and. (ndeep .ge. 0) ) then + call osd204(llr, apmask, ndeep, decoded, cw, nhardmin, dmin) + niterations=nhardmin + endif + n2err=0 + do i=1,N + if( cw(i)*(2*codeword(i)-1.0) .lt. 0 ) n2err=n2err+1 + enddo +!write(*,*) nerr,niterations,n2err + damp=0.75 + ndither=0 + if( niterations .lt. 0 ) then + do i=1, ndither + do in=1,N + dllr(in)=damp*gran() + enddo + llrd=llr+dllr + call bpdecode300(llrd, apmask, max_iterations, decoded, cw, nharderror, niterations) + if( niterations .ge. 0 ) exit + enddo + endif + +! If the decoder finds a valid codeword, niterations will be .ge. 0. + if( niterations .ge. 0 ) then +! Check the CRC + do ibyte=1,6 + itmp=0 + do ibit=1,8 + itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit)) + enddo + i1Dec8BitBytes(ibyte)=itmp + enddo + i1Dec8BitBytes(7)=decoded(49)*128+decoded(50)*64 +! Need to pack the received crc into bytes 8 and 9 for crc10_check + i1Dec8BitBytes(8)=decoded(51)*2+decoded(52) + i1Dec8BitBytes(9)=decoded(53)*128+decoded(54)*64+decoded(55)*32+decoded(56)*16 + i1Dec8BitBytes(9)=i1Dec8BitBytes(9)+decoded(57)*8+decoded(58)*4+decoded(59)*2+decoded(60)*1 + ncrcflag=0 + if( crc10_check( c_loc( i1Dec8BitBytes ), 9 ) ) ncrcflag=1 + + if( ncrcflag .ne. 1 ) then + nbadcrc=nbadcrc+1 + endif + nueflag=0 + + nerrmpc=0 + do i=1,K ! find number of errors in message+crc part of codeword + if( msgbits(i) .ne. decoded(i) ) then + nueflag=1 + nerrmpc=nerrmpc+1 + endif + enddo + if(nerrmpc.ge.1) nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1 ! This histogram should inform our selection of CRC poly + if( ncrcflag .eq. 1 .and. nueflag .eq. 0 ) then + ngood=ngood+1 + if(nerr.ge.1) nerrdec(nerr)=nerrdec(nerr)+1 + else if( ncrcflag .eq. 1 .and. nueflag .eq. 1 ) then + nue=nue+1; + endif + endif + enddo + snr2500=db+10*log10(1.389/2500.0) + 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,120 + 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,68 + write(25,'(i4,2x,i10)') i,nmpcbad(i) +enddo +close(25) + + + +end program ldpcsim204 diff --git a/lib/fsk4hf/osd204.f90 b/lib/fsk4hf/osd204.f90 new file mode 100644 index 000000000..acebbf64d --- /dev/null +++ b/lib/fsk4hf/osd204.f90 @@ -0,0 +1,365 @@ +subroutine osd204(llr,apmask,ndeep,decoded,cw,nhardmin,dmin) +! +! An ordered-statistics decoder for the (204,68) code. +! +include "ldpc_204_68_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),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 decoded(K) +integer indx(N) +real llr(N),rx(N),absrx(N) +logical first,reset +data first/.true./ +save first,gen + +if( first ) then ! fill the generator matrix + gen=0 + do i=1,M + do j=1,17 + read(g(i)(j:j),"(Z1)") istr + do jj=1, 4 + irow=(j-1)*4+jj + if( btest(istr,4-jj) ) gen(irow,i)=1 + enddo + enddo + enddo + do irow=1,K + gen(irow,M+irow)=1 + enddo +first=.false. +endif + +! Re-order received vector to place systematic msg bits at the end. +rx=llr(colorder+1) +apmaskr=apmask(colorder+1) + +! 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(rx) +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 + +! 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). +do id=1,K ! diagonal element indices + do icol=id,K+20 ! The 20 is ad hoc - beware + iflag=0 + if( genmrb(id,icol) .eq. 1 ) then + iflag=1 + if( icol .ne. id ) then ! reorder column + temp(1:K)=genmrb(1:K,id) + genmrb(1:K,id)=genmrb(1:K,icol) + genmrb(1:K,icol)=temp(1:K) + itmp=indices(id) + indices(id)=indices(icol) + indices(icol)=itmp + endif + do ii=1,K + if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then + genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N)) + endif + enddo + exit + endif + enddo +enddo + +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=absrx(indices) +rx=rx(indices) +apmaskr=apmaskr(indices) + +call mrbencode(m0,c0,g2,N,K) +nxor=ieor(c0,hdec) +nhardmin=sum(nxor) +dmin=sum(nxor*absrx) + +cw=c0 +ntotal=0 +nrejected=0 + +if(ndeep.eq.0) goto 998 ! norder=0 +if(ndeep.gt.5) ndeep=5 +if( ndeep.eq. 1) then + nord=1 + npre1=0 + npre2=0 + nt=40 + ntheta=12 +elseif(ndeep.eq.2) then + nord=1 + 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 + 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 + misub(1:K-iorder)=0 + misub(K-iorder+1:K)=1 + iflag=K-iorder+1 + do while(iflag .ge.0) + if(iorder.eq.nord .and. npre1.eq.0) then + iend=iflag + else + iend=1 + endif + do n1=iflag,iend,-1 + mi=misub + mi(n1)=1 + if(any(iand(apmaskr(1:K),mi).eq.1)) cycle + ntotal=ntotal+1 + me=ieor(m0,mi) + if(n1.eq.iflag) then + call mrbencode(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + e2=e2sub + nd1Kpt=sum(e2sub(1:nt))+1 + d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K)) + else + e2=ieor(e2sub,g2(K+1:N,n1)) + nd1Kpt=sum(e2(1:nt))+2 + endif + if(nd1Kpt .le. ntheta) then + call mrbencode(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + if(n1.eq.iflag) then + dd=d1+sum(e2sub*absrx(K+1:N)) + else + dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N)) + endif + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + nd1Kptbest=nd1Kpt + endif + else + nrejected=nrejected+1 + endif + enddo +! Get the next test error pattern, iflag will go negative +! when the last pattern with weight iorder has been generated. + call nextpat(misub,k,iorder,iflag) + enddo +enddo + +if(npre2.eq.1) then + reset=.true. + ntotal=0 + do i1=K,1,-1 + do i2=i1-1,1,-1 + ntotal=ntotal+1 + mi(1:ntau)=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 + misub(1:K-nord)=0 + misub(K-nord+1:K)=1 + iflag=K-nord+1 + 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. +cw(indices)=cw +hdec(indices)=hdec +decoded=cw(M+1:N) +cw(colorder+1)=cw ! put the codeword back into received-word order +return +end subroutine osd204 + +subroutine mrbencode(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 mrbencode + +subroutine nextpat(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 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 + + 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 +