diff --git a/CMakeLists.txt b/CMakeLists.txt index eaf45f1be..e1ccf1389 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -603,17 +603,17 @@ set (wsjt_FSRCS lib/fsk4hf/genwsprcpm.f90 lib/fsk4hf/encode204.f90 lib/fsk4hf/decode174_74.f90 - lib/fsk4hf/decode174_101.f90 + lib/fsk4hf/decode240_101.f90 lib/fsk4hf/ldpcsim174_91.f90 lib/fsk4hf/ldpcsim174_74.f90 - lib/fsk4hf/ldpcsim174_101.f90 + lib/fsk4hf/ldpcsim240_101.f90 lib/fsk4hf/get_crc24.f90 lib/fsk4hf/encode174_74.f90 - lib/fsk4hf/encode174_101.f90 + lib/fsk4hf/encode240_101.f90 lib/fsk4hf/bpdecode174_74.f90 - lib/fsk4hf/bpdecode174_101.f90 + lib/fsk4hf/bpdecode240_101.f90 lib/fsk4hf/osd174_74.f90 - lib/fsk4hf/osd174_101.f90 + lib/fsk4hf/osd240_101.f90 lib/fsk4hf/wspr4sim.f90 lib/fsk4hf/genwspr4.f90 lib/fsk4hf/gen_wspr4wave.f90 @@ -621,7 +621,9 @@ set (wsjt_FSRCS lib/fsk4hf/get_wspr4_bitmetrics.f90 lib/fsk4hf/ft4slowsim.f90 lib/fsk4hf/genft4slow.f90 - lib/fsk4hf/decode174_101.f90 + lib/fsk4hf/decode240_101.f90 + lib/fsk4hf/ft4sd.f90 + lib/fsk4hf/get_ft4s_bitmetrics.f90 ) # temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit @@ -1368,8 +1370,8 @@ target_link_libraries (ldpcsim174_91 wsjt_fort wsjt_cxx) add_executable (ldpcsim174_74 lib/fsk4hf/ldpcsim174_74.f90 wsjtx.rc) target_link_libraries (ldpcsim174_74 wsjt_fort wsjt_cxx) -add_executable (ldpcsim174_101 lib/fsk4hf/ldpcsim174_101.f90 wsjtx.rc) -target_link_libraries (ldpcsim174_101 wsjt_fort wsjt_cxx) +add_executable (ldpcsim240_101 lib/fsk4hf/ldpcsim240_101.f90 wsjtx.rc) +target_link_libraries (ldpcsim240_101 wsjt_fort wsjt_cxx) add_executable (wspr4sim lib/fsk4hf/wspr4sim.f90 wsjtx.rc) target_link_libraries (wspr4sim wsjt_fort wsjt_cxx) @@ -1380,6 +1382,9 @@ target_link_libraries (wspr4d wsjt_fort wsjt_cxx) add_executable (ft4slowsim lib/fsk4hf/ft4slowsim.f90 wsjtx.rc) target_link_libraries (ft4slowsim wsjt_fort wsjt_cxx) +add_executable (ft4sd lib/fsk4hf/ft4sd.f90 wsjtx.rc) +target_link_libraries (ft4sd wsjt_fort wsjt_cxx) + endif(WSJT_BUILD_UTILS) # build the main application diff --git a/lib/fsk4hf/bpdecode240_101.f90 b/lib/fsk4hf/bpdecode240_101.f90 new file mode 100644 index 000000000..1e2adbb68 --- /dev/null +++ b/lib/fsk4hf/bpdecode240_101.f90 @@ -0,0 +1,111 @@ +subroutine bpdecode240_101(llr,apmask,maxiterations,message101,cw,nharderror,iter,ncheck) +! +! A log-domain belief propagation decoder for the (240,101) code. +! + integer, parameter:: N=240, K=101, M=N-K + integer*1 cw(N),apmask(N) + integer*1 decoded(K) + integer*1 message101(101) + 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_240_101_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:101) + call get_crc24(decoded,101,nbadcrc) + nharderror=count( (2*cw-1)*llr .lt. 0.0 ) + if(nbadcrc.eq.0) then + message101=decoded(1:101) + 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 bpdecode240_101 diff --git a/lib/fsk4hf/decode240_101.f90 b/lib/fsk4hf/decode240_101.f90 new file mode 100644 index 000000000..71cbc8ff3 --- /dev/null +++ b/lib/fsk4hf/decode240_101.f90 @@ -0,0 +1,133 @@ +subroutine decode240_101(llr,Keff,ndeep,apmask,maxsuper,message101,cw,nharderror,iter,ncheck,dmin,isuper) +! +! A hybrid bp/osd decoder for the (240,101) code. +! + integer, parameter:: N=240, K=101, M=N-K + integer*1 cw(N),apmask(N) + integer*1 decoded(K) + integer*1 nxor(N),hdec(N) + integer*1 message101(101) + 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),zsum(N) + real llr(N) + real Tmn + + include "ldpc_240_101_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 + + maxiterations=1 + + zsum=0.0 + do isuper=1,maxsuper + + 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 + zsum=zsum+zn +! 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 - if crc is good, return it + decoded=cw(1:K) + call get_crc24(decoded,74,nbadcrc) + nharderror=count( (2*cw-1)*llr .lt. 0.0 ) + if(nbadcrc.eq.0) then + message101=decoded(1:101) + dmin=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. 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 ! bp iterations + + call osd240_101(zsum,Keff,apmask,ndeep,message101,cw,nharderror,dminosd) + if(nharderror.gt.0) then + hdec=0 + where(llr .ge. 0) hdec=1 + nxor=ieor(hdec,cw) + dmin=sum(nxor*abs(llr)) + return + endif + enddo ! super iterations + + nharderror=-1 + + return +end subroutine decode240_101 diff --git a/lib/fsk4hf/encode240_101.f90 b/lib/fsk4hf/encode240_101.f90 new file mode 100644 index 000000000..da0021df3 --- /dev/null +++ b/lib/fsk4hf/encode240_101.f90 @@ -0,0 +1,46 @@ +subroutine encode240_101(message,codeword) + use, intrinsic :: iso_c_binding + use iso_c_binding, only: c_loc,c_size_t + use crc + + integer, parameter:: N=240, K=101, M=N-K + character*24 c24 + integer*1 codeword(N) + integer*1 gen(M,K) + integer*1 message(K) + integer*1 pchecks(M) + integer*4 ncrc24 + include "ldpc_240_101_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,26 + read(g(i)(j:j),"(Z1)") istr + ibmax=4 + if(j.eq.26) ibmax=1 + 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 encode240_101 diff --git a/lib/fsk4hf/ft4s_params.f90 b/lib/fsk4hf/ft4s_params.f90 new file mode 100644 index 000000000..510d7e505 --- /dev/null +++ b/lib/fsk4hf/ft4s_params.f90 @@ -0,0 +1,16 @@ +! FT4A +! LDPC(240,101)/CRC24 code, four 4x4 Costas arrays for sync, ramp-up and ramp-down symbols + +parameter (KK=77) !Information bits (77 + CRC24) +parameter (ND=120) !Data symbols +parameter (NS=24) !Sync symbols +parameter (NN=NS+ND) !Sync and data symbols (144) +parameter (NN2=NS+ND+2) !Total channel symbols (146) +parameter (NSPS=9600) !Samples per symbol at 12000 S/s +parameter (NZ=NSPS*NN) !Sync and Data samples (1,382,400) +parameter (NZ2=NSPS*NN2) !Total samples in shaped waveform (1,397,760) +parameter (NMAX=408*3456) !Samples in iwave (1,410,048) +parameter (NFFT1=4*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra +parameter (NSTEP=NSPS) !Coarse time-sync step size +parameter (NHSYM=(NMAX-NFFT1)/NSTEP) !Number of symbol spectra (1/4-sym steps) +parameter (NDOWN=32) !Downsample factor diff --git a/lib/fsk4hf/ft4sd.f90 b/lib/fsk4hf/ft4sd.f90 new file mode 100644 index 000000000..e2c08e6ef --- /dev/null +++ b/lib/fsk4hf/ft4sd.f90 @@ -0,0 +1,459 @@ +program ft4sd + +! Decode ft4slow data read from *.c2 or *.wav files. + + use packjt77 + include 'ft4s_params.f90' + parameter (NSPS2=NSPS/32) + character arg*8,cbits*50,infile*80,fname*16,datetime*11 + character ch1*1,ch4*4,cseq*31 + character*22 decodes(100) + character*37 msg + character*120 data_dir + character*77 c77 + complex c2(0:NMAX/32-1) !Complex waveform + complex cframe(0:144*NSPS2-1) !Complex waveform + complex cd(0:144*20-1) !Complex waveform + real*8 fMHz + real llr(240),llra(240),llrb(240),llrc(240),llrd(240) + real candidates(100,2) + real bitmetrics(288,4) + integer ihdr(11) + integer*2 iwave(NMAX) !Generated full-length waveform + integer*1 apmask(240),cw(240) + integer*1 hbits(288) + integer*1 message101(101) + logical badsync,unpk77_success + + fs=12000.0/NDOWN !Sample rate + dt=1.0/fs !Sample interval (s) + tt=NSPS*dt !Duration of "itone" symbols (s) + txt=NZ*dt !Transmission length (s) + + nargs=iargc() + if(nargs.lt.1) then + print*,'Usage: ft4sd [-a ] [-f fMHz] file1 [file2 ...]' + go to 999 + endif + iarg=1 + data_dir="." + call getarg(iarg,arg) + if(arg(1:2).eq.'-a') then + call getarg(iarg+1,data_dir) + iarg=iarg+2 + endif + call getarg(iarg,arg) + if(arg(1:2).eq.'-f') then + call getarg(iarg+1,arg) + read(arg,*) fMHz + iarg=iarg+2 + endif + + ngood=0 + do ifile=iarg,nargs + call getarg(ifile,infile) + open(10,file=infile,status='old',access='stream') + j1=index(infile,'.c2') + j2=index(infile,'.wav') + if(j1.gt.0) then + read(10,end=999) fname,ntrmin,fMHz,c2 + read(fname(8:11),*) nutc + write(datetime,'(i11)') nutc + else if(j2.gt.0) then + read(10,end=999) ihdr,iwave + read(infile(j2-4:j2-1),*) nutc + datetime=infile(j2-11:j2-1) + call ft4s_downsample(iwave,c2) + else + print*,'Wrong file format?' + go to 999 + endif + close(10) + + fa=-100.0 + fb=100.0 + fs=12000.0/32.0 + npts=120*12000.0/32.0 + + call getcandidate_ft4s(c2,npts,fs,fa,fb,ncand,candidates) !First approx for freq + + ndecodes=0 + do icand=1,ncand + fc0=candidates(icand,1) + xsnr=candidates(icand,2) +!write(*,*) 'candidates ',icand,fc0,xsnr + do isync=0,1 + + del=1.5*fs/300.0 + if(isync.eq.0) then + fc1=fc0-del + is0=375 + ishw=350 + isst=30 + ifhw=10 + df=.1 + else if(isync.eq.1) then + fc1=fc2 + is0=isbest + ishw=100 + isst=10 + ifhw=10 + df=.02 + endif + smax=0.0 + do if=-ifhw,ifhw + fc=fc1+df*if + do istart=max(1,is0-ishw),is0+ishw,isst + call coherent_sync_ft4s(c2,istart,fc,1,sync) + if(sync.gt.smax) then + fc2=fc + isbest=istart + smax=sync + endif + enddo + enddo +! write(*,*) ifile,icand,isync,fc1+del,fc2+del,isbest,smax + enddo + +! if(smax .lt. 100.0 ) cycle +!isbest=375 +!fc2=-del + do ijitter=0,2 + if(ijitter.eq.0) ioffset=0 + if(ijitter.eq.1) ioffset=45 + if(ijitter.eq.2) ioffset=-45 + is0=isbest+ioffset + if(is0.lt.0) cycle + cframe=c2(is0:is0+144*300-1) + call downsample_ft4s(cframe,fc2,cd) + s2=sum(cd*conjg(cd))/(20*144) + cd=cd/sqrt(s2) + call get_ft4s_bitmetrics(cd,bitmetrics,badsync) + + hbits=0 + where(bitmetrics(:,1).ge.0) hbits=1 + ns1=count(hbits( 1: 8).eq.(/0,0,0,1,1,0,1,1/)) + ns2=count(hbits( 57: 64).eq.(/0,1,0,0,1,1,1,0/)) + ns3=count(hbits(113:120).eq.(/1,1,1,0,0,1,0,0/)) + ns4=count(hbits(169:176).eq.(/1,0,1,1,0,0,0,1/)) + ns5=count(hbits(225:232).eq.(/0,0,1,1,1,0,0,1/)) + ns6=count(hbits(281:288).eq.(/0,1,1,1,0,0,1,0/)) + nsync_qual=ns1+ns2+ns3+ns4+ns5+ns6 +! if(nsync_qual.lt. 20) cycle + + scalefac=2.83 + llra( 1: 48)=bitmetrics( 9: 56, 1) + llra( 49: 96)=bitmetrics( 65:112, 1) + llra( 97:144)=bitmetrics(121:168, 1) + llra(145:192)=bitmetrics(177:224, 1) + llra(193:240)=bitmetrics(233:280, 1) + llra=scalefac*llra + llrb( 1: 48)=bitmetrics( 9: 56, 2) + llrb( 49: 96)=bitmetrics( 65:112, 2) + llrb( 97:144)=bitmetrics(121:168, 2) + llrb(145:192)=bitmetrics(177:224, 2) + llrb(193:240)=bitmetrics(233:280, 2) + llrb=scalefac*llrb + llrc( 1: 48)=bitmetrics( 9: 56, 3) + llrc( 49: 96)=bitmetrics( 65:112, 3) + llrc( 97:144)=bitmetrics(121:168, 3) + llrc(145:192)=bitmetrics(177:224, 3) + llrc(193:240)=bitmetrics(233:280, 3) + llrc=scalefac*llrc + llrd( 1: 48)=bitmetrics( 9: 56, 4) + llrd( 49: 96)=bitmetrics( 65:112, 4) + llrd( 97:144)=bitmetrics(121:168, 4) + llrd(145:192)=bitmetrics(177:224, 4) + llrd(193:240)=bitmetrics(233:280, 4) + llrd=scalefac*llrd + apmask=0 + max_iterations=40 + + do itry=4,1,-1 + if(itry.eq.1) llr=llra + if(itry.eq.2) llr=llrb + if(itry.eq.3) llr=llrc + if(itry.eq.4) llr=llrd + nhardbp=0 + nhardosd=0 + dmin=0.0 + call bpdecode240_101(llr,apmask,max_iterations,message101,cw,nhardbp,niterations,nchecks) + Keff=91 +! if(nhardbp.lt.0) call osd240_101(llr,Keff,apmask,5,message101,cw,nhardosd,dmin) + maxsuperits=2 + ndeep=3 ! use ndeep=3 with Keff=91 + if(nhardbp.lt.0) then +! call osd240_101(llr,Keff,apmask,ndeep,message101,cw,nhardosd,dmin) + call decode240_101(llr,Keff,ndeep,apmask,maxsuperits,message101,cw,nhardosd,iter,ncheck,dmin,isuper) + endif + if(nhardbp.ge.0 .or. nhardosd.ge.0) then + write(c77,'(77i1)') message101(1:77) + call unpack77(c77,0,msg,unpk77_success) + if(unpk77_success .and. index(msg,'K9AN').gt.0) then + ngood=ngood+1 + write(*,1100) ifile,icand,xsnr,isbest/375.0-1.0,1500.0+fc2+del,msg(1:14),itry,nhardbp,nhardosd,dmin,ijitter +1100 format(i5,2x,i5,2x,f6.1,2x,f6.2,2x,f8.2,2x,a14,i4,i4,i4,f7.2,i6) + goto 2002 + else + cycle + endif + endif + enddo ! metrics + enddo ! istart jitter + enddo !candidate list +2002 continue + enddo !files + nfiles=nargs-iarg+1 + write(*,*) 'nfiles: ',nfiles,' ngood: ',ngood + write(*,1120) +1120 format("") + +999 end program ft4sd + +subroutine coherent_sync_ft4s(cd0,i0,f0,itwk,sync) + +! Compute sync power for a complex, downsampled FT4s signal. + + include 'ft4s_params.f90' + parameter(NP=NMAX/NDOWN,NSS=NSPS/NDOWN) + complex cd0(0:NP-1) + complex csynca(4*NSS),csyncb(4*NSS) + complex csyncc(4*NSS),csyncd(4*NSS) + complex csynce(4*NSS),csyncf(4*NSS) + complex csync2(4*NSS) + complex ctwk(4*NSS) + complex z1,z2,z3,z4,z5,z6 + logical first + integer icos4a(0:3),icos4b(0:3) + integer icos4c(0:3),icos4d(0:3) + integer icos4e(0:3),icos4f(0:3) + data icos4a/0,1,3,2/ + data icos4b/1,0,2,3/ + data icos4c/2,3,1,0/ + data icos4d/3,2,0,1/ + data icos4e/0,2,3,1/ + data icos4f/1,2,0,3/ + data first/.true./ + save first,twopi,csynca,csyncb,csyncc,csyncd,csynce,csyncf,fac + + p(z1)=(real(z1*fac)**2 + aimag(z1*fac)**2)**0.5 !Statement function for power + + if( first ) then + twopi=8.0*atan(1.0) + k=1 + phia=0.0 + phib=0.0 + phic=0.0 + phid=0.0 + phie=0.0 + phif=0.0 + do i=0,3 + dphia=twopi*icos4a(i)/real(NSS) + dphib=twopi*icos4b(i)/real(NSS) + dphic=twopi*icos4c(i)/real(NSS) + dphid=twopi*icos4d(i)/real(NSS) + dphie=twopi*icos4e(i)/real(NSS) + dphif=twopi*icos4f(i)/real(NSS) + do j=1,NSS + csynca(k)=cmplx(cos(phia),sin(phia)) + csyncb(k)=cmplx(cos(phib),sin(phib)) + csyncc(k)=cmplx(cos(phic),sin(phic)) + csyncd(k)=cmplx(cos(phid),sin(phid)) + csynce(k)=cmplx(cos(phie),sin(phie)) + csyncf(k)=cmplx(cos(phif),sin(phif)) + phia=mod(phia+dphia,twopi) + phib=mod(phib+dphib,twopi) + phic=mod(phic+dphic,twopi) + phid=mod(phid+dphid,twopi) + phie=mod(phie+dphie,twopi) + phif=mod(phif+dphif,twopi) + k=k+1 + enddo + enddo + first=.false. + fac=1.0/(4.0*NSS) + endif + + i1=i0 !four Costas arrays + i2=i0+28*NSS + i3=i0+56*NSS + i4=i0+84*NSS + i5=i0+112*NSS + i6=i0+140*NSS + + z1=0. + z2=0. + z3=0. + z4=0. + z5=0. + z6=0. + + if(itwk.eq.1) then + dt=1/(12000.0/32.0) + dphi=twopi*f0*dt + phi=0.0 + do i=1,4*NSS + ctwk(i)=cmplx(cos(phi),sin(phi)) + phi=mod(phi+dphi,twopi) + enddo + endif + + if(itwk.eq.1) csync2=ctwk*csynca !Tweak the frequency + if(i1.ge.0 .and. i1+4*NSS-1.le.NP-1) then + z1=sum(cd0(i1:i1+4*NSS-1)*conjg(csync2)) + elseif( i1.lt.0 ) then + npts=(i1+4*NSS-1)/2 + if(npts.le.40) then + z1=0. + else + z1=sum(cd0(0:i1+4*NSS-1)*conjg(csync2(4*NSS-npts:))) + endif + endif + + if(itwk.eq.1) csync2=ctwk*csyncb !Tweak the frequency + if(i2.ge.0 .and. i2+4*NSS-1.le.NP-1) then + z2=sum(cd0(i2:i2+4*NSS-1)*conjg(csync2)) + endif + + if(itwk.eq.1) csync2=ctwk*csyncc !Tweak the frequency + if(i3.ge.0 .and. i3+4*NSS-1.le.NP-1) then + z3=sum(cd0(i3:i3+4*NSS-1)*conjg(csync2)) + endif + + if(itwk.eq.1) csync2=ctwk*csyncd !Tweak the frequency + if(i4.ge.0 .and. i4+4*NSS-1.le.NP-1) then + z4=sum(cd0(i4:i4+4*NSS-1)*conjg(csync2)) + endif + + if(itwk.eq.1) csync2=ctwk*csynce !Tweak the frequency + if(i5.ge.0 .and. i5+4*NSS-1.le.NP-1) then + z5=sum(cd0(i5:i5+4*NSS-1)*conjg(csync2)) + endif + + if(itwk.eq.1) csync2=ctwk*csyncf !Tweak the frequency + if(i6.ge.0 .and. i6+4*NSS-1.le.NP-1) then + z6=sum(cd0(i6:i6+4*NSS-1)*conjg(csync2)) + elseif( i6+4*NSS-1.gt.NP-1 ) then + npts=(NP-1-i6+1) + if(npts.le.40) then + z6=0. + else + z6=sum(cd0(i6:i6+npts-1)*conjg(csync2(1:npts))) + endif + endif + + sync = p(z1) + p(z2) + p(z3) + p(z4) + p(z5) + p(z6) + + return +end subroutine coherent_sync_ft4s + +subroutine downsample_ft4s(ci,f0,co) + parameter(NI=144*300,NH=NI/2,NO=NI/15) ! downsample from 315 samples per symbol to 20 + complex ci(0:NI-1),ct(0:NI-1) + complex co(0:NO-1) + fs=12000.0/32.0 + df=fs/NI + ct=ci + call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain + i0=nint(f0/df) + ct=cshift(ct,i0) + co=0.0 + co(0)=ct(0) + b=16.0 + do i=1,NO/2 + arg=(i*df/b)**2 + filt=exp(-arg) + co(i)=ct(i)*filt + co(NO-i)=ct(NI-i)*filt + enddo + co=co/NO + call four2a(co,NO,1,1,1) !c2c FFT back to time domain + return +end subroutine downsample_ft4s + +subroutine getcandidate_ft4s(c,npts,fs,fa,fb,ncand,candidates) + parameter(NFFT1=120*12000/32,NH1=NFFT1/2,NFFT2=120*12000/320,NH2=NFFT2/2) + complex c(0:npts-1) !Complex waveform + complex cc(0:NFFT1-1) + complex csfil(0:NFFT2-1) + complex cwork(0:NFFT2-1) + real bigspec(0:NFFT2-1) + complex c2(0:NFFT1-1) !Short spectra + real s(-NH1+1:NH1) !Coarse spectrum + real ss(-NH1+1:NH1) !Smoothed coarse spectrum + real candidates(100,2) + integer indx(NFFT2-1) + logical first + data first/.true./ + save first,w,df,csfil + + if(first) then + df=10*fs/NFFT1 + csfil=cmplx(0.0,0.0) + do i=0,NFFT2-1 +! csfil(i)=exp(-((i-NH2)/32.0)**2) ! revisit this + csfil(i)=exp(-((i-NH2)/28.0)**2) ! revisit this + enddo + csfil=cshift(csfil,NH2) + call four2a(csfil,NFFT2,1,-1,1) + first=.false. + endif + + cc=cmplx(0.0,0.0) + cc(0:npts-1)=c; + call four2a(cc,NFFT1,1,-1,1) + cc=abs(cc)**2 + call four2a(cc,NFFT1,1,-1,1) + cwork(0:NH2)=cc(0:NH2)*conjg(csfil(0:NH2)) + cwork(NH2+1:NFFT2-1)=cc(NFFT1-NH2+1:NFFT1-1)*conjg(csfil(NH2+1:NFFT2-1)) + + call four2a(cwork,NFFT2,1,+1,1) + bigspec=cshift(real(cwork),-NH2) + il=NH2+fa/df + ih=NH2+fb/df + nnl=ih-il+1 + call indexx(bigspec(il:il+nnl-1),nnl,indx) + xn=bigspec(il-1+indx(nint(0.3*nnl))) + bigspec=bigspec/xn + ncand=0 + do i=il,ih + if((bigspec(i).gt.bigspec(i-1)).and. & + (bigspec(i).gt.bigspec(i+1)).and. & + (bigspec(i).gt.1.15).and.ncand.lt.100) then + ncand=ncand+1 + candidates(ncand,1)=df*(i-NH2) + candidates(ncand,2)=10*log10(bigspec(i)-1)-28.5 + endif + enddo + return +end subroutine getcandidate_ft4s + +subroutine ft4s_downsample(iwave,c) + +! Input: i*2 data in iwave() at sample rate 12000 Hz +! Output: Complex data in c(), sampled at 375 Hz + + include 'ft4s_params.f90' + parameter (NFFT2=NMAX/32) + integer*2 iwave(NMAX) + complex c(0:NMAX/32-1) + complex c1(0:NFFT2-1) + complex cx(0:NMAX/2) + real x(NMAX) + equivalence (x,cx) + + df=12000.0/NMAX + x=iwave + call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain + i0=nint(1500.0/df) + c1(0)=cx(i0) + do i=1,NFFT2/2 + c1(i)=cx(i0+i) + c1(NFFT2-i)=cx(i0-i) + enddo + c1=c1/NFFT2 + call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain + c=c1(0:NMAX/32-1) + return +end subroutine ft4s_downsample + diff --git a/lib/fsk4hf/ft4slowsim.f90 b/lib/fsk4hf/ft4slowsim.f90 index 603b679be..f0726868a 100644 --- a/lib/fsk4hf/ft4slowsim.f90 +++ b/lib/fsk4hf/ft4slowsim.f90 @@ -4,7 +4,7 @@ program ft4slowsim use wavhdr use packjt77 - include 'wspr4_params.f90' !Set various constants + include 'ft4s_params.f90' !Set various constants type(hdr) h !Header for .wav file character arg*12,fname*17 character msg37*37,msgsent37*37 @@ -83,8 +83,6 @@ program ft4slowsim if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NMAX,NZ,fs,delay,fspread) c=sig*c wave=real(c) - peak=maxval(abs(wave)) - nslots=1 if(snrdb.lt.90) then do i=1,NMAX !Add gaussian noise at specified SNR xnoise=gran() diff --git a/lib/fsk4hf/genft4slow.f90 b/lib/fsk4hf/genft4slow.f90 index 3f320b4c0..cd58094a3 100644 --- a/lib/fsk4hf/genft4slow.f90 +++ b/lib/fsk4hf/genft4slow.f90 @@ -8,13 +8,12 @@ subroutine genft4slow(msg0,ichk,msgsent,msgbits,i4tone) ! - i4tone array of audio tone values, {0,1,2,3} ! Frame structure: -! s16 + 87symbols + 2 ramp up/down = 105 total channel symbols -! r1 + s4 + d29 + s4 + d29 + s4 + d29 + s4 + r1 +! s4 d24 s4 d24 s4 d24 s4 d24 s4 d24 s4 -! Message duration: TxT = 105*13312/12000 = 116.48 s +! Message duration: TxT = 144*9600/12000 = 115.2 s use packjt77 - include 'wspr4_params.f90' + include 'ft4s_params.f90' character*37 msg0 character*37 message !Message to be generated character*37 msgsent !Message as it will be received @@ -23,13 +22,15 @@ subroutine genft4slow(msg0,ichk,msgsent,msgbits,i4tone) integer*4 i4tone(NN),itmp(ND) integer*1 codeword(2*ND) integer*1 msgbits(101),rvec(77) - integer icos4a(4),icos4b(4),icos4c(4),icos4d(4) + integer icos4a(4),icos4b(4),icos4c(4),icos4d(4),icos4e(4),icos4f(4) integer ncrc24 logical unpk77_success data icos4a/0,1,3,2/ data icos4b/1,0,2,3/ data icos4c/2,3,1,0/ data icos4d/3,2,0,1/ + data icos4e/0,2,3,1/ + data icos4f/1,2,0,3/ data rvec/0,1,0,0,1,0,1,0,0,1,0,1,1,1,1,0,1,0,0,0,1,0,0,1,1,0,1,1,0, & 1,0,0,1,0,1,1,0,0,0,0,1,0,0,0,1,0,1,0,0,1,1,1,1,0,0,1,0,1, & 0,1,0,1,0,1,1,0,1,1,1,1,1,0,0,0,1,0,1/ @@ -65,7 +66,7 @@ subroutine genft4slow(msg0,ichk,msgsent,msgbits,i4tone) entry get_ft4slow_tones_from_101bits(msgbits,i4tone) -2 call encode174_101(msgbits,codeword) +2 call encode240_101(msgbits,codeword) ! Grayscale mapping: ! bits tone @@ -82,12 +83,16 @@ entry get_ft4slow_tones_from_101bits(msgbits,i4tone) enddo i4tone(1:4)=icos4a - i4tone(5:33)=itmp(1:29) - i4tone(34:37)=icos4b - i4tone(38:66)=itmp(30:58) - i4tone(67:70)=icos4c - i4tone(71:99)=itmp(59:87) - i4tone(100:103)=icos4d + i4tone(5:28)=itmp(1:24) + i4tone(29:32)=icos4b + i4tone(33:56)=itmp(25:48) + i4tone(57:60)=icos4c + i4tone(61:84)=itmp(49:72) + i4tone(85:88)=icos4d + i4tone(89:112)=itmp(73:96) + i4tone(113:116)=icos4e + i4tone(117:140)=itmp(97:120) + i4tone(141:144)=icos4f 999 return end subroutine genft4slow diff --git a/lib/fsk4hf/get_ft4s_bitmetrics.f90 b/lib/fsk4hf/get_ft4s_bitmetrics.f90 new file mode 100644 index 000000000..9d6e3222f --- /dev/null +++ b/lib/fsk4hf/get_ft4s_bitmetrics.f90 @@ -0,0 +1,145 @@ +subroutine get_ft4s_bitmetrics(cd,bitmetrics,badsync) + + include 'ft4s_params.f90' + parameter (NSS=20) + complex cd(0:NN*NSS-1) + complex cs(0:3,NN) + complex csymb(NSS) + integer icos4a(0:3),icos4b(0:3) + integer icos4c(0:3),icos4d(0:3) + integer icos4e(0:3),icos4f(0:3) + integer graymap(0:3) + integer ip(1) + logical one(0:65535,0:15) ! 65536 8-symbol sequences, 16 bits + logical first + logical badsync + real bitmetrics(2*NN,4) + real s2(0:65535) + real s4(0:3,NN) + + data icos4a/0,1,3,2/ + data icos4b/1,0,2,3/ + data icos4c/2,3,1,0/ + data icos4d/3,2,0,1/ + data icos4e/0,2,3,1/ + data icos4f/1,2,0,3/ + data graymap/0,1,3,2/ + data first/.true./ + save first,one + + if(first) then + one=.false. + do i=0,65535 + do j=0,15 + if(iand(i,2**j).ne.0) one(i,j)=.true. + enddo + enddo + first=.false. + endif + + do k=1,NN + i1=(k-1)*NSS + csymb=cd(i1:i1+NSS-1) + call four2a(csymb,NSS,1,-1,1) + cs(0:3,k)=csymb(1:4) + s4(0:3,k)=abs(csymb(1:4)) + enddo + +! Sync quality check + is1=0 + is2=0 + is3=0 + is4=0 + is5=0 + is6=0 + badsync=.false. + ibmax=0 + + do k=1,4 + ip=maxloc(s4(:,k)) + if(icos4a(k-1).eq.(ip(1)-1)) is1=is1+1 + ip=maxloc(s4(:,k+28)) + if(icos4b(k-1).eq.(ip(1)-1)) is2=is2+1 + ip=maxloc(s4(:,k+56)) + if(icos4c(k-1).eq.(ip(1)-1)) is3=is3+1 + ip=maxloc(s4(:,k+84)) + if(icos4d(k-1).eq.(ip(1)-1)) is4=is4+1 + ip=maxloc(s4(:,k+112)) + if(icos4e(k-1).eq.(ip(1)-1)) is5=is5+1 + ip=maxloc(s4(:,k+140)) + if(icos4f(k-1).eq.(ip(1)-1)) is6=is6+1 + enddo + nsync=is1+is2+is3+is4+is5+is6 !Number of correct hard sync symbols, 0-24 + + badsync=.false. +! if(nsync .lt. 8) then +! badsync=.true. +! return +! endif + + do nseq=1,4 !Try coherent sequences of 1, 2, and 4 symbols + if(nseq.eq.1) nsym=1 + if(nseq.eq.2) nsym=2 + if(nseq.eq.3) nsym=4 + if(nseq.eq.4) nsym=8 + nt=2**(2*nsym) + do ks=1,NN-nsym+1,nsym !87+16=103 symbols. + amax=-1.0 + do i=0,nt-1 +! i1=i/64 +! i2=iand(i,63)/16 +! i3=iand(i,15)/4 +! i4=iand(i,3) + i1=i/16384 + i2=iand(i,16383)/4096 + i3=iand(i,4095)/1024 + i4=iand(i,1023)/256 + i5=iand(i,255)/64 + i6=iand(i,63)/16 + i7=iand(i,15)/4 + i8=iand(i,3) + if(nsym.eq.1) then + s2(i)=abs(cs(graymap(i8),ks)) + elseif(nsym.eq.2) then + s2(i)=abs(cs(graymap(i7),ks)+cs(graymap(i8),ks+1)) + elseif(nsym.eq.4) then + s2(i)=abs(cs(graymap(i5),ks ) + & + cs(graymap(i6),ks+1) + & + cs(graymap(i7),ks+2) + & + cs(graymap(i8),ks+3) & + ) + elseif(nsym.eq.8) then + s2(i)=abs(cs(graymap(i1),ks ) + & + cs(graymap(i2),ks+1) + & + cs(graymap(i3),ks+2) + & + cs(graymap(i4),ks+3) + & + cs(graymap(i5),ks+4) + & + cs(graymap(i6),ks+5) + & + cs(graymap(i7),ks+6) + & + cs(graymap(i8),ks+7) & + ) + else + print*,"Error - nsym must be 1, 2, 4, or 8." + endif + enddo + ipt=1+(ks-1)*2 + if(nsym.eq.1) ibmax=1 + if(nsym.eq.2) ibmax=3 + if(nsym.eq.4) ibmax=7 + if(nsym.eq.8) ibmax=15 + do ib=0,ibmax + bm=maxval(s2(0:nt-1),one(0:nt-1,ibmax-ib)) - & + maxval(s2(0:nt-1),.not.one(0:nt-1,ibmax-ib)) + if(ipt+ib.gt.2*NN) cycle + bitmetrics(ipt+ib,nseq)=bm + enddo + enddo + enddo + + call normalizebmet(bitmetrics(:,1),2*NN) + call normalizebmet(bitmetrics(:,2),2*NN) + call normalizebmet(bitmetrics(:,3),2*NN) + call normalizebmet(bitmetrics(:,4),2*NN) + return + +end subroutine get_ft4s_bitmetrics diff --git a/lib/fsk4hf/ldpcsim174_91.f90 b/lib/fsk4hf/ldpcsim174_91.f90 index 637a8a569..a753dd8ad 100644 --- a/lib/fsk4hf/ldpcsim174_91.f90 +++ b/lib/fsk4hf/ldpcsim174_91.f90 @@ -105,6 +105,7 @@ do idb = 20,-10,-1 ! max_iterations is max number of belief propagation iterations call bpdecode174_91(llr, apmask, max_iterations, message77, cw, nharderrors,niterations,ncheck) if( ndepth .ge. 0 .and. nharderrors .lt. 0 ) call osd174_91(llr, apmask, ndepth, message77, cw, nharderrors, dmin) +!call decode174_91(llr, apmask, max_iterations, message77, cw, nharderrors,niterations,ncheck,dmin,nsuper) ! If the decoder finds a valid codeword, nharderrors will be .ge. 0. if( nharderrors .ge. 0 ) then call extractmessage77(message77,msgreceived) diff --git a/lib/fsk4hf/ldpcsim240_101.f90 b/lib/fsk4hf/ldpcsim240_101.f90 new file mode 100644 index 000000000..d87241fa4 --- /dev/null +++ b/lib/fsk4hf/ldpcsim240_101.f90 @@ -0,0 +1,144 @@ +program ldpcsim240_101 + +! End-to-end test of the (240,101)/crc24 encoder and decoders. + + use packjt77 + + parameter(N=240, K=101, M=N-K) + character*8 arg + character*37 msg0,msg + character*77 c77 + character*24 c24 + integer*1 msgbits(101) + integer*1 apmask(240) + integer*1 cw(240) + integer*1 codeword(N),message101(101) + integer ncrc24 + real rxdata(N),llr(N) + real llrd(240) + logical first,unpk77_success + data first/.true./ + + nargs=iargc() + if(nargs.ne.5 .and. nargs.ne.6) then + print*,'Usage: ldpcsim niter ndeep #trials s K [msg]' + print*,'e.g. ldpcsim240_101 20 5 1000 0.85 91 "K9AN K1JT FN20"' + 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 [77,101]' + print*,'WSPR-format message is optional' + 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 + msg0='K9AN K1JT FN20 ' + if(nargs.eq.6) call getarg(6,msg0) + call pack77(msg0,i3,n3,c77) + + 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(c77,'(77i1)') msgbits(1:77) + write(*,*) 'message' + write(*,'(77i1)') msgbits(1:77) + + call get_crc24(msgbits,101,ncrc24) + write(c24,'(b24.24)') ncrc24 + read(c24,'(24i1)') msgbits(78:101) +write(*,'(24i1)') msgbits(78:101) + write(*,*) 'message with crc24' + write(*,'(101i1)') msgbits(1:101) + call encode240_101(msgbits,codeword) + call init_random_seed() + call sgran() + + write(*,*) 'codeword' + write(*,'(77i1,1x,24i1,1x,73i1)') 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 + nberr=nberr+nerr + + rxav=sum(rxdata)/N + rx2av=sum(rxdata*rxdata)/N + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig + 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 bpdecode240_101(llr,apmask,max_iterations,message101,cw,nharderror,niterations,nchecks) + dmin=0.0 + if( (nharderror .lt. 0) .and. (ndeep .ge. 0) ) then +! call osd240_101(llr, Keff, apmask, ndeep, message101, cw, nharderror, dmin) + maxsuper=2 + call decode240_101(llr, Keff, ndeep, apmask, maxsuper, message101, cw, nharderror, iterations, ncheck, dmin, isuper) + 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 + + if(first) then + write(c77,'(77i1)') message101(1:77) +write(*,'(101i1)') message101 + call unpack77(c77,0,msg,unpk77_success) + if(unpk77_success) then + write(*,1100) msg(1:37) +1100 format('Decoded message: ',a37) + else + print*,'Error unpacking message' + endif + first=.false. + endif + enddo + +end program ldpcsim240_101 diff --git a/lib/fsk4hf/osd240_101.f90 b/lib/fsk4hf/osd240_101.f90 new file mode 100644 index 000000000..3e2506805 --- /dev/null +++ b/lib/fsk4hf/osd240_101.f90 @@ -0,0 +1,403 @@ +subroutine osd240_101(llr,k,apmask,ndeep,message101,cw,nhardmin,dmin) +! +! An ordered-statistics decoder for the (240,101) code. +! Message payload is 77 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 77+p1. +! +! Valid values for k are in the range [77,101]. +! + character*24 c24 + integer, parameter:: N=240 + 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 message101(101) + integer indx(N) + real llr(N),rx(N),absrx(N) + + 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=101-k and p1+p2=24. +! +! The last p2 bits of the CRC24 are cascaded with the LDPC code. +! +! The first p1=k-77 CRC24 bits will be used for error detection. +! + allocate( gen(k,N) ) + gen=0 + do i=1,k + message101=0 + message101(i)=1 + if(i.le.77) then + call get_crc24(message101,101,ncrc24) + write(c24,'(b24.24)') ncrc24 + read(c24,'(24i1)') message101(78:101) + message101(78:k)=0 + endif + call encode240_101(message101,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 mrbencode101(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 mrbencode101(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 mrbencode101(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 nextpat101(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 boxit101(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 mrbencode101(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 fetchit101(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 mrbencode101(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 nextpat101(misub,k,nord,iflag) + enddo + endif + +998 continue +! Re-order the codeword to [message bits][parity bits] format. + cw(indices)=cw + hdec(indices)=hdec + message101=cw(1:101) + call get_crc24(message101,101,nbadcrc) + if(nbadcrc.ne.0) nhardmin=-nhardmin + + return +end subroutine osd240_101 + +subroutine mrbencode101(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 mrbencode101 + +subroutine nextpat101(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 nextpat101 + +subroutine boxit101(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 boxit101 + +subroutine fetchit101(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 fetchit101 +