WSJT-X/lib/bpdecode144.f90
Steven Franke 808259b7bf Add early termination code to the msk144 long-message decoder. Remove references to Radford Neal's ldpc code, which served us well but is no longer needed.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7051 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2016-09-01 21:50:52 +00:00

349 lines
7.9 KiB
Fortran

subroutine pltanh(x,y)
isign=+1
z=x
if( x.lt.0 ) then
isign=-1
z=abs(x)
endif
if( z.le. 0.8 ) then
y=0.83*x
return
elseif( z.le. 1.6 ) then
y=isign*(0.322*z+0.4064)
return
elseif( z.le. 3.0 ) then
y=isign*(0.0524*z+0.8378)
return
elseif( z.lt. 7.0 ) then
y=isign*(0.0012*z+0.9914)
return
else
y=isign*0.9998
return
endif
end subroutine pltanh
subroutine platanh(x,y)
isign=+1
z=x
if( x.lt.0 ) then
isign=-1
z=abs(x)
endif
if( z.le. 0.664 ) then
y=x/0.83
return
elseif( z.le. 0.9217 ) then
y=isign*(z-0.4064)/0.322
return
elseif( z.le. 0.9951 ) then
y=isign*(z-0.8378)/0.0524
return
elseif( z.le. 0.9998 ) then
y=isign*(z-0.9914)/0.0012
return
else
y=isign*7.0
return
endif
end subroutine platanh
subroutine bpdecode144(llr,maxiterations,decoded,niterations)
!
! A log-domain belief propagation decoder for the msk144 code.
! The code is a regular (128,80) code with column weight 3 and row weight 8.
! k9an August, 2016
!
integer, parameter:: N=128, K=80, M=N-K
integer*1 codeword(N),cw(N)
integer*1 colorder(N)
integer*1 decoded(K)
integer Nm(8,M) ! 8 bits per check
integer Mn(3,N) ! 3 checks per bit
integer synd(M)
real tov(3,N) ! single precision seems to be adequate in log-domain
real toc(8,M)
real tanhtoc(8,M)
real zn(N)
real llr(N)
real Tmn
data colorder/0,1,2,3,4,5,6,7,8,9, &
10,11,12,13,14,15,24,26,29,30, &
32,43,44,47,60,77,79,97,101,111, &
96,38,64,53,93,34,59,94,74,90, &
108,123,85,57,70,25,69,62,48,49, &
50,51,52,33,54,55,56,21,58,36, &
16,61,23,63,20,65,66,67,68,46, &
22,71,72,73,31,75,76,45,78,17, &
80,81,82,83,84,42,86,87,88,89, &
39,91,92,35,37,95,19,27,98,99, &
100,28,102,103,104,105,106,107,40,109, &
110,18,112,113,114,115,116,117,118,119, &
120,121,122,41,124,125,126,127/
data Mn/ &
1, 14, 38, &
2, 4, 41, &
3, 19, 39, &
5, 29, 34, &
6, 35, 40, &
7, 20, 45, &
8, 28, 48, &
9, 22, 25, &
10, 24, 36, &
11, 12, 37, &
13, 43, 44, &
15, 18, 46, &
16, 17, 47, &
21, 32, 33, &
23, 30, 31, &
26, 27, 42, &
1, 12, 46, &
2, 36, 38, &
3, 5, 10, &
4, 9, 23, &
6, 13, 39, &
7, 15, 17, &
8, 18, 27, &
11, 33, 40, &
14, 28, 44, &
16, 29, 31, &
19, 20, 22, &
21, 30, 42, &
24, 26, 47, &
25, 37, 48, &
32, 34, 45, &
8, 35, 41, &
12, 31, 43, &
1, 19, 21, &
2, 43, 45, &
3, 4, 11, &
5, 18, 33, &
6, 25, 47, &
7, 28, 30, &
9, 14, 34, &
10, 35, 42, &
13, 15, 22, &
16, 37, 38, &
17, 41, 44, &
20, 24, 29, &
18, 23, 39, &
12, 26, 32, &
27, 38, 40, &
15, 36, 48, &
2, 30, 46, &
1, 4, 13, &
3, 28, 32, &
5, 43, 47, &
6, 34, 46, &
7, 9, 40, &
8, 11, 45, &
10, 17, 23, &
14, 31, 35, &
16, 22, 42, &
19, 37, 44, &
20, 33, 48, &
21, 24, 41, &
25, 27, 29, &
26, 39, 48, &
19, 31, 36, &
1, 5, 7, &
2, 29, 39, &
3, 16, 46, &
4, 26, 37, &
6, 28, 45, &
8, 22, 33, &
9, 21, 43, &
10, 25, 38, &
11, 14, 24, &
12, 17, 40, &
13, 27, 30, &
15, 32, 35, &
18, 44, 47, &
20, 23, 36, &
34, 41, 42, &
1, 32, 48, &
2, 3, 33, &
4, 29, 42, &
5, 14, 37, &
6, 7, 36, &
8, 9, 39, &
10, 13, 19, &
11, 18, 30, &
12, 16, 20, &
15, 29, 44, &
17, 34, 38, &
6, 21, 22, &
23, 32, 40, &
24, 27, 46, &
25, 41, 45, &
7, 26, 43, &
28, 31, 47, &
20, 35, 38, &
1, 33, 41, &
2, 42, 44, &
3, 23, 48, &
4, 31, 45, &
5, 8, 30, &
9, 16, 36, &
10, 40, 47, &
11, 17, 46, &
12, 21, 34, &
13, 24, 28, &
14, 18, 43, &
15, 25, 26, &
19, 27, 35, &
22, 37, 39, &
1, 16, 18, &
2, 6, 20, &
3, 30, 43, &
4, 28, 33, &
5, 22, 23, &
7, 39, 42, &
8, 12, 38, &
9, 35, 46, &
10, 27, 32, &
11, 15, 34, &
13, 36, 37, &
14, 41, 47, &
17, 21, 25, &
19, 29, 45, &
24, 31, 48, &
26, 40, 44/
data Nm/ &
1, 17, 34, 51, 66, 81, 99, 113, &
2, 18, 35, 50, 67, 82, 100, 114, &
3, 19, 36, 52, 68, 82, 101, 115, &
2, 20, 36, 51, 69, 83, 102, 116, &
4, 19, 37, 53, 66, 84, 103, 117, &
5, 21, 38, 54, 70, 85, 92, 114, &
6, 22, 39, 55, 66, 85, 96, 118, &
7, 23, 32, 56, 71, 86, 103, 119, &
8, 20, 40, 55, 72, 86, 104, 120, &
9, 19, 41, 57, 73, 87, 105, 121, &
10, 24, 36, 56, 74, 88, 106, 122, &
10, 17, 33, 47, 75, 89, 107, 119, &
11, 21, 42, 51, 76, 87, 108, 123, &
1, 25, 40, 58, 74, 84, 109, 124, &
12, 22, 42, 49, 77, 90, 110, 122, &
13, 26, 43, 59, 68, 89, 104, 113, &
13, 22, 44, 57, 75, 91, 106, 125, &
12, 23, 37, 46, 78, 88, 109, 113, &
3, 27, 34, 60, 65, 87, 111, 126, &
6, 27, 45, 61, 79, 89, 98, 114, &
14, 28, 34, 62, 72, 92, 107, 125, &
8, 27, 42, 59, 71, 92, 112, 117, &
15, 20, 46, 57, 79, 93, 101, 117, &
9, 29, 45, 62, 74, 94, 108, 127, &
8, 30, 38, 63, 73, 95, 110, 125, &
16, 29, 47, 64, 69, 96, 110, 128, &
16, 23, 48, 63, 76, 94, 111, 121, &
7, 25, 39, 52, 70, 97, 108, 116, &
4, 26, 45, 63, 67, 83, 90, 126, &
15, 28, 39, 50, 76, 88, 103, 115, &
15, 26, 33, 58, 65, 97, 102, 127, &
14, 31, 47, 52, 77, 81, 93, 121, &
14, 24, 37, 61, 71, 82, 99, 116, &
4, 31, 40, 54, 80, 91, 107, 122, &
5, 32, 41, 58, 77, 98, 111, 120, &
9, 18, 49, 65, 79, 85, 104, 123, &
10, 30, 43, 60, 69, 84, 112, 123, &
1, 18, 43, 48, 73, 91, 98, 119, &
3, 21, 46, 64, 67, 86, 112, 118, &
5, 24, 48, 55, 75, 93, 105, 128, &
2, 32, 44, 62, 80, 95, 99, 124, &
16, 28, 41, 59, 80, 83, 100, 118, &
11, 33, 35, 53, 72, 96, 109, 115, &
11, 25, 44, 60, 78, 90, 100, 128, &
6, 31, 35, 56, 70, 95, 102, 126, &
12, 17, 50, 54, 68, 94, 106, 120, &
13, 29, 38, 53, 78, 97, 105, 124, &
7, 30, 49, 61, 64, 81, 101, 127/
nrw=8
ncw=3
toc=0
tov=0
tanhtoc=0
! initial messages to checks
do j=1,M
do i=1,nrw
toc(i,j)=llr((Nm(i,j)))
enddo
enddo
ncnt=0
do iter=0,maxiterations
! Update bit log likelihood ratios
do i=1,N
zn(i)=llr(i)+sum(tov(1:ncw,i))
enddo
! Check to see if we have a codeword
cw=0
where( zn .gt. 0. ) cw=1
ncheck=0
do i=1,M
synd(i)=sum(cw(Nm(:,i)))
if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
enddo
if( ncheck .eq. 0 ) then ! we have a codeword
niterations=iter
codeword=cw(colorder+1)
decoded=codeword(M+1:N)
return
endif
if( iter.gt.0 ) then ! this code block implements an early stopping criterion
nd=ncheck-nclast
if( nd .lt. 0 ) then ! # of unsatisfied parity checks decreased
ncnt=0 ! reset counter
else
ncnt=ncnt+1
endif
! write(*,*) iter,ncheck,nd,ncnt
if( ncnt .ge. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) then
niterations=-1
return
endif
endif
nclast=ncheck
! Send messages from bits to check nodes
do j=1,M
do i=1,nrw
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 ! Mn(3,128)
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:nrw,i)=tanh(-toc(1:nrw,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(:,ichk),mask=Nm(:,ichk).ne.j)
call platanh(-Tmn,y)
tov(i,j)=2*y
enddo
enddo
enddo
niterations=-1
end subroutine bpdecode144