1subroutine pltanh(x,y)
2  isign=+1
3  z=x
4  if( x.lt.0 ) then
5    isign=-1
6    z=abs(x)
7  endif
8  if( z.le. 0.8 ) then
9    y=0.83*x
10    return
11  elseif( z.le. 1.6 ) then
12    y=isign*(0.322*z+0.4064)
13    return
14  elseif( z.le. 3.0 ) then
15    y=isign*(0.0524*z+0.8378)
16    return
17  elseif( z.lt. 7.0 ) then
18    y=isign*(0.0012*z+0.9914)
19    return
20  else
21    y=isign*0.9998
22    return
23  endif
24end subroutine pltanh
25
26subroutine platanh(x,y)
27  isign=+1
28  z=x
29  if( x.lt.0 ) then
30    isign=-1
31    z=abs(x)
32  endif
33  if( z.le. 0.664 ) then
34    y=x/0.83
35    return
36  elseif( z.le. 0.9217 ) then
37    y=isign*(z-0.4064)/0.322
38    return
39  elseif( z.le. 0.9951 ) then
40    y=isign*(z-0.8378)/0.0524
41    return
42  elseif( z.le. 0.9998 ) then
43    y=isign*(z-0.9914)/0.0012
44    return
45  else
46    y=isign*7.0
47    return
48  endif
49end subroutine platanh
50
51subroutine bpdecode144(llr,maxiterations,decoded,niterations)
52!
53! A log-domain belief propagation decoder for the msk144 code.
54! The code is a regular (128,80) code with column weight 3 and row weight 8.
55! k9an August, 2016
56!
57integer, parameter:: N=128, K=80, M=N-K
58integer*1 codeword(N),cw(N)
59integer*1 colorder(N)
60integer*1 decoded(K)
61integer Nm(8,M)  ! 8 bits per check
62integer Mn(3,N)  ! 3 checks per bit
63integer synd(M)
64real tov(3,N)    ! single precision seems to be adequate in log-domain
65real toc(8,M)
66real tanhtoc(8,M)
67real zn(N)
68real llr(N)
69real Tmn
70
71data colorder/0,1,2,3,4,5,6,7,8,9, &
72              10,11,12,13,14,15,24,26,29,30, &
73              32,43,44,47,60,77,79,97,101,111, &
74              96,38,64,53,93,34,59,94,74,90, &
75              108,123,85,57,70,25,69,62,48,49, &
76              50,51,52,33,54,55,56,21,58,36, &
77              16,61,23,63,20,65,66,67,68,46, &
78              22,71,72,73,31,75,76,45,78,17, &
79              80,81,82,83,84,42,86,87,88,89, &
80              39,91,92,35,37,95,19,27,98,99, &
81              100,28,102,103,104,105,106,107,40,109, &
82              110,18,112,113,114,115,116,117,118,119, &
83              120,121,122,41,124,125,126,127/
84
85data Mn/               &
86   1,  14,  38, &
87   2,   4,  41, &
88   3,  19,  39, &
89   5,  29,  34, &
90   6,  35,  40, &
91   7,  20,  45, &
92   8,  28,  48, &
93   9,  22,  25, &
94  10,  24,  36, &
95  11,  12,  37, &
96  13,  43,  44, &
97  15,  18,  46, &
98  16,  17,  47, &
99  21,  32,  33, &
100  23,  30,  31, &
101  26,  27,  42, &
102   1,  12,  46, &
103   2,  36,  38, &
104   3,   5,  10, &
105   4,   9,  23, &
106   6,  13,  39, &
107   7,  15,  17, &
108   8,  18,  27, &
109  11,  33,  40, &
110  14,  28,  44, &
111  16,  29,  31, &
112  19,  20,  22, &
113  21,  30,  42, &
114  24,  26,  47, &
115  25,  37,  48, &
116  32,  34,  45, &
117   8,  35,  41, &
118  12,  31,  43, &
119   1,  19,  21, &
120   2,  43,  45, &
121   3,   4,  11, &
122   5,  18,  33, &
123   6,  25,  47, &
124   7,  28,  30, &
125   9,  14,  34, &
126  10,  35,  42, &
127  13,  15,  22, &
128  16,  37,  38, &
129  17,  41,  44, &
130  20,  24,  29, &
131  18,  23,  39, &
132  12,  26,  32, &
133  27,  38,  40, &
134  15,  36,  48, &
135   2,  30,  46, &
136   1,   4,  13, &
137   3,  28,  32, &
138   5,  43,  47, &
139   6,  34,  46, &
140   7,   9,  40, &
141   8,  11,  45, &
142  10,  17,  23, &
143  14,  31,  35, &
144  16,  22,  42, &
145  19,  37,  44, &
146  20,  33,  48, &
147  21,  24,  41, &
148  25,  27,  29, &
149  26,  39,  48, &
150  19,  31,  36, &
151   1,   5,   7, &
152   2,  29,  39, &
153   3,  16,  46, &
154   4,  26,  37, &
155   6,  28,  45, &
156   8,  22,  33, &
157   9,  21,  43, &
158  10,  25,  38, &
159  11,  14,  24, &
160  12,  17,  40, &
161  13,  27,  30, &
162  15,  32,  35, &
163  18,  44,  47, &
164  20,  23,  36, &
165  34,  41,  42, &
166   1,  32,  48, &
167   2,   3,  33, &
168   4,  29,  42, &
169   5,  14,  37, &
170   6,   7,  36, &
171   8,   9,  39, &
172  10,  13,  19, &
173  11,  18,  30, &
174  12,  16,  20, &
175  15,  29,  44, &
176  17,  34,  38, &
177   6,  21,  22, &
178  23,  32,  40, &
179  24,  27,  46, &
180  25,  41,  45, &
181   7,  26,  43, &
182  28,  31,  47, &
183  20,  35,  38, &
184   1,  33,  41, &
185   2,  42,  44, &
186   3,  23,  48, &
187   4,  31,  45, &
188   5,   8,  30, &
189   9,  16,  36, &
190  10,  40,  47, &
191  11,  17,  46, &
192  12,  21,  34, &
193  13,  24,  28, &
194  14,  18,  43, &
195  15,  25,  26, &
196  19,  27,  35, &
197  22,  37,  39, &
198   1,  16,  18, &
199   2,   6,  20, &
200   3,  30,  43, &
201   4,  28,  33, &
202   5,  22,  23, &
203   7,  39,  42, &
204   8,  12,  38, &
205   9,  35,  46, &
206  10,  27,  32, &
207  11,  15,  34, &
208  13,  36,  37, &
209  14,  41,  47, &
210  17,  21,  25, &
211  19,  29,  45, &
212  24,  31,  48, &
213  26,  40,  44/
214
215data Nm/               &
216   1,  17,  34,  51,  66,  81,  99, 113, &
217   2,  18,  35,  50,  67,  82, 100, 114, &
218   3,  19,  36,  52,  68,  82, 101, 115, &
219   2,  20,  36,  51,  69,  83, 102, 116, &
220   4,  19,  37,  53,  66,  84, 103, 117, &
221   5,  21,  38,  54,  70,  85,  92, 114, &
222   6,  22,  39,  55,  66,  85,  96, 118, &
223   7,  23,  32,  56,  71,  86, 103, 119, &
224   8,  20,  40,  55,  72,  86, 104, 120, &
225   9,  19,  41,  57,  73,  87, 105, 121, &
226  10,  24,  36,  56,  74,  88, 106, 122, &
227  10,  17,  33,  47,  75,  89, 107, 119, &
228  11,  21,  42,  51,  76,  87, 108, 123, &
229   1,  25,  40,  58,  74,  84, 109, 124, &
230  12,  22,  42,  49,  77,  90, 110, 122, &
231  13,  26,  43,  59,  68,  89, 104, 113, &
232  13,  22,  44,  57,  75,  91, 106, 125, &
233  12,  23,  37,  46,  78,  88, 109, 113, &
234   3,  27,  34,  60,  65,  87, 111, 126, &
235   6,  27,  45,  61,  79,  89,  98, 114, &
236  14,  28,  34,  62,  72,  92, 107, 125, &
237   8,  27,  42,  59,  71,  92, 112, 117, &
238  15,  20,  46,  57,  79,  93, 101, 117, &
239   9,  29,  45,  62,  74,  94, 108, 127, &
240   8,  30,  38,  63,  73,  95, 110, 125, &
241  16,  29,  47,  64,  69,  96, 110, 128, &
242  16,  23,  48,  63,  76,  94, 111, 121, &
243   7,  25,  39,  52,  70,  97, 108, 116, &
244   4,  26,  45,  63,  67,  83,  90, 126, &
245  15,  28,  39,  50,  76,  88, 103, 115, &
246  15,  26,  33,  58,  65,  97, 102, 127, &
247  14,  31,  47,  52,  77,  81,  93, 121, &
248  14,  24,  37,  61,  71,  82,  99, 116, &
249   4,  31,  40,  54,  80,  91, 107, 122, &
250   5,  32,  41,  58,  77,  98, 111, 120, &
251   9,  18,  49,  65,  79,  85, 104, 123, &
252  10,  30,  43,  60,  69,  84, 112, 123, &
253   1,  18,  43,  48,  73,  91,  98, 119, &
254   3,  21,  46,  64,  67,  86, 112, 118, &
255   5,  24,  48,  55,  75,  93, 105, 128, &
256   2,  32,  44,  62,  80,  95,  99, 124, &
257  16,  28,  41,  59,  80,  83, 100, 118, &
258  11,  33,  35,  53,  72,  96, 109, 115, &
259  11,  25,  44,  60,  78,  90, 100, 128, &
260   6,  31,  35,  56,  70,  95, 102, 126, &
261  12,  17,  50,  54,  68,  94, 106, 120, &
262  13,  29,  38,  53,  78,  97, 105, 124, &
263   7,  30,  49,  61,  64,  81, 101, 127/
264
265nrw=8
266ncw=3
267
268toc=0
269tov=0
270tanhtoc=0
271
272! initial messages to checks
273do j=1,M
274  do i=1,nrw
275    toc(i,j)=llr((Nm(i,j)))
276  enddo
277enddo
278
279ncnt=0
280
281do iter=0,maxiterations
282
283! Update bit log likelihood ratios
284  do i=1,N
285    zn(i)=llr(i)+sum(tov(1:ncw,i))
286  enddo
287
288! Check to see if we have a codeword
289  cw=0
290  where( zn .gt. 0. ) cw=1
291  ncheck=0
292  do i=1,M
293    synd(i)=sum(cw(Nm(:,i)))
294    if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
295  enddo
296
297  if( ncheck .eq. 0 ) then ! we have a codeword
298    niterations=iter
299    codeword=cw(colorder+1)
300    decoded=codeword(M+1:N)
301    return
302  endif
303
304  if( iter.gt.0 ) then  ! this code block implements an early stopping criterion
305    nd=ncheck-nclast
306    if( nd .lt. 0 ) then ! # of unsatisfied parity checks decreased
307      ncnt=0  ! reset counter
308    else
309      ncnt=ncnt+1
310    endif
311!    write(*,*) iter,ncheck,nd,ncnt
312    if( ncnt .ge. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) then
313      niterations=-1
314      return
315    endif
316  endif
317  nclast=ncheck
318
319! Send messages from bits to check nodes
320  do j=1,M
321    do i=1,nrw
322      ibj=Nm(i,j)
323      toc(i,j)=zn(ibj)
324      do kk=1,ncw ! subtract off what the bit had received from the check
325        if( Mn(kk,ibj) .eq. j ) then  ! Mn(3,128)
326          toc(i,j)=toc(i,j)-tov(kk,ibj)
327        endif
328      enddo
329    enddo
330  enddo
331
332! send messages from check nodes to variable nodes
333  do i=1,M
334    tanhtoc(1:nrw,i)=tanh(-toc(1:nrw,i)/2)
335  enddo
336
337  do j=1,N
338    do i=1,ncw
339      ichk=Mn(i,j)  ! Mn(:,j) are the checks that include bit j
340      Tmn=product(tanhtoc(:,ichk),mask=Nm(:,ichk).ne.j)
341      call platanh(-Tmn,y)
342      tov(i,j)=2*y
343    enddo
344  enddo
345
346enddo
347niterations=-1
348end subroutine bpdecode144
349