diff --git a/lib/fsk4hf/bpdecode174_91.f90 b/lib/fsk4hf/bpdecode174_91.f90 new file mode 100644 index 000000000..740e3175a --- /dev/null +++ b/lib/fsk4hf/bpdecode174_91.f90 @@ -0,0 +1,393 @@ +subroutine bpdecode174_91(llr,apmask,maxiterations,decoded,cw,nharderror,iter) +! +! A log-domain belief propagation decoder for the (174,91) code. +! +integer, parameter:: N=174, K=91, M=N-K +integer*1 codeword(N),cw(N),apmask(N) +integer colorder(N) +integer*1 decoded(K) +integer Nm(7,M) +integer Mn(3,N) ! 3 checks per bit +integer synd(M) +real tov(3,N) +real toc(7,M) +real tanhtoc(7,M) +real zn(N) +real llr(N) +real Tmn +integer nrw(M) + +data colorder/ & + 0, 1, 2, 3, 28, 4, 5, 6, 7, 8, 9, 10, 11, 34, 12, 32, 13, 14, 15, 16,& + 17, 18, 36, 29, 40, 19, 20, 38, 21, 41, 30, 42, 22, 44, 37, 47, 48, 23, 33, 43,& + 49, 45, 56, 39, 25, 26, 46, 50, 51, 52, 24, 57, 58, 61, 31, 54, 64, 35, 27, 62,& + 59, 53, 60, 63, 55, 70, 66, 67, 68, 65, 71, 74, 72, 73, 77, 75, 69, 76, 79, 82,& + 83, 78, 81, 80, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,& + 100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,& + 120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,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/ + +data Mn/ & + 1, 24, 66, & + 2, 5, 70, & + 3, 31, 65, & + 4, 49, 58, & + 6, 60, 67, & + 7, 32, 75, & + 8, 48, 82, & + 9, 35, 41, & + 10, 39, 62, & + 11, 14, 61, & + 12, 71, 74, & + 13, 23, 78, & + 15, 16, 79, & + 17, 54, 63, & + 18, 50, 57, & + 19, 30, 47, & + 20, 64, 80, & + 21, 28, 69, & + 22, 25, 43, & + 26, 34, 72, & + 27, 36, 37, & + 29, 40, 44, & + 33, 52, 53, & + 38, 55, 83, & + 42, 51, 59, & + 45, 76, 81, & + 46, 68, 77, & + 56, 67, 73, & + 1, 4, 5, & + 2, 47, 51, & + 3, 46, 82, & + 6, 24, 76, & + 7, 9, 16, & + 8, 10, 78, & + 11, 35, 55, & + 12, 38, 64, & + 13, 42, 83, & + 14, 27, 54, & + 15, 21, 34, & + 17, 44, 53, & + 18, 25, 28, & + 19, 33, 57, & + 20, 22, 73, & + 23, 40, 81, & + 26, 49, 68, & + 29, 71, 75, & + 30, 65, 79, & + 31, 36, 60, & + 32, 43, 77, & + 37, 62, 70, & + 39, 69, 74, & + 41, 52, 66, & + 45, 50, 61, & + 48, 63, 80, & + 56, 59, 72, & + 58, 64, 65, & + 1, 13, 28, & + 2, 48, 75, & + 3, 53, 69, & + 4, 11, 44, & + 5, 73, 79, & + 6, 12, 17, & + 7, 57, 60, & + 8, 15, 61, & + 9, 39, 59, & + 10, 19, 49, & + 14, 43, 52, & + 16, 54, 68, & + 18, 41, 63, & + 20, 36, 45, & + 21, 67, 77, & + 10, 22, 55, & + 23, 65, 72, & + 24, 27, 82, & + 25, 26, 29, & + 30, 35, 37, & + 31, 51, 66, & + 17, 32, 78, & + 33, 42, 76, & + 34, 70, 83, & + 38, 46, 81, & + 40, 62, 80, & + 45, 47, 74, & + 50, 56, 71, & + 7, 37, 58, & + 1, 16, 71, & + 2, 6, 61, & + 3, 22, 50, & + 4, 59, 77, & + 5, 41, 81, & + 8, 58, 74, & + 9, 20, 26, & + 11, 21, 31, & + 12, 66, 79, & + 13, 14, 57, & + 15, 33, 40, & + 18, 44, 82, & + 19, 69, 83, & + 23, 49, 63, & + 24, 29, 39, & + 25, 47, 56, & + 27, 55, 72, & + 28, 64, 70, & + 30, 48, 77, & + 32, 34, 45, & + 35, 68, 80, & + 36, 38, 52, & + 42, 43, 62, & + 46, 60, 78, & + 51, 54, 67, & + 53, 73, 75, & + 14, 73, 76, & + 1, 22, 30, & + 2, 35, 43, & + 3, 47, 63, & + 4, 25, 76, & + 5, 33, 78, & + 6, 20, 83, & + 7, 12, 72, & + 8, 54, 70, & + 9, 61, 65, & + 10, 34, 51, & + 11, 46, 75, & + 13, 39, 68, & + 15, 17, 56, & + 16, 23, 36, & + 18, 32, 55, & + 19, 31, 81, & + 21, 37, 71, & + 24, 57, 64, & + 26, 38, 48, & + 27, 49, 50, & + 28, 52, 59, & + 29, 41, 58, & + 40, 60, 74, & + 42, 44, 79, & + 51, 53, 80, & + 62, 67, 82, & + 23, 66, 69, & + 1, 53, 61, & + 2, 18, 39, & + 3, 4, 12, & + 5, 26, 74, & + 6, 30, 52, & + 7, 82, 83, & + 8, 35, 73, & + 9, 19, 67, & + 10, 64, 75, & + 11, 20, 33, & + 13, 45, 48, & + 3, 14, 40, & + 15, 43, 49, & + 16, 55, 76, & + 17, 62, 65, & + 21, 47, 78, & + 22, 59, 81, & + 24, 34, 63, & + 25, 37, 66, & + 27, 79, 80, & + 28, 60, 79, & + 29, 31, 70, & + 32, 58, 69, & + 10, 36, 77, & + 38, 50, 51, & + 13, 41, 56, & + 42, 63, 71, & + 44, 47, 68, & + 1, 46, 72, & + 54, 57, 75, & + 2, 33, 58, & + 4, 17, 83, & + 5, 14, 55, & + 6, 23, 48, & + 7, 52, 56/ + +data Nm/ & + 1, 29, 57, 86, 113, 140, 168, & + 2, 30, 58, 87, 114, 141, 170, & + 3, 31, 59, 88, 115, 142, 151, & + 4, 29, 60, 89, 116, 142, 171, & + 2, 29, 61, 90, 117, 143, 172, & + 5, 32, 62, 87, 118, 144, 173, & + 6, 33, 63, 85, 119, 145, 174, & + 7, 34, 64, 91, 120, 146, 0, & + 8, 33, 65, 92, 121, 147, 0, & + 9, 34, 66, 72, 122, 148, 163, & + 10, 35, 60, 93, 123, 149, 0, & + 11, 36, 62, 94, 119, 142, 0, & + 12, 37, 57, 95, 124, 150, 165, & + 10, 38, 67, 95, 112, 151, 172, & + 13, 39, 64, 96, 125, 152, 0, & + 13, 33, 68, 86, 126, 153, 0, & + 14, 40, 62, 78, 125, 154, 171, & + 15, 41, 69, 97, 127, 141, 0, & + 16, 42, 66, 98, 128, 147, 0, & + 17, 43, 70, 92, 118, 149, 0, & + 18, 39, 71, 93, 129, 155, 0, & + 19, 43, 72, 88, 113, 156, 0, & + 12, 44, 73, 99, 126, 139, 173, & + 1, 32, 74, 100, 130, 157, 0, & + 19, 41, 75, 101, 116, 158, 0, & + 20, 45, 75, 92, 131, 143, 0, & + 21, 38, 74, 102, 132, 159, 0, & + 18, 41, 57, 103, 133, 160, 0, & + 22, 46, 75, 100, 134, 161, 0, & + 16, 47, 76, 104, 113, 144, 0, & + 3, 48, 77, 93, 128, 161, 0, & + 6, 49, 78, 105, 127, 162, 0, & + 23, 42, 79, 96, 117, 149, 170, & + 20, 39, 80, 105, 122, 157, 0, & + 8, 35, 76, 106, 114, 146, 0, & + 21, 48, 70, 107, 126, 163, 0, & + 21, 50, 76, 85, 129, 158, 0, & + 24, 36, 81, 107, 131, 164, 0, & + 9, 51, 65, 100, 124, 141, 0, & + 22, 44, 82, 96, 135, 151, 0, & + 8, 52, 69, 90, 134, 165, 0, & + 25, 37, 79, 108, 136, 166, 0, & + 19, 49, 67, 108, 114, 152, 0, & + 22, 40, 60, 97, 136, 167, 0, & + 26, 53, 70, 83, 105, 150, 0, & + 27, 31, 81, 109, 123, 168, 0, & + 16, 30, 83, 101, 115, 155, 167, & + 7, 54, 58, 104, 131, 150, 173, & + 4, 45, 66, 99, 132, 152, 0, & + 15, 53, 84, 88, 132, 164, 0, & + 25, 30, 77, 110, 122, 137, 164, & + 23, 52, 67, 107, 133, 144, 174, & + 23, 40, 59, 111, 137, 140, 0, & + 14, 38, 68, 110, 120, 169, 0, & + 24, 35, 72, 102, 127, 153, 172, & + 28, 55, 84, 101, 125, 165, 174, & + 15, 42, 63, 95, 130, 169, 0, & + 4, 56, 85, 91, 134, 162, 170, & + 25, 55, 65, 89, 133, 156, 0, & + 5, 48, 63, 109, 135, 160, 0, & + 10, 53, 64, 87, 121, 140, 0, & + 9, 50, 82, 108, 138, 154, 0, & + 14, 54, 69, 99, 115, 157, 166, & + 17, 36, 56, 103, 130, 148, 0, & + 3, 47, 56, 73, 121, 154, 0, & + 1, 52, 77, 94, 139, 158, 0, & + 5, 28, 71, 110, 138, 147, 0, & + 27, 45, 68, 106, 124, 167, 0, & + 18, 51, 59, 98, 139, 162, 0, & + 2, 50, 80, 103, 120, 161, 0, & + 11, 46, 84, 86, 129, 166, 0, & + 20, 55, 73, 102, 119, 168, 0, & + 28, 43, 61, 111, 112, 146, 0, & + 11, 51, 83, 91, 135, 143, 0, & + 6, 46, 58, 111, 123, 148, 169, & + 26, 32, 79, 112, 116, 153, 0, & + 27, 49, 71, 89, 104, 163, 0, & + 12, 34, 78, 109, 117, 155, 0, & + 13, 47, 61, 94, 136, 159, 160, & + 17, 54, 82, 106, 137, 159, 0, & + 26, 44, 81, 90, 128, 156, 0, & + 7, 31, 74, 97, 138, 145, 0, & + 24, 37, 80, 98, 118, 145, 171/ + +data nrw/ & + 7,7,7,7,7,7,7,6,6,7,6,6,7,7,6,6,7,6, & + 6,6,6,6,7,6,6,6,6,6,6,6,6,6,7,6,6,6, & + 6,6,6,6,6,6,6,6,6,6,7,7,6,6,7,7,6,6, & + 7,7,6,7,6,6,6,6,7,6,6,6,6,6,6,6,6,6, & + 6,6,7,6,6,6,7,6,6,6,7/ + +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 bpdecode174_91 diff --git a/lib/fsk4hf/bpdecode174b.f90 b/lib/fsk4hf/bpdecode174b.f90 new file mode 100644 index 000000000..8291e7603 --- /dev/null +++ b/lib/fsk4hf/bpdecode174b.f90 @@ -0,0 +1,393 @@ +subroutine bpdecode174b(llr,apmask,maxiterations,decoded,cw,nharderror,iter) +! +! A log-domain belief propagation decoder for the (174,91) code. +! +integer, parameter:: N=174, K=91, M=N-K +integer*1 codeword(N),cw(N),apmask(N) +integer colorder(N) +integer*1 decoded(K) +integer Nm(7,M) +integer Mn(3,N) ! 3 checks per bit +integer synd(M) +real tov(3,N) +real toc(7,M) +real tanhtoc(7,M) +real zn(N) +real llr(N) +real Tmn +integer nrw(M) + +data colorder/ & + 0, 1, 2, 3, 28, 4, 5, 6, 7, 8, 9, 10, 11, 34, 12, 32, 13, 14, 15, 16,& + 17, 18, 36, 29, 40, 19, 20, 38, 21, 41, 30, 42, 22, 44, 37, 47, 48, 23, 33, 43,& + 49, 45, 56, 39, 25, 26, 46, 50, 51, 52, 24, 57, 58, 61, 31, 54, 64, 35, 27, 62,& + 59, 53, 60, 63, 55, 70, 66, 67, 68, 65, 71, 74, 72, 73, 77, 75, 69, 76, 79, 82,& + 83, 78, 81, 80, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,& + 100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,& + 120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,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/ + +data Mn/ & + 1, 24, 66, & + 2, 5, 70, & + 3, 31, 65, & + 4, 49, 58, & + 6, 60, 67, & + 7, 32, 75, & + 8, 48, 82, & + 9, 35, 41, & + 10, 39, 62, & + 11, 14, 61, & + 12, 71, 74, & + 13, 23, 78, & + 15, 16, 79, & + 17, 54, 63, & + 18, 50, 57, & + 19, 30, 47, & + 20, 64, 80, & + 21, 28, 69, & + 22, 25, 43, & + 26, 34, 72, & + 27, 36, 37, & + 29, 40, 44, & + 33, 52, 53, & + 38, 55, 83, & + 42, 51, 59, & + 45, 76, 81, & + 46, 68, 77, & + 56, 67, 73, & + 1, 4, 5, & + 2, 47, 51, & + 3, 46, 82, & + 6, 24, 76, & + 7, 9, 16, & + 8, 10, 78, & + 11, 35, 55, & + 12, 38, 64, & + 13, 42, 83, & + 14, 27, 54, & + 15, 21, 34, & + 17, 44, 53, & + 18, 25, 28, & + 19, 33, 57, & + 20, 22, 73, & + 23, 40, 81, & + 26, 49, 68, & + 29, 71, 75, & + 30, 65, 79, & + 31, 36, 60, & + 32, 43, 77, & + 37, 62, 70, & + 39, 69, 74, & + 41, 52, 66, & + 45, 50, 61, & + 48, 63, 80, & + 56, 59, 72, & + 58, 64, 65, & + 1, 13, 28, & + 2, 48, 75, & + 3, 53, 69, & + 4, 11, 44, & + 5, 73, 79, & + 6, 12, 17, & + 7, 57, 60, & + 8, 15, 61, & + 9, 39, 59, & + 10, 19, 49, & + 14, 43, 52, & + 16, 54, 68, & + 18, 41, 63, & + 20, 36, 45, & + 21, 67, 77, & + 10, 22, 55, & + 23, 65, 72, & + 24, 27, 82, & + 25, 26, 29, & + 30, 35, 37, & + 31, 51, 66, & + 17, 32, 78, & + 33, 42, 76, & + 34, 70, 83, & + 38, 46, 81, & + 40, 62, 80, & + 45, 47, 74, & + 50, 56, 71, & + 7, 37, 58, & + 1, 16, 71, & + 2, 6, 61, & + 3, 22, 50, & + 4, 59, 77, & + 5, 41, 81, & + 8, 58, 74, & + 9, 20, 26, & + 11, 21, 31, & + 12, 66, 79, & + 13, 14, 57, & + 15, 33, 40, & + 18, 44, 82, & + 19, 69, 83, & + 23, 49, 63, & + 24, 29, 39, & + 25, 47, 56, & + 27, 55, 72, & + 28, 64, 70, & + 30, 48, 77, & + 32, 34, 45, & + 35, 68, 80, & + 36, 38, 52, & + 42, 43, 62, & + 46, 60, 78, & + 51, 54, 67, & + 53, 73, 75, & + 14, 73, 76, & + 1, 22, 30, & + 2, 35, 43, & + 3, 47, 63, & + 4, 25, 76, & + 5, 33, 78, & + 6, 20, 83, & + 7, 12, 72, & + 8, 54, 70, & + 9, 61, 65, & + 10, 34, 51, & + 11, 46, 75, & + 13, 39, 68, & + 15, 17, 56, & + 16, 23, 36, & + 18, 32, 55, & + 19, 31, 81, & + 21, 37, 71, & + 24, 57, 64, & + 26, 38, 48, & + 27, 49, 50, & + 28, 52, 59, & + 29, 41, 58, & + 40, 60, 74, & + 42, 44, 79, & + 51, 53, 80, & + 62, 67, 82, & + 23, 66, 69, & + 1, 53, 61, & + 2, 18, 39, & + 3, 4, 12, & + 5, 26, 74, & + 6, 30, 52, & + 7, 82, 83, & + 8, 35, 73, & + 9, 19, 67, & + 10, 64, 75, & + 11, 20, 33, & + 13, 45, 48, & + 3, 14, 40, & + 15, 43, 49, & + 16, 55, 76, & + 17, 62, 65, & + 21, 47, 78, & + 22, 59, 81, & + 24, 34, 63, & + 25, 37, 66, & + 27, 79, 80, & + 28, 60, 79, & + 29, 31, 70, & + 32, 58, 69, & + 10, 36, 77, & + 38, 50, 51, & + 13, 41, 56, & + 42, 63, 71, & + 44, 47, 68, & + 1, 46, 72, & + 54, 57, 75, & + 2, 33, 58, & + 4, 17, 83, & + 5, 14, 55, & + 6, 23, 48, & + 7, 52, 56/ + +data Nm/ & + 1, 29, 57, 86, 113, 140, 168, & + 2, 30, 58, 87, 114, 141, 170, & + 3, 31, 59, 88, 115, 142, 151, & + 4, 29, 60, 89, 116, 142, 171, & + 2, 29, 61, 90, 117, 143, 172, & + 5, 32, 62, 87, 118, 144, 173, & + 6, 33, 63, 85, 119, 145, 174, & + 7, 34, 64, 91, 120, 146, 0, & + 8, 33, 65, 92, 121, 147, 0, & + 9, 34, 66, 72, 122, 148, 163, & + 10, 35, 60, 93, 123, 149, 0, & + 11, 36, 62, 94, 119, 142, 0, & + 12, 37, 57, 95, 124, 150, 165, & + 10, 38, 67, 95, 112, 151, 172, & + 13, 39, 64, 96, 125, 152, 0, & + 13, 33, 68, 86, 126, 153, 0, & + 14, 40, 62, 78, 125, 154, 171, & + 15, 41, 69, 97, 127, 141, 0, & + 16, 42, 66, 98, 128, 147, 0, & + 17, 43, 70, 92, 118, 149, 0, & + 18, 39, 71, 93, 129, 155, 0, & + 19, 43, 72, 88, 113, 156, 0, & + 12, 44, 73, 99, 126, 139, 173, & + 1, 32, 74, 100, 130, 157, 0, & + 19, 41, 75, 101, 116, 158, 0, & + 20, 45, 75, 92, 131, 143, 0, & + 21, 38, 74, 102, 132, 159, 0, & + 18, 41, 57, 103, 133, 160, 0, & + 22, 46, 75, 100, 134, 161, 0, & + 16, 47, 76, 104, 113, 144, 0, & + 3, 48, 77, 93, 128, 161, 0, & + 6, 49, 78, 105, 127, 162, 0, & + 23, 42, 79, 96, 117, 149, 170, & + 20, 39, 80, 105, 122, 157, 0, & + 8, 35, 76, 106, 114, 146, 0, & + 21, 48, 70, 107, 126, 163, 0, & + 21, 50, 76, 85, 129, 158, 0, & + 24, 36, 81, 107, 131, 164, 0, & + 9, 51, 65, 100, 124, 141, 0, & + 22, 44, 82, 96, 135, 151, 0, & + 8, 52, 69, 90, 134, 165, 0, & + 25, 37, 79, 108, 136, 166, 0, & + 19, 49, 67, 108, 114, 152, 0, & + 22, 40, 60, 97, 136, 167, 0, & + 26, 53, 70, 83, 105, 150, 0, & + 27, 31, 81, 109, 123, 168, 0, & + 16, 30, 83, 101, 115, 155, 167, & + 7, 54, 58, 104, 131, 150, 173, & + 4, 45, 66, 99, 132, 152, 0, & + 15, 53, 84, 88, 132, 164, 0, & + 25, 30, 77, 110, 122, 137, 164, & + 23, 52, 67, 107, 133, 144, 174, & + 23, 40, 59, 111, 137, 140, 0, & + 14, 38, 68, 110, 120, 169, 0, & + 24, 35, 72, 102, 127, 153, 172, & + 28, 55, 84, 101, 125, 165, 174, & + 15, 42, 63, 95, 130, 169, 0, & + 4, 56, 85, 91, 134, 162, 170, & + 25, 55, 65, 89, 133, 156, 0, & + 5, 48, 63, 109, 135, 160, 0, & + 10, 53, 64, 87, 121, 140, 0, & + 9, 50, 82, 108, 138, 154, 0, & + 14, 54, 69, 99, 115, 157, 166, & + 17, 36, 56, 103, 130, 148, 0, & + 3, 47, 56, 73, 121, 154, 0, & + 1, 52, 77, 94, 139, 158, 0, & + 5, 28, 71, 110, 138, 147, 0, & + 27, 45, 68, 106, 124, 167, 0, & + 18, 51, 59, 98, 139, 162, 0, & + 2, 50, 80, 103, 120, 161, 0, & + 11, 46, 84, 86, 129, 166, 0, & + 20, 55, 73, 102, 119, 168, 0, & + 28, 43, 61, 111, 112, 146, 0, & + 11, 51, 83, 91, 135, 143, 0, & + 6, 46, 58, 111, 123, 148, 169, & + 26, 32, 79, 112, 116, 153, 0, & + 27, 49, 71, 89, 104, 163, 0, & + 12, 34, 78, 109, 117, 155, 0, & + 13, 47, 61, 94, 136, 159, 160, & + 17, 54, 82, 106, 137, 159, 0, & + 26, 44, 81, 90, 128, 156, 0, & + 7, 31, 74, 97, 138, 145, 0, & + 24, 37, 80, 98, 118, 145, 171/ + +data nrw/ & + 7,7,7,7,7,7,7,6,6,7,6,6,7,7,6,6,7,6, & + 6,6,6,6,7,6,6,6,6,6,6,6,6,6,7,6,6,6, & + 6,6,6,6,6,6,6,6,6,6,7,7,6,6,7,7,6,6, & + 7,7,6,7,6,6,6,6,7,6,6,6,6,6,6,6,6,6, & + 6,6,7,6,6,6,7,6,6,6,7/ + +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 bpdecode174b diff --git a/lib/fsk4hf/encode174_91.f90 b/lib/fsk4hf/encode174_91.f90 new file mode 100644 index 000000000..afe193fa5 --- /dev/null +++ b/lib/fsk4hf/encode174_91.f90 @@ -0,0 +1,47 @@ +subroutine encode174_91(message,codeword) +! Encode an 91-bit message and return a 174-bit codeword. +! The generator matrix has dimensions (83,91). +! The code is a (174,91) regular ldpc code with column weight 3. +! + +include "ldpc_174_91_a_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,23 + read(g(i)(j:j),"(Z1)") istr + ibmax=4 + if(j.eq.23) ibmax=3 + 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 +itmp(1:M)=pchecks +itmp(M+1:N)=message(1:K) +codeword(colorder+1)=itmp(1:N) + +return +end subroutine encode174_91 diff --git a/lib/fsk4hf/ldpc_174_91_a_params.f90 b/lib/fsk4hf/ldpc_174_91_a_params.f90 new file mode 100644 index 000000000..4001d5746 --- /dev/null +++ b/lib/fsk4hf/ldpc_174_91_a_params.f90 @@ -0,0 +1,100 @@ +integer, parameter:: N=174, K=91, M=N-K +character*22 g(83) +integer colorder(N) + +data g/ & + "2a6a9ab98f661b797baa21a", & + "5fda604488977fdc30ff630", & + "8a450007032409e8797b13a", & + "0f7923bfa5d559590ef6efe", & + "44dc14bc4461e645f847c78", & + "f8c66224febd3e60f5e3708", & + "6d60329e83fb0e1b1bdac2e", & + "66435a472a7837a0e5d4e12", & + "9d0feced8745a66e328c310", & + "1791b0e7c5eaa43710c4276", & + "e5cbd2d5b4d65d3a432d97e", & + "c2241f8795e0e5bc6bd9052", & + "222d861201a4697c2689576", & + "aa2ee5d6d462e206f59cbe8", & + "e486eb73894e6a0964d8c40", & + "4099d5b42d36301cff6dbd6", & + "40c50b9341f7b5ea08dabde", & + "c90359074895363d428f072", & + "ca819cb6569fbfe26b68ef8", & + "4d983341fb56b8e1dae3450", & + "2dce341bc8fd0e5de04fa52", & + "3e7b01b376e3e5f6080de0e", & + "6c8b0813ca2394c08564f94", & + "c322ca8ea866784adc9451a", & + "6378aa1a03fab3e163aa4b0", & + "3c92ea8df0003883a021d70", & + "c793729067176eca26b83c2", & + "d3fae76046a36dff711207a", & + "bc9bf3ef57137fda1c325da", & + "a4eabe2df65a083ea6387c8", & + "650e3da3a0c0349154131d8", & + "1fb4c59ffc11c648ad06760", & + "1471f9599543f13fd7eb6ae", & + "6111012405186e84cba67ce", & + "c4da3574edafefff976fc08", & + "953f854e40701063115c0f2", & + "1f7ae6982f9a5733c44fb70", & + "83e101fe5e80c1b8541728e", & + "50375654edd53054f81e228", & + "1bb03a21a6cde34dff7ec96", & + "b0b279a934342aa0e188b3a", & + "e1989846a20a09cd77b1f64", & + "4eb68e01cb07fdbc83edee2", & + "f33ac4ec36a7c8e6ea8364c", & + "99b03a21a6c5e34dfffec96", & + "e50e3de3a1c034915413158", & + "fda09f8b05b8fb80ac78600", & + "ca8709be6b193204dd25ab0", & + "35701ff0cc3a03f213a93d2", & + "c2bfdec67f7b5a4c5ee7544", & + "dc184fe7e93a65c1b4b7cd2", & + "8cf8aac820f107d6ec6b30a", & + "e74b3da5a3e43d593d680e2", & + "c1e51f79db6124243fceadc", & + "29237d5d05dc1a4cca2ddd0", & + "050e76be4749b3b279d6414", & + "dd163959ae739673cde18c6", & + "03e100fe5e81c1b85417a8e", & + "06b2b17f70e75fc365bed20", & + "6df9e72abecd3e03e4b77fa", & + "4fa5370361b4bf3cf6b1296", & + "eabbf88f0a88307629bfd1c", & + "190674f88cf69989c8b8a40", & + "37740c13cfad07f61dcac3a", & + "4e7923bfa5d579590ef4ede", & + "fe74d37b8e5a63a2905da28", & + "2101e7a95979b2c5c44257e", & + "841f3ec7a4585a159fb5796", & + "aa7ff31d4b7f859c21254c2", & + "6e69229ba0cdb7ddcd50930", & + "29cfc4288af223bea58b96e", & + "5d03eba9f51956176b87abe", & + "399cbc33a7498b31d9f79e4", & + "034967e48ab80135b1c7fca", & + "721ad006ac715928df9775e", & + "37210b395327446ac7108f8", & + "52acf6de27477ea937e5330", & + "1f3a8549435c198b68231c8", & + "ef6809edb4a3557cd173d0a", & + "09a31639fef9c7a8b6fcae2", & + "03bc87c137eeec711c68d36", & + "b09347742319f90131d3146", & + "a723c9cef1de8c97f34c94c"/ + +data colorder/ & + 0, 1, 2, 3, 28, 4, 5, 6, 7, 8, 9, 10, 11, 34, 12, 32, 13, 14, 15, 16,& + 17, 18, 36, 29, 40, 19, 20, 38, 21, 41, 30, 42, 22, 44, 37, 47, 48, 23, 33, 43,& + 49, 45, 56, 39, 25, 26, 46, 50, 51, 52, 24, 57, 58, 61, 31, 54, 64, 35, 27, 62,& + 59, 53, 60, 63, 55, 70, 66, 67, 68, 65, 71, 74, 72, 73, 77, 75, 69, 76, 79, 82,& + 83, 78, 81, 80, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,& + 100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,& + 120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,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/ + diff --git a/lib/fsk4hf/ldpcsim174_91.f90 b/lib/fsk4hf/ldpcsim174_91.f90 new file mode 100644 index 000000000..5276fd8c3 --- /dev/null +++ b/lib/fsk4hf/ldpcsim174_91.f90 @@ -0,0 +1,234 @@ +program ldpcsim174_91 +! End to end test of the (174,77)/crc14 encoder and decoder. +use crc +use packjt + +character*22 msg,msgsent,msgreceived +character*8 arg +character*6 grid +integer*1, allocatable :: codeword(:), decoded(:), message(:) +integer*1, target:: i1Msg8BitBytes(11) +integer*1 msgbits(91) +integer*1 apmask(174), cw(174) +integer*2 checksum +integer*4 i4Msg6BitWords(13) +integer colorder(174) +integer nerrtot(174),nerrdec(174),nmpcbad(91) +logical checksumok,fsk,bpsk +real*8, allocatable :: rxdata(:) +real, allocatable :: llr(:) + +data colorder/ & + 0, 1, 2, 3, 28, 4, 5, 6, 7, 8, 9, 10, 11, 34, 12, 32, 13, 14, 15, 16,& + 17, 18, 36, 29, 40, 19, 20, 38, 21, 41, 30, 42, 22, 44, 37, 47, 48, 23, 33, 43,& + 49, 45, 56, 39, 25, 26, 46, 50, 51, 52, 24, 57, 58, 61, 31, 54, 64, 35, 27, 62,& + 59, 53, 60, 63, 55, 70, 66, 67, 68, 65, 71, 74, 72, 73, 77, 75, 69, 76, 79, 82,& + 83, 78, 81, 80, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,& + 100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,& + 120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,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/ + +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 ndepth #trials s ' + print*,'eg: ldpcsim 10 2 1000 0.84' + print*,'belief propagation iterations: niter, ordered-statistics depth: ndepth' + 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,*) ndepth +call getarg(3,arg) +read(arg,*) ntrials +call getarg(4,arg) +read(arg,*) s + +fsk=.false. +bpsk=.true. + +! don't count crc bits as data bits +N=174 +K=91 +! scale Eb/No for a (174,91) code +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) ) + + msg="K1JT K9AN EN50" +! msg="G4WJS K9AN EN50" + call packmsg(msg,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes + call unpackmsg(i4Msg6BitWords,msgsent,.false.,grid) !Unpack to get msgsent + write(*,*) "message sent ",msgsent + + i4=0 + ik=0 + im=0 + do i=1,12 + nn=i4Msg6BitWords(i) + do j=1, 6 + ik=ik+1 + i4=i4+i4+iand(1,ishft(nn,j-6)) + i4=iand(i4,255) + if(ik.eq.8) then + im=im+1 +! if(i4.gt.127) i4=i4-256 + i1Msg8BitBytes(im)=i4 + ik=0 + endif + enddo + enddo + + i1Msg8BitBytes(10:11)=0 + checksum = crc14 (c_loc (i1Msg8BitBytes), 11) +! For reference, the next 3 lines show how to check the CRC + i1Msg8BitBytes(10)=checksum/256 + i1Msg8BitBytes(11)=iand (checksum,255) + checksumok = crc14_check(c_loc (i1Msg8BitBytes), 11) + if( checksumok ) write(*,*) 'Good checksum' + +! K=87, For now: +! msgbits(1:72) JT message bits +! msgbits(73:75) 3 free message bits (set to 0) +! msgbits(76:87) CRC12 + mbit=0 + do i=1, 9 + i1=i1Msg8BitBytes(i) + do ibit=1,8 + mbit=mbit+1 + msgbits(mbit)=iand(1,ishft(i1,ibit-8)) + enddo + enddo + msgbits(73:75)=0 ! the three extra message bits go here + i1=i1Msg8BitBytes(10) ! First 4 bits of crc12 are LSB of this byte + do ibit=1,4 + msgbits(75+ibit)=iand(1,ishft(i1,ibit-4)) + enddo + i1=i1Msg8BitBytes(11) ! Now shift in last 8 bits of the CRC + do ibit=1,8 + msgbits(79+ibit)=iand(1,ishft(i1,ibit-8)) + enddo + + write(*,*) 'message' + write(*,'(11(8i1,1x))') msgbits + + call encode174(msgbits,codeword) + call init_random_seed() +! call sgran() + + write(*,*) 'codeword' + write(*,'(22(8i1,1x))') codeword + +write(*,*) "Es/N0 SNR2500 ngood nundetected nbadcrc sigma" +do idb = 20,-10,-1 +!do idb = -3,-3,-1 + db=idb/2.0-1.0 + sigma=1/sqrt( 2*(10**(db/10.0)) ) + 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) + nap=0 ! number of AP bits + llr(colorder(174-87+1:174-87+nap)+1)=5*(2.0*msgbits(1:nap)-1.0) + apmask=0 + apmask(colorder(174-87+1:174-87+nap)+1)=1 + +! max_iterations is max number of belief propagation iterations + call bpdecode174(llr, apmask, max_iterations, decoded, cw, nharderrors,niterations) + if( ndepth .ge. 0 .and. nharderrors .lt. 0 ) call osd174(llr, apmask, ndepth, decoded, cw, nharderrors, dmin) +! If the decoder finds a valid codeword, nharderrors will be .ge. 0. + if( nharderrors .ge. 0 ) then + call extractmessage174(decoded,msgreceived,ncrcflag) + 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 + if( ncrcflag .eq. 1 ) then + if( nueflag .eq. 0 ) then + ngood=ngood+1 + if(nerr.ge.1) nerrdec(nerr)=nerrdec(nerr)+1 + else if( nueflag .eq. 1 ) then + nue=nue+1; + endif + endif + endif + enddo + baud=12000/1920 + snr2500=db+10.0*log10((baud/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,174 + 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,87 + write(25,'(i4,2x,i10)') i,nmpcbad(i) +enddo +close(25) + +end program ldpcsim174_91 diff --git a/lib/fsk4hf/osd174_91.f90 b/lib/fsk4hf/osd174_91.f90 new file mode 100644 index 000000000..091214444 --- /dev/null +++ b/lib/fsk4hf/osd174_91.f90 @@ -0,0 +1,365 @@ +subroutine osd174_91(llr,apmask,ndeep,decoded,cw,nhardmin,dmin) +! +! An ordered-statistics decoder for the (174,91) code. +! +include "ldpc_174_91_a_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,22 + 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 osd174 + +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 +