diff --git a/lib/bpdecode128_90.f90 b/lib/bpdecode128_90.f90 index 763b91963..c48ba1fa8 100644 --- a/lib/bpdecode128_90.f90 +++ b/lib/bpdecode128_90.f90 @@ -8,18 +8,18 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter) integer*1 cw(N),apmask(N) integer*1 decoded(K) integer*1 message77(77) - integer Nm(12,M) - integer Mn(4,N) - integer nrw(M),ncw(N) + integer Nm(11,M) + integer Mn(3,N) + integer nrw(M) integer synd(M) real tov(4,N) - real toc(12,M) - real tanhtoc(12,M) + real toc(11,M) + real tanhtoc(11,M) real zn(N) real llr(N) real Tmn - include "ldpc_128_90_b_reordered_parity.f90" + include "ldpc_128_90_reordered_parity.f90" decoded=0 toc=0 @@ -39,7 +39,7 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter) ! 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),i)) + zn(i)=llr(i)+sum(tov(1:ncw,i)) else zn(i)=llr(i) endif @@ -74,7 +74,7 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter) ncnt=ncnt+1 endif ! write(*,*) iter,ncheck,nd,ncnt - if( ncnt .ge. 5 .and. iter .ge. 10 .and. ncheck .gt. 15) then + if( ncnt .ge. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) then nharderror=-1 return endif @@ -86,7 +86,7 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter) do i=1,nrw(j) ibj=Nm(i,j) toc(i,j)=zn(ibj) - do kk=1,4 ! subtract off what the bit had received from the check + 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 @@ -96,11 +96,11 @@ subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter) ! send messages from check nodes to variable nodes do i=1,M - tanhtoc(1:12,i)=tanh(-toc(1:12,i)/2) + tanhtoc(1:11,i)=tanh(-toc(1:11,i)/2) enddo do j=1,N - do i=1,ncw(j) + 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) diff --git a/lib/bpdecode128_90.f90.save b/lib/bpdecode128_90.f90.save new file mode 100644 index 000000000..69537afbb --- /dev/null +++ b/lib/bpdecode128_90.f90.save @@ -0,0 +1,116 @@ +subroutine bpdecode128_90(llr,apmask,maxiterations,message77,cw,nharderror,iter) +! +! A log-domain belief propagation decoder for the (128,90) code. +! + use iso_c_binding, only: c_loc,c_size_t + use crc + integer, parameter:: N=128, K=90, M=N-K + integer*1 cw(N),apmask(N) + integer*1 decoded(K) + integer*1 message77(77) + integer Nm(12,M) + integer Mn(4,N) + integer nrw(M),ncw(N) + integer synd(M) + real tov(4,N) + real toc(12,M) + real tanhtoc(12,M) + real zn(N) + real llr(N) + real Tmn + + include "ldpc_128_90_b_reordered_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 + + 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),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 + decoded=cw(1:K) + call chkcrc13a(decoded,nbadcrc) + if(nbadcrc.eq.0) then + message77=decoded(1:77) + nharderror=count( (2*cw-1)*llr .lt. 0.0 ) + 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. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) 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,4 ! 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:12,i)=tanh(-toc(1:12,i)/2) + enddo + + do j=1,N + do i=1,ncw(j) + 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 bpdecode128_90 diff --git a/lib/encode_128_90.f90 b/lib/encode_128_90.f90 new file mode 100644 index 000000000..e59cf11b6 --- /dev/null +++ b/lib/encode_128_90.f90 @@ -0,0 +1,58 @@ +subroutine encode_128_90(message77,codeword) +! +! Add a 13-bit CRC to a 77-bit message and return a 128-bit codeword +! +use, intrinsic :: iso_c_binding +use iso_c_binding, only: c_loc,c_size_t +use crc + +integer, parameter:: N=128, K=90, M=N-K +character*90 tmpchar +integer*1 codeword(N) +integer*1 gen(M,K) +integer*1 message77(77),message(K) +integer*1 pchecks(M) +integer*1, target :: i1MsgBytes(12) +include "ldpc_128_90_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,23 + read(g(i)(j:j),"(Z1)") istr + ibmax=4 + if(j.eq.23) 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 + +! Add 13 bit CRC to form 90-bit message+CRC13 +write(tmpchar,'(77i1)') message77 +tmpchar(78:80)='000' +i1MsgBytes=0 +read(tmpchar,'(10b8)') i1MsgBytes(1:10) +ncrc13 = crc13 (c_loc (i1MsgBytes), 12) +write(tmpchar(78:90),'(b13)') ncrc13 +read(tmpchar,'(90i1)') message + +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 encode_128_90 diff --git a/lib/genmsk_128_90.f90 b/lib/genmsk_128_90.f90 new file mode 100644 index 000000000..8e6e4ba20 --- /dev/null +++ b/lib/genmsk_128_90.f90 @@ -0,0 +1,114 @@ +subroutine genmsk_128_90(msg0,mygrid,ichk,bcontest,msgsent,i4tone,itype) +! s8 + 48bits + s8 + 80 bits = 144 bits (72ms message duration) +! +! Encode an MSK144 message +! Input: +! - msg0 requested message to be transmitted +! - ichk if ichk=1, return only msgsent +! if ichk.ge.10000, set imsg=ichk-10000 for short msg +! - msgsent message as it will be decoded +! - i4tone array of audio tone values, 0 or 1 +! - itype message type +! 1 = standard message "Call_1 Call_2 Grid/Rpt" +! 2 = type 1 prefix +! 3 = type 1 suffix +! 4 = type 2 prefix +! 5 = type 2 suffix +! 6 = free text (up to 13 characters) +! 7 = short message " Rpt" + + use iso_c_binding, only: c_loc,c_size_t + use packjt77 + character*37 msg0 + character*37 message !Message to be generated + character*37 msgsent !Message as it will be received + character*77 c77 + character*6 mygrid + integer*4 i4tone(144) + integer*1 codeword(128) + integer*1 msgbits(77) + integer*1 bitseq(144) !Tone #s, data and sync (values 0-1) + integer*1 s8(8) + logical bcontest + real*8 pp(12) + real*8 xi(864),xq(864),pi,twopi + data s8/0,1,1,1,0,0,1,0/ + equivalence (ihash,i1hash) + logical first + data first/.true./ + save + + if(first) then + first=.false. + nsym=128 + pi=4.0*atan(1.0) + twopi=8.*atan(1.0) + do i=1,12 + pp(i)=sin((i-1)*pi/12) + enddo + endif + + if(msg0(1:1).eq.'@') then !Generate a fixed tone + read(msg0(2:5),*,end=1,err=1) nfreq !at specified frequency + go to 2 +1 nfreq=1000 +2 i4tone(1)=nfreq + else + message=msg0 + do i=1,22 + if(ichar(message(i:i)).eq.0) then + message(i:)=' ' + exit + endif + enddo + + do i=1,22 !Strip leading blanks + if(message(1:1).ne.' ') exit + message=message(i+1:) + enddo + + if(message(1:1).eq.'<') then + call genmsk40(message,msgsent,ichk,i4tone,itype) + if(itype.lt.0) go to 999 + i4tone(41)=-40 + go to 999 + endif + + call pack77(message,i3,n3,c77) + call unpack77(c77,msgsent) !Unpack to get msgsent + if(ichk.eq.1) go to 999 + read(c77,"77i1") msgbits + call encode_128_90(msgbits,codeword) + +!Create 144-bit channel vector: +!8-bit sync word + 48 bits + 8-bit sync word + 80 bits + bitseq=0 + bitseq(1:8)=s8 + bitseq(9:56)=codeword(1:48) + bitseq(57:64)=s8 + bitseq(65:144)=codeword(49:128) + bitseq=2*bitseq-1 + + xq(1:6)=bitseq(1)*pp(7:12) !first bit is mapped to 1st half-symbol on q + do i=1,71 + is=(i-1)*12+7 + xq(is:is+11)=bitseq(2*i+1)*pp + enddo + xq(864-5:864)=bitseq(1)*pp(1:6) !last half symbol + do i=1,72 + is=(i-1)*12+1 + xi(is:is+11)=bitseq(2*i)*pp + enddo +! Map I and Q to tones. + i4tone=0 + do i=1,72 + i4tone(2*i-1)=(bitseq(2*i)*bitseq(2*i-1)+1)/2; + i4tone(2*i)=-(bitseq(2*i)*bitseq(mod(2*i,144)+1)-1)/2; + enddo + endif + +! Flip polarity + i4tone=-i4tone+1 + +999 return +end subroutine genmsk_128_90 diff --git a/lib/ldpcsim128_90.f90 b/lib/ldpcsim128_90.f90 index 6ddfda22d..ebe759806 100644 --- a/lib/ldpcsim128_90.f90 +++ b/lib/ldpcsim128_90.f90 @@ -1,132 +1,134 @@ -program ldpcsim +program ldpcsim128_90 -use packjt -integer, parameter:: NRECENT=10, N=128, K=90, M=N-K -character*12 recent_calls(NRECENT) -character*22 msg,msgsent,msgreceived -character*96 tmpchar -character*8 arg -integer*1 codeword(N), message77(77) -integer*1 apmask(N),cw(N) -integer*1 msgbits(77) -integer*4 i4Msg6BitWords(13) -integer nerrtot(0:N),nerrdec(0:N) -real*8 rxdata(N), rxavgd(N) -real llr(N) +! Simulate the performance of the (128,90) code that is used in +! the second incarnation of MSK144. -do i=1,NRECENT - recent_calls(i)=' ' -enddo -nerrtot=0 -nerrdec=0 + use packjt77 + integer, parameter:: NRECENT=10, N=128, K=90, M=N-K + character*12 recent_calls(NRECENT) + character*37 msg,msgsent,msgreceived + character*77 c77 + character*8 arg + integer*1 codeword(N), message77(77) + integer*1 apmask(N),cw(N) + integer*1 msgbits(77) + integer*4 i4Msg6BitWords(13) + integer nerrtot(0:N),nerrdec(0:N) + real*8 rxdata(N), rxavgd(N) + real llr(N) -nargs=iargc() -if(nargs.ne.4) then - print*,'Usage: ldpcsim niter navg #trials s ' - print*,'eg: ldpcsim 10 1 1000 0.75' - return -endif -call getarg(1,arg) -read(arg,*) max_iterations -call getarg(2,arg) -read(arg,*) navg -call getarg(3,arg) -read(arg,*) ntrials -call getarg(4,arg) -read(arg,*) s + do i=1,NRECENT + recent_calls(i)=' ' + enddo + nerrtot=0 + nerrdec=0 -rate=real(K)/real(N) + nargs=iargc() + if(nargs.ne.4) then + print*,'Usage: ldpcsim niter navg #trials s ' + print*,'eg: ldpcsim 10 1 1000 0.75' + return + endif + call getarg(1,arg) + read(arg,*) max_iterations + call getarg(2,arg) + read(arg,*) navg + call getarg(3,arg) + read(arg,*) ntrials + call getarg(4,arg) + read(arg,*) s -write(*,*) "rate: ",rate -write(*,*) "niter= ",max_iterations," navg= ",navg," s= ",s + rate=real(77)/real(N) + + write(*,*) "rate: ",rate + write(*,*) "niter= ",max_iterations," navg= ",navg," s= ",s + + msg="K1ABC RR73; W9XYZ -12" + i3=0 + n3=1 + call pack77(msg,i3,n3,c77) + call unpack77(c77,msgsent) + read(c77,'(77i1)') msgbits -!msg="K9AN K1JT EN50" -msg="G4WJS K1JT FN20" - call packmsg(msg,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes - call unpackmsg(i4Msg6BitWords,msgsent,.false.,' ') !Unpack to get msgsent write(*,*) "message sent ",msgsent - tmpchar=' ' - write(tmpchar,'(12b6.6)') i4Msg6BitWords(1:12) - tmpchar(73:77)="00000" !i5bit - read(tmpchar,'(77i1)') msgbits(1:77) - write(*,*) 'msgbits' - write(*,'(28i1,1x,28i1,1x,16i1,1x,5i1)') msgbits + write(*,'(77i1)') msgbits ! msgbits is the 77-bit message, codeword is 128 bits - call encode128_90(msgbits,codeword) + call encode_128_90(msgbits,codeword) - call init_random_seed() +! call init_random_seed() -write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr" -do idb = 14,-6,-1 - db=idb/2.0-1.0 - sigma=1/sqrt( 2*rate*(10**(db/10.0)) ) - ngood=0 - nue=0 - nbadcrc=0 - nsumerr=0 + write(*,*) "Eb/N0 SNR2500 ngood nundetected sigma psymerr" + do idb = 14,0,-1 + db=idb/2.0-1.0 + sigma=1/sqrt( 2*rate*(10**(db/10.0)) ) + ngood=0 + nue=0 + nbadcrc=0 + nsumerr=0 - do itrial=1, ntrials - rxavgd=0d0 - do iav=1,navg - call sgran() + do itrial=1, ntrials + rxavgd=0d0 + do iav=1,navg +! call sgran() ! Create a realization of a noisy received word - do i=1,N - rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran() + do i=1,N + rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran() + enddo + rxavgd=rxavgd+rxdata enddo - rxavgd=rxavgd+rxdata - enddo - rxdata=rxavgd - 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 + rxdata=rxavgd + 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 - rxav=sum(rxdata)/N - rx2av=sum(rxdata*rxdata)/N - rxsig=sqrt(rx2av-rxav*rxav) - rxdata=rxdata/rxsig + rxav=sum(rxdata)/N + rx2av=sum(rxdata*rxdata)/N + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig ! 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 + if( s .lt. 0 ) then + ss=sigma + else + ss=s + endif - llr=2.0*rxdata/(ss*ss) + llr=2.0*rxdata/(ss*ss) - apmask=0 + apmask=0 ! max_iterations is max number of belief propagation iterations - call bpdecode128_90(llr, apmask, max_iterations, message77, cw, nharderrors, niterations) + call bpdecode128_90(llr, apmask, max_iterations, message77, cw, nharderrors, niterations) ! If the decoder finds a valid codeword, nharderrors will be .ge. 0. - if( nharderrors .ge. 0 ) then - call extractmessage77(message77,msgreceived) - nhw=count(cw.ne.codeword) - if(nhw.eq.0) then ! this is a good decode + if( nharderrors .ge. 0 ) then + write(c77,'(77i1)') message77 + call unpack77(c77,msgreceived) + nhw=count(cw.ne.codeword) + if(nhw.eq.0) then ! this is a good decode ngood=ngood+1 nerrdec(nerr)=nerrdec(nerr)+1 - else ! this is an undetected error + else ! this is an undetected error nue=nue+1 - endif - endif - nsumerr=nsumerr+nerr + endif + endif + nsumerr=nsumerr+nerr + enddo + + snr2500=db+10*log10(rate*2000.0/2500.0) ! symbol rate is 2000 s^-1 and ref BW is 2500 Hz. + pberr=real(nsumerr)/real(ntrials*N) + write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,7x,f5.2,3x,e10.3)") db,snr2500,ngood,nue,ss,pberr + enddo - snr2500=db-2.5 - pberr=real(nsumerr)/real(ntrials*N) - write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,7x,f5.2,3x,e10.3)") db,snr2500,ngood,nue,ss,pberr - -enddo + open(unit=23,file='nerrhisto.dat',status='unknown') + do i=0,N + 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=23,file='nerrhisto.dat',status='unknown') -do i=0,N - 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 ldpcsim +end program ldpcsim128_90 diff --git a/lib/ldpcsim144.f90 b/lib/ldpcsim144.f90 index 31271d5a9..20ce88c18 100644 --- a/lib/ldpcsim144.f90 +++ b/lib/ldpcsim144.f90 @@ -41,7 +41,7 @@ call getarg(4,arg) read(arg,*) s ! don't count hash bits as data bits -K=80 +K=72 N=128 rate=real(K)/real(N) @@ -161,7 +161,7 @@ do idb = 14,-6,-1 endif endif enddo - snr2500=db-3.0 + snr2500=db-3.5 write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2)") db,snr2500,ngood,nue,nbadhash,ss enddo