From 65fa6aed6cd33bc51e7c64744f3053cfa2ddb3ee Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Mon, 13 Apr 2020 10:23:02 -0500 Subject: [PATCH] Add a new (174,74) LDPC code (regular, column weight 3) and supporting files. --- lib/fsk4hf/bpdecode174_74.f90 | 113 ++++++++ lib/fsk4hf/encode174_74.f90 | 47 ++++ lib/fsk4hf/get_crc24.f90 | 25 ++ lib/fsk4hf/ldpc_174_74_generator.f90 | 105 +++++++ lib/fsk4hf/ldpc_174_74_parity.f90 | 288 +++++++++++++++++++ lib/fsk4hf/ldpcsim174_74.f90 | 144 ++++++++++ lib/fsk4hf/osd174_74.f90 | 405 +++++++++++++++++++++++++++ 7 files changed, 1127 insertions(+) create mode 100644 lib/fsk4hf/bpdecode174_74.f90 create mode 100644 lib/fsk4hf/encode174_74.f90 create mode 100644 lib/fsk4hf/get_crc24.f90 create mode 100644 lib/fsk4hf/ldpc_174_74_generator.f90 create mode 100644 lib/fsk4hf/ldpc_174_74_parity.f90 create mode 100644 lib/fsk4hf/ldpcsim174_74.f90 create mode 100644 lib/fsk4hf/osd174_74.f90 diff --git a/lib/fsk4hf/bpdecode174_74.f90 b/lib/fsk4hf/bpdecode174_74.f90 new file mode 100644 index 000000000..66ab28141 --- /dev/null +++ b/lib/fsk4hf/bpdecode174_74.f90 @@ -0,0 +1,113 @@ +subroutine bpdecode174_74(llr,apmask,maxiterations,message50,cw,nharderror,iter) +! +! A log-domain belief propagation decoder for the (174,74) code. +! + + integer, parameter:: N=174, K=74, M=N-K + integer*1 cw(N),apmask(N) + integer*1 decoded(K) + integer*1 message50(50) + integer nrw(M),ncw + integer Nm(6,M) + 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 + + include "ldpc_174_74_parity.f90" + + 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 + nclast=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 + if( ncheck .eq. 0 ) then ! we have a codeword - if crc is good, return it + decoded=cw(1:74) + call get_crc24(decoded,nbadcrc) + nharderror=count( (2*cw-1)*llr .lt. 0.0 ) + if(nbadcrc.eq.0) then + message50=decoded(1:50) + return + endif + 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_74 diff --git a/lib/fsk4hf/encode174_74.f90 b/lib/fsk4hf/encode174_74.f90 new file mode 100644 index 000000000..ba8ab18b1 --- /dev/null +++ b/lib/fsk4hf/encode174_74.f90 @@ -0,0 +1,47 @@ +subroutine encode174_74(message,codeword) + use, intrinsic :: iso_c_binding + use iso_c_binding, only: c_loc,c_size_t + use crc + + integer, parameter:: N=174, K=74, M=N-K + character*24 c24 + integer*1 codeword(N) + integer*1 gen(M,K) + integer*1 message(K) + integer*1 pchecks(M) + integer*1, target :: i1MsgBytes(10) + integer*4 ncrc24 + include "ldpc_174_74_generator.f90" + logical first + data first/.true./ + save first,gen + + if( first ) then ! fill the generator matrix + gen=0 + do i=1,M + do j=1,19 + read(g(i)(j:j),"(Z1)") istr + ibmax=4 + if(j.eq.19) ibmax=2 + 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 + + codeword(1:K)=message + codeword(K+1:N)=pchecks + + return +end subroutine encode174_74 diff --git a/lib/fsk4hf/get_crc24.f90 b/lib/fsk4hf/get_crc24.f90 new file mode 100644 index 000000000..d73310792 --- /dev/null +++ b/lib/fsk4hf/get_crc24.f90 @@ -0,0 +1,25 @@ +subroutine get_crc24(mc,ncrc) +! +! 1. To calculate 24-bit CRC, mc(1:50) is the message and mc(51:74) are zero. +! 2. To check a received CRC, mc(1:74) is the received message plus CRC. +! ncrc will be zero if the received message/CRC are consistent. +! + character c24*24 + + integer*1 mc(74),r(25),p(25) + integer ncrc +! polynomial for 24-bit CRC 0x100065b + data p/1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,1,1,0,1,1/ + +! divide by polynomial + r=mc(1:25) + do i=0,49 + r(25)=mc(i+25) + r=mod(r+r(1)*p,2) + r=cshift(r,1) + enddo + + write(c24,'(24b1)') r(1:24) + read(c24,'(b24.24)') ncrc + +end subroutine get_crc24 diff --git a/lib/fsk4hf/ldpc_174_74_generator.f90 b/lib/fsk4hf/ldpc_174_74_generator.f90 new file mode 100644 index 000000000..116ccfddc --- /dev/null +++ b/lib/fsk4hf/ldpc_174_74_generator.f90 @@ -0,0 +1,105 @@ +! generator matrix for regular column weight 3 (174,74) LDPC code. + +character*19 g(100) + +data g/ & + "b190e319bd45882ed74", & + "b159282d395467cabe4", & + "f502387ec63db738358", & + "a6c7911729277b2a178", & + "05d812d04122ff2842c", & + "eb040701b66b26d12ec", & + "c617358e34398e73c6c", & + "37b62b499bbb84aeb0c", & + "60257b5d4e41594a250", & + "e8ac26253c0268ba33c", & + "0f243baf67353230318", & + "521d5eb1268bed86854", & + "53bab7dbe89962bba00", & + "417abaf10e0912604b8", & + "c0d371dbb301f49aae8", & + "e7014cf533a2cb9fd24", & + "948175882e16ecd8cc8", & + "87db37137999fb15504", & + "2557139852451e678c8", & + "aaaf0b6b1f70db8e5ac", & + "f5be069b0a41fd5bb28", & + "e7789f2237b2175d494", & + "94554737d22b00d5980", & + "525e935db67c1af214c", & + "9c57c640427a2c2e33c", & + "9a82e00fb570e371cac", & + "39ebbdd43570f690818", & + "a037514614e0d5cc2a4", & + "ff19fc0eee4376f6de0", & + "f8853aad262b1a14cf0", & + "f5687424fe7c5156ee0", & + "fba0aa4876b79e45d78", & + "dfdfb60046769dec900", & + "600b4517a14560fad64", & + "39c618d3f629809c064", & + "5c821087e8c365869f8", & + "a4f26e15e3ef8264c04", & + "ac230e4147016f5bf98", & + "e11a6981f5257957d84", & + "b9dd003c09cf2abc5b0", & + "326ff2588a1bfa6a310", & + "a84e8e04722185f23ac", & + "8a66abe81aff313f9a8", & + "f6047ea2cad01957e08", & + "f14b63fdff262eb74bc", & + "7588be7336de21f7680", & + "312d0e1d5d1c3666fec", & + "5ab69333712cdbf9c38", & + "3c8e8c949be183939f4", & + "ed3b36d068e55ef76d8", & + "8193b051800415c06e0", & + "bc8e88949be18393bd4", & + "a7db37037999fb15404", & + "c5ec69ecc57ed7800b8", & + "475b645148268e10afc", & + "1fe90fea7ae941c04a8", & + "513b196d2a6e43c9504", & + "ffc27ceba420d04f468", & + "972c2cc31e578dba968", & + "7fc874b734a8188a2f8", & + "3d0327a801275734cc4", & + "b1e77d50857f56b6a40", & + "f25389644e47dae2384", & + "4e7e815e6b3c20507a0", & + "27d63d2e80a23f057cc", & + "381388a8a6fad77ec50", & + "b785abf747ea18bc350", & + "40a2a8214e2bed48090", & + "0e891f175b06fed80d4", & + "dbd155acd9fbec5b4a4", & + "9d3476e615c702f8e60", & + "050ea06fbd1f532d164", & + "d03767bca8394f31628", & + "455d568ff3047e9d5ac", & + "6b343bcf7378e1283f4", & + "a0d371dbb311f49aae8", & + "36ea237c911eb2ac27c", & + "54a636ec612a744f368", & + "5cabf5c9d5a0d2d9ba4", & + "00d632bffc3dac0d548", & + "d86bf5593c70dcb91fc", & + "bada10bb78be8c219c0", & + "b98028f926fed2beab0", & + "c0347b3cc45c2888094", & + "0662d6a3c2974e0a910", & + "8036b9e83c9fdc2cda8", & + "e5db38aad196024c21c", & + "746c8af5783b5daedcc", & + "1dc47211c27e39ec5dc", & + "6b98898e40559a2e128", & + "52d9077dbfa44c6d75c", & + "9ca1e6bd4515559a054", & + "7b2dd815e5991f88d14", & + "bfde5ebc6e09940460c", & + "487f5ffeaf139c209f4", & + "08d6b3c9686cc0f6ff4", & + "e198f5466141f53ab84", & + "0a7c7af0ac612d14f40", & + "a4192113ec53f4d165c", & + "1423ae72e003614be88"/ diff --git a/lib/fsk4hf/ldpc_174_74_parity.f90 b/lib/fsk4hf/ldpc_174_74_parity.f90 new file mode 100644 index 000000000..12187cc32 --- /dev/null +++ b/lib/fsk4hf/ldpc_174_74_parity.f90 @@ -0,0 +1,288 @@ +! parity check matrix for regular column weight 3 (174,74) LDPC code + +data Mn/ & + 28, 32, 98, & + 1, 94, 95, & + 70, 71, 94, & + 3, 9, 39, & + 4, 22, 84, & + 5, 25, 85, & + 6, 55, 100, & + 7, 41, 67, & + 8, 62, 77, & + 10, 37, 40, & + 11, 16, 36, & + 12, 29, 47, & + 14, 57, 91, & + 15, 49, 59, & + 17, 18, 52, & + 19, 48, 58, & + 20, 34, 72, & + 21, 38, 87, & + 23, 46, 79, & + 24, 43, 83, & + 26, 74, 78, & + 27, 51, 98, & + 30, 35, 42, & + 31, 63, 88, & + 25, 33, 82, & + 44, 53, 96, & + 50, 75, 90, & + 46, 54, 71, & + 56, 76, 81, & + 60, 65, 69, & + 61, 95, 97, & + 64, 89, 99, & + 66, 73, 76, & + 68, 80, 92, & + 74, 86, 87, & + 52, 55, 93, & + 1, 3, 21, & + 2, 57, 64, & + 4, 11, 82, & + 5, 60, 81, & + 6, 13, 66, & + 7, 43, 59, & + 8, 27, 85, & + 9, 34, 94, & + 10, 28, 88, & + 12, 19, 53, & + 14, 33, 65, & + 15, 75, 84, & + 16, 56, 68, & + 17, 44, 90, & + 18, 23, 73, & + 20, 26, 83, & + 22, 42, 91, & + 24, 70, 79, & + 29, 96, 97, & + 30, 48, 77, & + 31, 37, 67, & + 32, 35, 78, & + 36, 80, 89, & + 38, 62, 93, & + 39, 54, 72, & + 40, 58, 61, & + 41, 51, 86, & + 9, 45, 63, & + 46, 47, 92, & + 49, 50, 80, & + 69, 75, 98, & + 14, 71, 99, & + 21, 85, 100, & + 1, 6, 50, & + 2, 74, 98, & + 3, 73, 92, & + 4, 32, 79, & + 5, 63, 96, & + 1, 29, 79, & + 2, 6, 84, & + 3, 36, 78, & + 4, 59, 69, & + 5, 72, 80, & + 2, 41, 44, & + 7, 38, 90, & + 8, 58, 99, & + 9, 42, 49, & + 10, 46, 74, & + 11, 14, 73, & + 12, 85, 88, & + 13, 25, 93, & + 5, 11, 61, & + 15, 16, 94, & + 8, 15, 32, & + 17, 64, 76, & + 18, 60, 71, & + 19, 34, 57, & + 20, 77, 96, & + 21, 31, 82, & + 22, 23, 48, & + 16, 17, 22, & + 24, 35, 86, & + 13, 51, 64, & + 26, 33, 39, & + 4, 27, 47, & + 28, 54, 56, & + 1, 62, 76, & + 30, 53, 100, & + 21, 46, 57, & + 32, 55, 61, & + 33, 92, 95, & + 19, 67, 68, & + 14, 35, 75, & + 3, 84, 99, & + 37, 75, 92, & + 7, 30, 63, & + 23, 39, 50, & + 27, 40, 44, & + 41, 66, 91, & + 9, 10, 69, & + 12, 43, 52, & + 6, 20, 88, & + 43, 45, 89, & + 24, 31, 60, & + 47, 50, 70, & + 26, 38, 45, & + 25, 42, 98, & + 40, 79, 87, & + 51, 68, 83, & + 52, 97, 98, & + 48, 71, 78, & + 28, 29, 80, & + 53, 59, 73, & + 49, 54, 95, & + 34, 36, 90, & + 55, 82, 89, & + 4, 77, 90, & + 49, 58, 81, & + 18, 56, 77, & + 62, 63, 65, & + 60, 91, 100, & + 2, 11, 62, & + 15, 24, 93, & + 37, 66, 83, & + 65, 85, 86, & + 5, 10, 48, & + 1, 69, 89, & + 67, 81, 87, & + 13, 56, 75, & + 70, 72, 97, & + 8, 50, 57, & + 6, 19, 33, & + 7, 25, 99, & + 74, 94, 96, & + 14, 52, 87, & + 16, 30, 40, & + 20, 42, 79, & + 17, 21, 72, & + 9, 12, 41, & + 18, 61, 67, & + 3, 83, 97, & + 26, 80, 91, & + 23, 55, 65, & + 27, 31, 95, & + 28, 84, 86, & + 29, 34, 93, & + 46, 51, 63, & + 35, 39, 76, & + 13, 44, 78, & + 32, 37, 38, & + 22, 43, 92, & + 45, 53, 54, & + 58, 73, 74, & + 47, 64, 100, & + 59, 68, 85, & + 66, 82, 96, & + 36, 81, 88, & + 2, 45, 70/ + +data Nm/ & + 2, 37, 70, 75, 103, 143, & + 38, 71, 76, 80, 138, 174, & + 4, 37, 72, 77, 110, 157, & + 5, 39, 73, 78, 101, 133, & + 6, 40, 74, 79, 88, 142, & + 7, 41, 70, 76, 118, 148, & + 8, 42, 81, 112, 149, 0, & + 9, 43, 82, 90, 147, 0, & + 4, 44, 64, 83, 116, 155, & + 10, 45, 84, 116, 142, 0, & + 11, 39, 85, 88, 138, 0, & + 12, 46, 86, 117, 155, 0, & + 41, 87, 99, 145, 165, 0, & + 13, 47, 68, 85, 109, 151, & + 14, 48, 89, 90, 139, 0, & + 11, 49, 89, 97, 152, 0, & + 15, 50, 91, 97, 154, 0, & + 15, 51, 92, 135, 156, 0, & + 16, 46, 93, 108, 148, 0, & + 17, 52, 94, 118, 153, 0, & + 18, 37, 69, 95, 105, 154, & + 5, 53, 96, 97, 167, 0, & + 19, 51, 96, 113, 159, 0, & + 20, 54, 98, 120, 139, 0, & + 6, 25, 87, 123, 149, 0, & + 21, 52, 100, 122, 158, 0, & + 22, 43, 101, 114, 160, 0, & + 1, 45, 102, 128, 161, 0, & + 12, 55, 75, 128, 162, 0, & + 23, 56, 104, 112, 152, 0, & + 24, 57, 95, 120, 160, 0, & + 1, 58, 73, 90, 106, 166, & + 25, 47, 100, 107, 148, 0, & + 17, 44, 93, 131, 162, 0, & + 23, 58, 98, 109, 164, 0, & + 11, 59, 77, 131, 173, 0, & + 10, 57, 111, 140, 166, 0, & + 18, 60, 81, 122, 166, 0, & + 4, 61, 100, 113, 164, 0, & + 10, 62, 114, 124, 152, 0, & + 8, 63, 80, 115, 155, 0, & + 23, 53, 83, 123, 153, 0, & + 20, 42, 117, 119, 167, 0, & + 26, 50, 80, 114, 165, 0, & + 64, 119, 122, 168, 174, 0, & + 19, 28, 65, 84, 105, 163, & + 12, 65, 101, 121, 170, 0, & + 16, 56, 96, 127, 142, 0, & + 14, 66, 83, 130, 134, 0, & + 27, 66, 70, 113, 121, 147, & + 22, 63, 99, 125, 163, 0, & + 15, 36, 117, 126, 151, 0, & + 26, 46, 104, 129, 168, 0, & + 28, 61, 102, 130, 168, 0, & + 7, 36, 106, 132, 159, 0, & + 29, 49, 102, 135, 145, 0, & + 13, 38, 93, 105, 147, 0, & + 16, 62, 82, 134, 169, 0, & + 14, 42, 78, 129, 171, 0, & + 30, 40, 92, 120, 137, 0, & + 31, 62, 88, 106, 156, 0, & + 9, 60, 103, 136, 138, 0, & + 24, 64, 74, 112, 136, 163, & + 32, 38, 91, 99, 170, 0, & + 30, 47, 136, 141, 159, 0, & + 33, 41, 115, 140, 172, 0, & + 8, 57, 108, 144, 156, 0, & + 34, 49, 108, 125, 171, 0, & + 30, 67, 78, 116, 143, 0, & + 3, 54, 121, 146, 174, 0, & + 3, 28, 68, 92, 127, 0, & + 17, 61, 79, 146, 154, 0, & + 33, 51, 72, 85, 129, 169, & + 21, 35, 71, 84, 150, 169, & + 27, 48, 67, 109, 111, 145, & + 29, 33, 91, 103, 164, 0, & + 9, 56, 94, 133, 135, 0, & + 21, 58, 77, 127, 165, 0, & + 19, 54, 73, 75, 124, 153, & + 34, 59, 66, 79, 128, 158, & + 29, 40, 134, 144, 173, 0, & + 25, 39, 95, 132, 172, 0, & + 20, 52, 125, 140, 157, 0, & + 5, 48, 76, 110, 161, 0, & + 6, 43, 69, 86, 141, 171, & + 35, 63, 98, 141, 161, 0, & + 18, 35, 124, 144, 151, 0, & + 24, 45, 86, 118, 173, 0, & + 32, 59, 119, 132, 143, 0, & + 27, 50, 81, 131, 133, 0, & + 13, 53, 115, 137, 158, 0, & + 34, 65, 72, 107, 111, 167, & + 36, 60, 87, 139, 162, 0, & + 2, 3, 44, 89, 150, 0, & + 2, 31, 107, 130, 160, 0, & + 26, 55, 74, 94, 150, 172, & + 31, 55, 126, 146, 157, 0, & + 1, 22, 67, 71, 123, 126, & + 32, 68, 82, 110, 149, 0, & + 7, 69, 104, 137, 170, 0/ + +data nrw/ & + 6,6,6,6,6,6,5,5,6,5,5,5,5,6,5,5,5,5,5,5, & + 6,5,5,5,5,5,5,5,5,5,5,6,5,5,5,5,5,5,5,5, & + 5,5,5,5,5,6,5,5,5,6,5,5,5,5,5,5,5,5,5,5, & + 5,5,6,5,5,5,5,5,5,5,5,5,6,6,6,5,5,5,6,6, & + 5,5,5,5,6,5,5,5,5,5,5,6,5,5,5,6,5,6,5,5/ + +ncw=3 diff --git a/lib/fsk4hf/ldpcsim174_74.f90 b/lib/fsk4hf/ldpcsim174_74.f90 new file mode 100644 index 000000000..51ca6eb84 --- /dev/null +++ b/lib/fsk4hf/ldpcsim174_74.f90 @@ -0,0 +1,144 @@ +program ldpcsim174_74 + +! End-to-end test of the (174,74)/crc24 encoder and decoders. + + use crc + use packjt + + parameter(N=174, K=74, M=N-K) + character*8 arg + character*50 cmsg + character*24 c24 + integer*1 msgbits(74) + integer*1 apmask(174) + integer*1 cw(174) + integer*1 codeword(N),message(50) + integer ncrc24 + integer nerrtot(174),nerrdec(174),nmpcbad(74) + real rxdata(N),llr(N) + real dllr(174),llrd(174) + + data cmsg/'11111111000000001111111100000000111111110000000011'/ + + 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.5) then + print*,'Usage: ldpcsim niter ndeep #trials s K' + print*,'e.g. ldpcsim174_74 20 5 1000 0.85 64' + print*,'s : if negative, then value is ignored and sigma is calculated from SNR.' + print*,'niter: is the number of BP iterations.' + print*,'ndeep: -1 is BP only, ndeep>=0 is OSD order' + print*,'K :is the number of message+CRC bits and must be in the range [50,74]' + 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 + call getarg(5,arg) + read(arg,*) Keff + + rate=real(Keff)/real(N) + + write(*,*) "code rate: ",rate + write(*,*) "niter : ",max_iterations + write(*,*) "ndeep : ",ndeep + write(*,*) "s : ",s + write(*,*) "K : ",Keff + + msgbits=0 + read(cmsg,'(50i1)') msgbits(1:50) + write(*,*) 'message' + write(*,'(74i1)') msgbits + + call get_crc24(msgbits,ncrc24) + write(c24,'(b24.24)') ncrc24 + read(c24,'(24i1)') msgbits(51:74) + call encode174_74(msgbits,codeword) + call init_random_seed() + call sgran() + + write(*,*) 'codeword' + write(*,'(50i1,1x,24i1,1x,100i1)') codeword + + write(*,*) "Eb/N0 Es/N0 ngood nundetected sigma symbol error rate" + do idb = 8,-3,-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 + nberr=0 + do itrial=1, ntrials +! Create a realization of a noisy received word + do i=1,N + rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran() + 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 + + 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 bpdecode174_74(llr,apmask,max_iterations,message,cw,nharderror,niterations) + dmin=0.0 + if( (nharderror .lt. 0) .and. (ndeep .ge. 0) ) then + call osd174_74(llr, Keff, apmask, ndeep, message, cw, nharderror, dmin) + endif + + if(nharderror.ge.0) then + n2err=0 + do i=1,N + if( cw(i)*(2*codeword(i)-1.0) .lt. 0 ) n2err=n2err+1 + enddo + if(n2err.eq.0) then + ngood=ngood+1 + else + nue=nue+1 + endif + endif + enddo +! snr2500=db+10*log10(200.0/116.0/2500.0) + esn0=db+10*log10(rate) + pberr=real(nberr)/(real(ntrials*N)) + write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,8x,f5.2,8x,e10.3)") db,esn0,ngood,nue,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 ldpcsim174_74 diff --git a/lib/fsk4hf/osd174_74.f90 b/lib/fsk4hf/osd174_74.f90 new file mode 100644 index 000000000..7f87b50bf --- /dev/null +++ b/lib/fsk4hf/osd174_74.f90 @@ -0,0 +1,405 @@ +subroutine osd174_74(llr,k,apmask,ndeep,message74,cw,nhardmin,dmin) +! +! An ordered-statistics decoder for the (174,74) code. +! Message payload is 50 bits. Any or all of a 24-bit CRC can be +! used for detecting incorrect codewords. The remaining CRC bits are +! cascaded with the LDPC code for the purpose of improving the +! distance spectrum of the code. +! +! If p1 (0.le.p1.le.24) is the number of CRC24 bits that are +! to be used for bad codeword detection, then the argument k should +! be set to 50+p1. +! +! Valid values for k are in the range [50,74]. +! + character*24 c24 + integer, parameter:: N=174 + integer*1 apmask(N),apmaskr(N) + integer*1, allocatable, save :: gen(:,:) + integer*1, allocatable :: genmrb(:,:),g2(:,:) + integer*1, allocatable :: temp(:),m0(:),me(:),mi(:),misub(:),e2sub(:),e2(:),ui(:) + integer*1, allocatable :: r2pat(:) + integer indices(N),nxor(N) + integer*1 cw(N),ce(N),c0(N),hdec(N) + integer*1, allocatable :: decoded(:) + integer*1 message50(50),message74(74) + integer indx(N) + real llr(N),rx(N),absrx(N) + +!include "ldpc_174_74_generator.f90" + + logical first,reset + data first/.true./ + save first + + allocate( genmrb(k,N), g2(N,k) ) + allocate( temp(k), m0(k), me(k), mi(k), misub(k), e2sub(N-k), e2(N-k), ui(N-k) ) + allocate( r2pat(N-k), decoded(k) ) + + if( first ) then ! fill the generator matrix +! +! Create generator matrix for partial CRC cascaded with LDPC code. +! +! Let p2=74-k and p1+p2=24. +! +! The last p2 bits of the CRC24 are cascaded with the LDPC code. +! +! The first p1=k-50 CRC24 bits will be used for error detection. +! + allocate( gen(k,N) ) + gen=0 + do i=1,k + message74=0 + message74(i)=1 + if(i.le.50) then + call get_crc24(message74,ncrc24) + write(c24,'(b24.24)') ncrc24 + read(c24,'(24i1)') message74(51:74) + message74(51:k)=0 + endif + call encode174_74(message74,cw) + gen(i,:)=cw + enddo + + first=.false. + endif + + rx=llr + apmaskr=apmask + +! 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 mrbencode74(m0,c0,g2,N,k) + nxor=ieor(c0,hdec) + nhardmin=sum(nxor) + dmin=sum(nxor*absrx) + + cw=c0 + ntotal=0 + nrejected=0 + npre1=0 + npre2=0 + + if(ndeep.eq.0) goto 998 ! norder=0 + if(ndeep.gt.6) ndeep=6 + 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=1 + nt=40 + ntheta=12 + ntau=19 + elseif(ndeep.eq.5) then + nord=3 + npre1=1 + npre2=1 + nt=40 + ntheta=12 + ntau=19 + elseif(ndeep.eq.6) then + nord=4 + 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 + d1=0. + 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 mrbencode74(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 mrbencode74(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 nextpat74(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 boxit74(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 mrbencode74(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 fetchit74(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 mrbencode74(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 nextpat74(misub,k,nord,iflag) + enddo + endif + +998 continue +! Re-order the codeword to [message bits][parity bits] format. + cw(indices)=cw + hdec(indices)=hdec + message74=cw(1:74) + call get_crc24(message74,nbadcrc) + if(nbadcrc.ne.0) nhardmin=-nhardmin + + return +end subroutine osd174_74 + +subroutine mrbencode74(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 mrbencode74 + +subroutine nextpat74(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 nextpat74 + +subroutine boxit74(reset,e2,ntau,npindex,i1,i2) + integer*1 e2(1:ntau) + integer indexes(5000,2),fp(0:525000),np(5000) + 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 boxit74 + +subroutine fetchit74(reset,e2,ntau,i1,i2) + integer indexes(5000,2),fp(0:525000),np(5000) + 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 fetchit74 +