From 4565b182ff37f70630b4fdb70320578faffa257c Mon Sep 17 00:00:00 2001
From: Steven Franke <s.j.franke@icloud.com>
Date: Tue, 28 Mar 2017 12:05:49 +0000
Subject: [PATCH] Add some experimental routines.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7621 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
---
 CMakeLists.txt                   |   8 +
 lib/fsk4hf/bpdecode168.f90       | 380 +++++++++++++++++++++++++++++++
 lib/fsk4hf/encode168.f90         | 141 ++++++++++++
 lib/fsk4hf/extractmessage168.f90 |  48 ++++
 lib/fsk4hf/ldpcsim168.f90        | 219 ++++++++++++++++++
 5 files changed, 796 insertions(+)
 create mode 100644 lib/fsk4hf/bpdecode168.f90
 create mode 100644 lib/fsk4hf/encode168.f90
 create mode 100644 lib/fsk4hf/extractmessage168.f90
 create mode 100644 lib/fsk4hf/ldpcsim168.f90

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6e5e8d27d..996960637 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -309,6 +309,7 @@ set (wsjt_FSRCS
   lib/badmsg.f90
   lib/bpdecode40.f90
   lib/bpdecode144.f90
+  lib/fsk4hf/bpdecode168.f90
   lib/baddata.f90
   lib/ccf2.f90
   lib/ccf65.f90
@@ -316,6 +317,7 @@ set (wsjt_FSRCS
   lib/chkmsg.f90
   lib/chkss2.f90
   lib/coord.f90
+  lib/crc.f90
   lib/db.f90
   lib/decode4.f90
   lib/decode65a.f90
@@ -332,11 +334,13 @@ set (wsjt_FSRCS
   lib/encode4.f90
   lib/encode_msk40.f90
   lib/encode_msk144.f90
+  lib/fsk4hf/encode168.f90
   lib/entail.f90
   lib/ephem.f90
   lib/extract.f90
   lib/extract4.f90
   lib/extractmessage144.f90
+  lib/fsk4hf/extractmessage168.f90
   lib/fano232.f90
   lib/fast9.f90
   lib/fast_decode.f90
@@ -396,6 +400,7 @@ set (wsjt_FSRCS
   lib/jt9_decode.f90
   lib/jt9fano.f90
   lib/ldpcsim144.f90
+  lib/fsk4hf/ldpcsim168.f90
   lib/ldpcsim40.f90
   lib/libration.f90
   lib/lorentzian.f90
@@ -1091,6 +1096,9 @@ target_link_libraries (ldpcsim40 wsjt_fort wsjt_cxx)
 add_executable (ldpcsim144 lib/ldpcsim144.f90 wsjtx.rc)
 target_link_libraries (ldpcsim144 wsjt_fort wsjt_cxx)
 
+add_executable (ldpcsim168 lib/fsk4hf/ldpcsim168.f90 lib/crc12.cpp wsjtx.rc)
+target_link_libraries (ldpcsim168 wsjt_fort wsjt_cxx)
+
 add_executable (msk144sim lib/msk144sim.f90 wsjtx.rc)
 target_link_libraries (msk144sim wsjt_fort wsjt_cxx)
 
diff --git a/lib/fsk4hf/bpdecode168.f90 b/lib/fsk4hf/bpdecode168.f90
new file mode 100644
index 000000000..d527f71d2
--- /dev/null
+++ b/lib/fsk4hf/bpdecode168.f90
@@ -0,0 +1,380 @@
+subroutine bpdecode168(llr,apmask,maxiterations,decoded,niterations)
+!
+! A log-domain belief propagation decoder for the (168,84) code.
+!
+integer, parameter:: N=168, K=84, M=N-K
+integer*1 codeword(N),cw(N),apmask(N)
+integer  colorder(N)
+integer*1 decoded(K)
+integer Nm(7,M)  ! 5, 6, or 7 bits per check 
+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,42,31,20,21,41,40,30,38,22,19,47,37,46,35,44,33,49,24, &
+   43,51,25,26,27,50,52,57,69,54,55,45,59,58,56,61,60,53,48,23,62, &
+   63,64,67,66,65,68,39,70,71,72,74,73,75,76,77,80,81,78,82,79,83, &
+   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/
+
+data Mn/               &
+  1,24,67, &
+  2,5,71, &
+  3,31,66, &
+  4,50,58, &
+  6,60,65, &
+  7,32,76, &
+  8,49,83, &
+  9,36,41, &
+  10,40,63, &
+  11,14,62, &
+  12,72,75, &
+  13,23,78, &
+  15,16,80, &
+  17,54,64, &
+  18,51,59, &
+  19,30,48, &
+  20,68,81, &
+  21,29,70, &
+  22,25,43, &
+  26,34,73, &
+  27,35,37, &
+  28,39,44, &
+  33,53,55, &
+  38,52,84, &
+  42,56,57, &
+  45,74,82, &
+  46,69,79, &
+  47,61,77, &
+  1,4,5, &
+  2,48,52, &
+  3,47,82, &
+  6,26,76, &
+  7,9,16, &
+  8,10,78, &
+  11,36,56, &
+  12,38,65, &
+  13,43,81, &
+  14,33,68, &
+  15,18,44, &
+  17,59,77, &
+  19,27,69, &
+  20,21,58, &
+  22,45,79, &
+  23,34,54, &
+  24,28,40, &
+  25,80,84, &
+  29,37,51, &
+  30,42,83, &
+  31,63,72, &
+  32,50,66, &
+  35,67,73, &
+  39,55,74, &
+  41,61,71, &
+  46,60,62, &
+  49,70,74, &
+  53,64,75, &
+  25,57,67, &
+  1,46,64, &
+  2,51,63, &
+  3,14,80, &
+  4,15,78, &
+  5,27,74, &
+  6,13,70, &
+  7,19,20, &
+  8,38,77, &
+  9,75,83, &
+  10,36,69, &
+  11,22,29, &
+  12,58,82, &
+  16,35,60, &
+  17,32,43, &
+  18,42,45, &
+  21,53,84, &
+  23,39,48, &
+  24,52,68, &
+  26,33,61, &
+  28,56,76, &
+  30,65,66, &
+  31,34,49, &
+  37,47,81, &
+  16,40,54, &
+  41,44,65, &
+  50,73,79, &
+  55,59,60, &
+  54,57,71, &
+  23,62,72, &
+  1,36,47, &
+  2,32,70, &
+  3,28,69, &
+  4,7,33, &
+  5,20,26, &
+  6,14,63, &
+  8,22,68, &
+  9,13,67, &
+  10,55,71, &
+  11,15,19, &
+  12,51,56, &
+  17,27,52, &
+  18,34,46, &
+  21,41,42, &
+  24,50,80, &
+  25,39,75, &
+  29,54,76, &
+  30,40,84, &
+  31,35,58, &
+  37,79,83, &
+  38,43,73, &
+  44,72,81, &
+  7,45,62, &
+  47,48,49, &
+  53,57,78, &
+  20,59,66, &
+  28,61,64, &
+  11,75,77, &
+  33,54,82, &
+  1,14,44, &
+  2,62,73, &
+  3,9,26, &
+  4,37,84, &
+  5,56,80, &
+  6,45,71, &
+  8,67,72, &
+  10,76,81, &
+  12,32,78, &
+  13,59,82, &
+  15,17,79, &
+  16,42,69, &
+  18,61,70, &
+  19,31,64, &
+  21,39,63, &
+  22,30,58, &
+  23,27,66, &
+  24,41,49, &
+  25,36,60, &
+  29,65,67, &
+  34,36,53, &
+  35,48,76, &
+  15,38,55, &
+  40,43,74, &
+  46,52,57, &
+  50,63,77, &
+  51,68,69, &
+  2,44,83, &
+  1,30,55, &
+  3,29,78, &
+  4,34,65, &
+  5,31,38, &
+  6,52,58, &
+  7,25,51, &
+  8,16,66, &
+  9,46,74, &
+  10,70,75, &
+  11,32,84, &
+  12,48,79, &
+  13,50,64, &
+  14,37,57, &
+  17,42,72, &
+  18,43,48, &
+  19,24,60, &
+  20,54,83, &
+  21,47,62, &
+  22,28,59, &
+  23,61,80, &
+  8,26,39, &
+  27,44,53, &
+  33,49,56, &
+  35,68,71, &
+  12,26,40/
+
+data Nm/               &
+   1,29,58,87,116,144,0,&
+   2,30,59,88,117,143,0,&
+   3,31,60,89,118,145,0,&
+   4,29,61,90,119,146,0,&
+   2,29,62,91,120,147,0,&
+   5,32,63,92,121,148,0,&
+   6,33,64,90,109,149,0,&
+   7,34,65,93,122,150,164,&
+   8,33,66,94,118,151,0,&
+   9,34,67,95,123,152,0,&
+   10,35,68,96,114,153,0,&
+   11,36,69,97,124,154,168,&
+   12,37,63,94,125,155,0,&
+   10,38,60,92,116,156,0,&
+   13,39,61,96,126,138,0,&
+   13,33,70,81,127,150,0,&
+   14,40,71,98,126,157,0,&
+   15,39,72,99,128,158,0,&
+   16,41,64,96,129,159,0,&
+   17,42,64,91,112,160,0,&
+   18,42,73,100,130,161,0,&
+   19,43,68,93,131,162,0,&
+   12,44,74,86,132,163,0,&
+   1,45,75,101,133,159,0,&
+   19,46,57,102,134,149,0,&
+   20,32,76,91,118,164,168,&
+   21,41,62,98,132,165,0,&
+   22,45,77,89,113,162,0,&
+   18,47,68,103,135,145,0,&
+   16,48,78,104,131,144,0,&
+   3,49,79,105,129,147,0,&
+   6,50,71,88,124,153,0,&
+   23,38,76,90,115,166,0,&
+   20,44,79,99,136,146,0,&
+   21,51,70,105,137,167,0,&
+   8,35,67,87,134,136,0,&
+   21,47,80,106,119,156,0,&
+   24,36,65,107,138,147,0,&
+   22,52,74,102,130,164,0,&
+   9,45,81,104,139,168,0,&
+   8,53,82,100,133,0,0,&
+   25,48,72,100,127,157,0,&
+   19,37,71,107,139,158,0,&
+   22,39,82,108,116,143,165,&
+   26,43,72,109,121,0,0,&
+   27,54,58,99,140,151,0,&
+   28,31,80,87,110,161,0,&
+   16,30,74,110,137,154,158,&
+   7,55,79,110,133,166,0,&
+   4,50,83,101,141,155,0,&
+   15,47,59,97,142,149,0,&
+   24,30,75,98,140,148,0,&
+   23,56,73,111,136,165,0,&
+   14,44,81,85,103,115,160,&
+   23,52,84,95,138,144,0,&
+   25,35,77,97,120,166,0,&
+   25,57,85,111,140,156,0,&
+   4,42,69,105,131,148,0,&
+   15,40,84,112,125,162,0,&
+   5,54,70,84,134,159,0,&
+   28,53,76,113,128,163,0,&
+   10,54,86,109,117,161,0,&
+   9,49,59,92,130,141,0,&
+   14,56,58,113,129,155,0,&
+   5,36,78,82,135,146,0,&
+   3,50,78,112,132,150,0,&
+   1,51,57,94,122,135,0,&
+   17,38,75,93,142,167,0,&
+   27,41,67,89,127,142,0,&
+   18,55,63,88,128,152,0,&
+   2,53,85,95,121,167,0,&
+   11,49,86,108,122,157,0,&
+   20,51,83,107,117,0,0,&
+   26,52,55,62,139,151,0,&
+   11,56,66,102,114,152,0,&
+   6,32,77,103,123,137,0,&
+   28,40,65,114,141,0,0,&
+   12,34,61,111,124,145,0,&
+   27,43,83,106,126,154,0,&
+   13,46,60,101,120,163,0,&
+   17,37,80,108,123,0,0,&
+   26,31,69,115,125,0,0,&
+   7,48,66,106,143,160,0,&
+   24,46,73,104,119,153,0/
+
+data nrw/   &
+6,6,6,6,6,6,6,7,6,6,6,7,6,6,6,6,6,6,6,6,6, &
+6,6,6,6,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,6, &
+6,7,5,6,6,7,6,6,6,6,6,7,6,6,6,6,6,6,6,6,6, &
+6,6,6,6,6,6,6,6,6,5,6,6,6,5,6,6,6,5,5,6,6/
+
+ncw=3
+
+toc=0
+tov=0
+tanhtoc=0
+!write(*,*) llr
+! 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
+    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(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
+niterations=-1
+return
+end subroutine bpdecode168
diff --git a/lib/fsk4hf/encode168.f90 b/lib/fsk4hf/encode168.f90
new file mode 100644
index 000000000..bf1212154
--- /dev/null
+++ b/lib/fsk4hf/encode168.f90
@@ -0,0 +1,141 @@
+subroutine encode168(message,codeword)
+! Encode an 84-bit message and return a 168-bit codeword. 
+! The generator matrix has dimensions (84,84). 
+! The code is a (168,84) 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 
+!
+character*21 g(84)
+integer*1 codeword(168)
+integer colorder(168)
+integer*1 gen(84,168)
+integer*1 itmp(168)
+integer*1 message(84)
+integer*1 pchecks(84)
+logical first
+data first/.true./
+data g/                  &  !parity generator matrix for (168,84) code
+ "25c5bf31ef6710fde9a5a", &
+ "18038ef7899cd97a77d96", &
+ "270dde504dad076c02b1f", &
+ "ed37fe12616565bd7d500", &
+ "12b99aa49b5367aff3838", &
+ "41cc27f2fac8b228aac21", &
+ "2265b233a3cff0b9cee24", &
+ "292760cd4f7f4a526a2f1", &
+ "2b3db4c8bd831911680cc", &
+ "cef2b24ce203bdc60b266", &
+ "5045a24f9340915d807ab", &
+ "3592b7fc60ba85139502e", &
+ "9318023145637bd798f0e", &
+ "ad796023c3d58d1e6509c", &
+ "3da5eab57f040e75d7413", &
+ "27466d1d2734d0ff64830", &
+ "2ed50bb1ce313bbfb1ab0", &
+ "9a616bda01b25b7e6eeaf", &
+ "a84c8c1e9df103169d10d", &
+ "a40da29b4aca9234a8942", &
+ "dd258d02d79a5f209d3d0", &
+ "bdfdc06713511997b5621", &
+ "25c58f12f4096cd8ead1a", &
+ "b2638a478f21e10fe97de", &
+ "4051020f43c605d458156", &
+ "f651aad14322a526dae35", &
+ "a1c147e31bcc9d87330bf", &
+ "7524b53d996d48284647b", &
+ "a72e7d25ce31b27282e56", &
+ "a97f53b019022350b7519", &
+ "56106c6340c0810790984", &
+ "c63b8e03a57208635992b", &
+ "43a3de2aa3a2b1afb65dc", &
+ "9baa64847ead03b77fecc", &
+ "251cbd1895c8839c46b0d", &
+ "2858107dde2d173e13530", &
+ "20096f6a870f636b704e7", &
+ "7f833ccbceec52dd6eb79", &
+ "a9108dd77b8015b75242a", &
+ "689666a79e5579c916236", &
+ "aa5dff46459787f69911f", &
+ "794558c13138d08171089", &
+ "c937042857b291cee8dfd", &
+ "6f0bf3248bb9a231366b8", &
+ "1c09e756ef1656c96f2d2", &
+ "073b875b6774e71fba549", &
+ "f7d840aafc037febd2d5c", &
+ "dcc0e7d0da5fe17c99ad3", &
+ "98238ef7819cd97a77d94", &
+ "177c2594743477421a262", &
+ "7d01a833c19374fbaaa6e", &
+ "7bb800216660482ffd1c4", &
+ "39a92e2dba0d4cfda98d2", &
+ "44b8d88622698816456a8", &
+ "791db2334d6d86639229b", &
+ "ba6004b086bd38559ea48", &
+ "f94558e13138d18170089", &
+ "08ba145302cfbed7845ae", &
+ "fb8e64b6da3602168ed38", &
+ "1045a2cf1340915d8072b", &
+ "7592b6fc64ba85139582e", &
+ "3eb238a11bc6654452bae", &
+ "b69d8d23b1ea170f70214", &
+ "0123dfae84fb20462a614", &
+ "4131066ad52a339b3c0d7", &
+ "fd2cc26850951c43ed737", &
+ "a644d4eb7e56c40f0d050", &
+ "0c3bd9d5dab7c9ee2c8fc", &
+ "4a198b37af56d7ceffb56", &
+ "b6e946c429294cf0eed8b", &
+ "98384d75e758774f5ff3b", &
+ "5c58e5d9a4d0531d37384", &
+ "7a0af02719afed521fd06", &
+ "8cd5b2e694e7854abbc70", &
+ "1a2f061912d0ea19702d3", &
+ "6ffbce557d8fa691a50e8", &
+ "d43438e2e2ed5d9f14011", &
+ "8d502106083b809adba00", &
+ "67e22f9b9983aa715964d", &
+ "b31f3a3f3c1f406b1fd58", &
+ "529f60ac291f827d97331", &
+ "476a815424f2e2cbe641f", &
+ "81c82c89bcc3feec42458", &
+ "2c882d0e281b178e80364"/
+
+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,42,31,20,21,41,40,30,38,22,19,47,37,46,35,44,33,49,24, &
+  43,51,25,26,27,50,52,57,69,54,55,45,59,58,56,61,60,53,48,23,62, &
+  63,64,67,66,65,68,39,70,71,72,74,73,75,76,77,80,81,78,82,79,83, &
+  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/
+
+save first,gen
+
+if( first ) then ! fill the generator matrix
+  gen=0
+  do i=1,84
+    do j=1,21
+      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, 84 
+  nsum=0
+  do j=1, 84 
+    nsum=nsum+message(j)*gen(i,j)
+  enddo
+  pchecks(i)=mod(nsum,2)
+enddo
+itmp(1:84)=pchecks
+itmp(85:168)=message(1:84)
+codeword(colorder+1)=itmp(1:168)
+
+return
+end subroutine encode168
diff --git a/lib/fsk4hf/extractmessage168.f90 b/lib/fsk4hf/extractmessage168.f90
new file mode 100644
index 000000000..466d2a638
--- /dev/null
+++ b/lib/fsk4hf/extractmessage168.f90
@@ -0,0 +1,48 @@
+subroutine extractmessage168(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
+  use iso_c_binding, only: c_loc,c_size_t
+  use crc
+  use packjt
+
+  character*22 msgreceived
+  character*12 call1,call2
+  character*12 recent_calls(nrecent)
+  integer*1 decoded(84)
+  integer*1, target::  i1Dec8BitBytes(11)
+  integer*4 i4Dec6BitWords(12)
+
+! Collapse 84 decoded bits to 11 bytes. Bytes 1-9 are the message, byte 10 and first half of byte 11 is the crc
+  do ibyte=1,9
+      itmp=0
+    do ibit=1,8
+      itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
+    enddo
+    i1Dec8BitBytes(ibyte)=itmp
+  enddo
+! Need to pack the crc into bytes 10 and 11 for crc12_check
+  i1Dec8BitBytes(10)=decoded(73)*8+decoded(74)*4+decoded(75)*2+decoded(76)
+  i1Dec8BitBytes(11)=decoded(77)*128+decoded(78)*64+decoded(79)*2*32+decoded(80)*16
+  i1Dec8BitBytes(11)=i1Dec8BitBytes(11)+decoded(81)*8+decoded(82)*4+decoded(83)*2+decoded(84)
+
+  if( crc12_check(c_loc (i1Dec8BitBytes), 11) ) then
+! CRC12 checks out --- unpack 72-bit message
+    do ibyte=1,12
+      itmp=0
+      do ibit=1,6
+        itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
+      enddo
+      i4Dec6BitWords(ibyte)=itmp
+    enddo
+    call unpackmsg144(i4Dec6BitWords,msgreceived,call1,call2)
+    ncrcflag=1
+    if( call1(1:2) .ne. 'CQ' .and. call1(1:2) .ne. '  ' ) then
+      call update_recent_calls(call1,recent_calls,nrecent)
+    endif
+    if( call2(1:2) .ne. '  ' ) then
+      call update_recent_calls(call2,recent_calls,nrecent)
+    endif
+  else
+    msgreceived=' '
+    ncrcflag=-1
+  endif 
+  return
+  end subroutine extractmessage168
diff --git a/lib/fsk4hf/ldpcsim168.f90 b/lib/fsk4hf/ldpcsim168.f90
new file mode 100644
index 000000000..8d0219a41
--- /dev/null
+++ b/lib/fsk4hf/ldpcsim168.f90
@@ -0,0 +1,219 @@
+program ldpcsim168
+! End to end test of the (168,84)/crc12 encoder and decoder.
+use crc
+use packjt
+
+parameter(NRECENT=10)
+character*12 recent_calls(NRECENT)
+character*22 msg,msgsent,msgreceived
+character*8 arg
+integer*1, allocatable ::  codeword(:), decoded(:), message(:)
+integer*1, target:: i1Msg8BitBytes(11)
+integer*1 msgbits(84)
+integer*1 apmask(168)
+integer*2 checksum
+integer*4 i4Msg6BitWords(13)
+integer colorder(168)
+integer nerrtot(168),nerrdec(168)
+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,42,31,20,21,41,40,30,38,22,19,47,37,46,35,44,33,49,24, &
+  43,51,25,26,27,50,52,57,69,54,55,45,59,58,56,61,60,53,48,23,62, &
+  63,64,67,66,65,68,39,70,71,72,74,73,75,76,77,80,81,78,82,79,83, &
+  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/
+
+do i=1,NRECENT
+  recent_calls(i)='            '
+enddo
+nerrtot=0
+nerrdec=0
+
+nargs=iargc()
+if(nargs.ne.3) then
+   print*,'Usage: ldpcsim  niter  #trials  s '
+   print*,'eg:    ldpcsim    10   1000    0.75'
+   return
+endif
+call getarg(1,arg)
+read(arg,*) max_iterations 
+call getarg(2,arg)
+read(arg,*) ntrials 
+call getarg(3,arg)
+read(arg,*) s
+
+fsk=.true.
+bpsk=.false.
+
+! don't count crc bits as data bits
+N=168
+K=84
+! scale Eb/No for a (168,72) code
+rate=real(72)/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"
+  call packmsg(msg,i4Msg6BitWords,itype)  !Pack into 12 6-bit bytes
+  call unpackmsg(i4Msg6BitWords,msgsent)  !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 = crc12 (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 = crc12_check(c_loc (i1Msg8BitBytes), 11)
+  if( checksumok ) write(*,*) 'Good checksum'
+
+  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
+  i1=i1Msg8BitBytes(10) ! First 4 bits of crc12 are LSB of this byte
+  do ibit=1,4
+    msgbits(72+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(76+ibit)=iand(1,ishft(i1,ibit-8))
+  enddo
+
+  write(*,*) 'message'
+  write(*,'(11(8i1,1x))') msgbits
+
+  call encode168(msgbits,codeword)
+  call init_random_seed()
+  call sgran()
+
+  write(*,*) 'codeword' 
+  write(*,'(21(8i1,1x))') codeword
+
+write(*,*) "Es/N0  SNR2500   ngood  nundetected nbadcrc   sigma"
+do idb = -10, 24 
+  db=idb/2.0-1.0
+!  sigma=1/sqrt( 2*rate*(10**(db/10.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
+    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(168-84+1:168-84+nap)+1)=5*(2.0*msgbits(1:nap)-1.0)
+    apmask=0
+    apmask(colorder(168-84+1:168-84+nap)+1)=1
+
+! max_iterations is max number of belief propagation iterations
+    call bpdecode168(llr, apmask, max_iterations, decoded, niterations)
+! If the decoder finds a valid codeword, niterations will be .ge. 0.
+    if( niterations .ge. 0 ) then
+      call extractmessage168(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
+      if( ncrcflag .eq. 1 ) then
+        ncrcflag=1
+      else 
+        ncrcflag=0
+      endif
+      if( ncrcflag .ne. 1 ) then
+        nbadcrc=nbadcrc+1
+      endif
+      nueflag=0
+
+! Check the message plus crc against what was sent.
+      do i=1,K
+        if( msgbits(i) .ne. decoded(i) ) then
+          nueflag=1
+        endif
+      enddo
+      if( ncrcflag .eq. 1 .and. nueflag .eq. 0 ) then
+        ngood=ngood+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(8.333/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,128
+  write(23,'(i4,2x,i10,i10,f10.2)') i,nerrdec(i),nerrtot(i),real(nerrdec(i))/real(nerrtot(i)+1e-10)
+enddo
+close(23)
+
+end program ldpcsim168