Add routines for ft4slow based on (280,101) code.

This commit is contained in:
Steven Franke 2020-05-02 08:10:35 -05:00
parent e82edf2365
commit b0ef40cb1d
14 changed files with 2295 additions and 5 deletions

View File

@ -604,26 +604,37 @@ set (wsjt_FSRCS
lib/fsk4hf/encode204.f90
lib/fsk4hf/decode174_74.f90
lib/fsk4hf/decode240_101.f90
lib/fsk4hf/decode280_101.f90
lib/fsk4hf/ldpcsim174_91.f90
lib/fsk4hf/ldpcsim174_74.f90
lib/fsk4hf/ldpcsim240_101.f90
lib/fsk4hf/ldpcsim280_101.f90
lib/fsk4hf/get_crc24.f90
lib/fsk4hf/encode174_74.f90
lib/fsk4hf/encode240_101.f90
lib/fsk4hf/encode280_101.f90
lib/fsk4hf/bpdecode174_74.f90
lib/fsk4hf/bpdecode240_101.f90
lib/fsk4hf/bpdecode280_101.f90
lib/fsk4hf/osd174_74.f90
lib/fsk4hf/osd240_101.f90
lib/fsk4hf/osd280_101.f90
lib/fsk4hf/wspr4sim.f90
lib/fsk4hf/genwspr4.f90
lib/fsk4hf/gen_wspr4wave.f90
lib/fsk4hf/wspr4d.f90
lib/fsk4hf/get_wspr4_bitmetrics.f90
lib/fsk4hf/get_ft280_bitmetrics.f90
lib/fsk4hf/ft4slowsim.f90
lib/fsk4hf/ft280sim.f90
lib/fsk4hf/genft4slow.f90
lib/fsk4hf/genft280.f90
lib/fsk4hf/decode240_101.f90
lib/fsk4hf/decode280_101.f90
lib/fsk4hf/ft4sd.f90
lib/fsk4hf/ft280d.f90
lib/fsk4hf/get_ft4s_bitmetrics.f90
lib/fsk4hf/get_ft280_bitmetrics.f90
)
# temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit
@ -1373,6 +1384,9 @@ target_link_libraries (ldpcsim174_74 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 (ldpcsim280_101 lib/fsk4hf/ldpcsim280_101.f90 wsjtx.rc)
target_link_libraries (ldpcsim280_101 wsjt_fort wsjt_cxx)
add_executable (wspr4sim lib/fsk4hf/wspr4sim.f90 wsjtx.rc)
target_link_libraries (wspr4sim wsjt_fort wsjt_cxx)
@ -1382,9 +1396,15 @@ 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 (ft280sim lib/fsk4hf/ft280sim.f90 wsjtx.rc)
target_link_libraries (ft280sim wsjt_fort wsjt_cxx)
add_executable (ft4sd lib/fsk4hf/ft4sd.f90 wsjtx.rc)
target_link_libraries (ft4sd wsjt_fort wsjt_cxx)
add_executable (ft280d lib/fsk4hf/ft280d.f90 wsjtx.rc)
target_link_libraries (ft280d wsjt_fort wsjt_cxx)
endif(WSJT_BUILD_UTILS)
# build the main application

View File

@ -0,0 +1,111 @@
subroutine bpdecode280_101(llr,apmask,maxiterations,message101,cw,nharderror,iter,ncheck)
!
! A log-domain belief propagation decoder for the (280,101) code.
!
integer, parameter:: N=280, 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_280_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 bpdecode280_101

View File

@ -0,0 +1,133 @@
subroutine decode280_101(llr,Keff,ndeep,apmask,maxsuper,message101,cw,nharderror,iter,ncheck,dmin,isuper)
!
! A hybrid bp/osd decoder for the (280,101) code.
!
integer, parameter:: N=280, 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_280_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 osd280_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 decode280_101

View File

@ -0,0 +1,46 @@
subroutine encode280_101(message,codeword)
use, intrinsic :: iso_c_binding
use iso_c_binding, only: c_loc,c_size_t
use crc
integer, parameter:: N=280, 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_280_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 encode280_101

427
lib/fsk4hf/ft280d.f90 Normal file
View File

@ -0,0 +1,427 @@
program ft280d
! Decode ft280 data read from *.c2 or *.wav files.
use packjt77
include 'ft4s280_params.f90'
parameter (NSPS2=NSPS/NDOWN)
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/NDOWN-1) !Complex waveform
complex cframe(0:164*NSPS2-1) !Complex waveform
complex cd(0:164*20-1) !Complex waveform
real*8 fMHz
real llr(280),llra(280),llrb(280),llrc(280),llrd(280)
real candidates(100,2)
real bitmetrics(328,4)
integer ihdr(11)
integer*2 iwave(NMAX) !Generated full-length waveform
integer*1 apmask(280),cw(280)
integer*1 hbits(328)
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)
hmod=1.0
Keff=91
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: ft280d [-a <data_dir>] [-f fMHz] [-h hmod] [-k Keff] 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
call getarg(iarg,arg)
endif
if(arg(1:2).eq.'-f') then
call getarg(iarg+1,arg)
read(arg,*) fMHz
iarg=iarg+2
call getarg(iarg,arg)
endif
if(arg(1:2).eq.'-h') then
call getarg(iarg+1,arg)
read(arg,*) hmod
iarg=iarg+2
call getarg(iarg,arg)
endif
if(arg(1:2).eq.'-k') then
call getarg(iarg+1,arg)
read(arg,*) Keff
iarg=iarg+2
call getarg(iarg,arg)
endif
if(arg(1:2).eq.'-d') then
call getarg(iarg+1,arg)
read(arg,*) ndeep
iarg=iarg+2
endif
ngood=0
ngoodsync=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 ft280_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_ft280(c2,npts,hmod,fs,fa,fb,ncand,candidates) !First approx for freq
del=1.5*hmod*fs/300.0
ndecodes=0
do icand=1,ncand
! do icand=1,1
fc0=candidates(icand,1)
xsnr=candidates(icand,2)
!write(*,*) 'candidates ',icand,fc0,xsnr
do isync=0,1
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_ft280(c2,istart,hmod,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(abs((isbest-429)/429.0) .lt. 0.07 .and. abs(fc2+del).lt.0.2) ngoodsync=ngoodsync+1
!cycle
! if(smax .lt. 100.0 ) cycle
!isbest=429
!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+164*300-1)
call downsample_ft280(cframe,fc2+del,hmod,cd)
s2=sum(cd*conjg(cd))/(20*144)
cd=cd/sqrt(s2)
call get_ft280_bitmetrics(cd,hmod,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( 9: 16).eq.(/0,1,0,0,1,1,1,0/))
ns3=count(hbits(157:164).eq.(/0,0,0,1,1,0,1,1/))
ns4=count(hbits(165:172).eq.(/0,1,0,0,1,1,1,0/))
ns5=count(hbits(313:320).eq.(/0,0,0,1,1,0,1,1/))
ns6=count(hbits(321:328).eq.(/0,1,0,0,1,1,1,0/))
nsync_qual=ns1+ns2+ns3+ns4+ns5+ns6
! if(nsync_qual.lt. 20) cycle
scalefac=2.83
llra( 1:140)=bitmetrics( 17:156, 1)
llra(141:280)=bitmetrics(173:312, 1)
llra=scalefac*llra
llrb( 1:140)=bitmetrics( 17:156, 2)
llrb(141:280)=bitmetrics(173:312, 2)
llrb=scalefac*llrb
llrc( 1:140)=bitmetrics( 17:156, 3)
llrc(141:280)=bitmetrics(173:312, 3)
llrc=scalefac*llrc
llrd( 1:140)=bitmetrics( 17:156, 4)
llrd(141:280)=bitmetrics(173:312, 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 bpdecode280_101(llr,apmask,max_iterations,message101,cw,nhardbp,niterations,nchecks)
! if(nhardbp.lt.0) call osd280_101(llr,Keff,apmask,5,message101,cw,nhardosd,dmin)
maxsuperits=2
if(nhardbp.lt.0) then
! call osd280_101(llr,Keff,apmask,ndeep,message101,cw,nhardosd,dmin)
call decode280_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-2,icand,xsnr,isbest/375.0-1.0,1500.0+fc2+del,msg(1:20),itry,nhardbp,nhardosd,dmin,ijitter
1100 format(i5,2x,i5,2x,f6.1,2x,f6.2,2x,f8.2,2x,a20,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,' ngoodsync: ',ngoodsync
write(*,1120)
1120 format("<DecodeFinished>")
999 end program ft280d
subroutine coherent_sync_ft280(cd0,i0,hmod,f0,itwk,sync)
! Compute sync power for a complex, downsampled FT4s signal.
include 'ft4s280_params.f90'
parameter(NP=NMAX/NDOWN,NSS=NSPS/NDOWN)
complex cd0(0:NP-1)
complex csynca(8*NSS)
complex csync2(8*NSS)
complex ctwk(8*NSS)
complex z1,z2,z3,z4,z5,z6
logical first
integer icos4(0:7)
data icos4/0,1,3,2,1,0,2,3/
data first/.true./
save first,twopi,csynca,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
do i=0,7
dphia=twopi*hmod*icos4(i)/real(NSS)
do j=1,NSS
csynca(k)=cmplx(cos(phia),sin(phia))
phia=mod(phia+dphia,twopi)
k=k+1
enddo
enddo
first=.false.
fac=1.0/(8.0*NSS)
endif
i1=i0 !four Costas arrays
i2=i0+78*NSS
i3=i0+156*NSS
z1=0.
z2=0.
z3=0.
if(itwk.eq.1) then
dt=1/(12000.0/32.0)
dphi=twopi*f0*dt
phi=0.0
do i=1,8*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+8*NSS-1.le.NP-1) then
z1=sum(cd0(i1:i1+8*NSS-1)*conjg(csync2))
! z1=abs(sum(cd0(i1:i1+4*NSS-1)*conjg(csync2(1:4*NSS))))**2
! z1=z1+abs(sum(cd0(i1+4*NSS:i1+8*NSS-1)*conjg(csync2(4*NSS+1:8*NSS))))**2
elseif( i1.lt.0 ) then
npts=(i1+8*NSS-1)/2
if(npts.le.40) then
z1=0.
else
z1=sum(cd0(0:i1+8*NSS-1)*conjg(csync2(8*NSS-npts:)))
endif
endif
if(i2.ge.0 .and. i2+8*NSS-1.le.NP-1) then
z2=sum(cd0(i2:i2+8*NSS-1)*conjg(csync2))
! z2=abs(sum(cd0(i2:i2+4*NSS-1)*conjg(csync2(1:4*NSS))))**2
! z2=z2+abs(sum(cd0(i2+4*NSS:i2+8*NSS-1)*conjg(csync2(4*NSS+1:8*NSS))))**2
endif
if(i3.ge.0 .and. i3+8*NSS-1.le.NP-1) then
z3=sum(cd0(i3:i3+8*NSS-1)*conjg(csync2))
! z3=abs(sum(cd0(i3:i3+4*NSS-1)*conjg(csync2(1:4*NSS))))**2
! z3=z3+abs(sum(cd0(i3+4*NSS:i3+8*NSS-1)*conjg(csync2(4*NSS+1:8*NSS))))**2
elseif( i3+8*NSS-1.gt.NP-1 ) then
npts=(NP-1-i3+1)
if(npts.le.40) then
z3=0.
else
z3=sum(cd0(i3:i3+npts-1)*conjg(csync2(1:npts)))
endif
endif
sync = p(z1) + p(z2) + p(z3)
!sync=z1+z2+z3
return
end subroutine coherent_sync_ft280
subroutine downsample_ft280(ci,f0,hmod,co)
parameter(NI=164*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/28.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*hmod
b=10.0*hmod
icutoff=nint(24.0/df)
do i=1,NO/2
! arg=(i*df/b)**2
! filt=exp(-arg)
filt=0
if(i.le.icutoff) filt=1
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_ft280
subroutine getcandidate_ft280(c,npts,hmod,fs,fa,fb,ncand,candidates)
parameter(NFFT1=120*12000/28,NH1=NFFT1/2,NFFT2=120*12000/300,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)/(hmod*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)-26.5
endif
enddo
return
end subroutine getcandidate_ft280
subroutine ft280_downsample(iwave,c)
! Input: i*2 data in iwave() at sample rate 12000 Hz
! Output: Complex data in c(), sampled at 375 Hz
include 'ft4s280_params.f90'
parameter (NFFT2=NMAX/28)
integer*2 iwave(NMAX)
complex c(0:NMAX/28-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/28-1)
return
end subroutine ft280_downsample

113
lib/fsk4hf/ft280sim.f90 Normal file
View File

@ -0,0 +1,113 @@
program ft280sim
! Generate simulated signals for experimental slow FT4 mode
use wavhdr
use packjt77
include 'ft4s280_params.f90' !Set various constants
type(hdr) h !Header for .wav file
character arg*12,fname*17
character msg37*37,msgsent37*37
character c77*77
complex c0(0:NMAX-1)
complex c(0:NMAX-1)
real wave(NMAX)
integer itone(NN)
integer*1 msgbits(101)
integer*2 iwave(NMAX) !Generated full-length waveform
! Get command-line argument(s)
nargs=iargc()
if(nargs.ne.8) then
print*,'Usage: ft280sim "message" f0 DT h fdop del nfiles snr'
print*,'Examples: ft280sim "K1JT K9AN EN50" 1500 0.0 1.0 0.1 1.0 10 -15'
go to 999
endif
call getarg(1,msg37) !Message to be transmitted
call getarg(2,arg)
read(arg,*) f0 !Frequency (only used for single-signal)
call getarg(3,arg)
read(arg,*) xdt !Time offset from nominal (s)
call getarg(4,arg)
read(arg,*) hmod !Modulation index, h
call getarg(5,arg)
read(arg,*) fspread !Watterson frequency spread (Hz)
call getarg(6,arg)
read(arg,*) delay !Watterson delay (ms)
call getarg(7,arg)
read(arg,*) nfiles !Number of files
call getarg(8,arg)
read(arg,*) snrdb !SNR_2500
nfiles=abs(nfiles)
twopi=8.0*atan(1.0)
fs=12000.0 !Sample rate (Hz)
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of symbols (s)
baud=1.0/tt !Keying rate (baud)
txt=NZ2*dt !Transmission length (s)
bandwidth_ratio=2500.0/(fs/2.0)
sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
call genft280(msg37,0,msgsent37,msgbits,itone)
write(*,*)
write(*,'(a9,a37,3x,a7,i1,a1,i1)') 'Message: ',msgsent37,'i3.n3: ',i3,'.',n3
write(*,1000) f0,xdt,hmod,txt,snrdb
1000 format('f0:',f9.3,' DT:',f6.2,' hmod:',f6.3,' TxT:',f6.1,' SNR:',f6.1)
write(*,*)
if(i3.eq.1) then
write(*,*) ' mycall hiscall hisgrid'
write(*,'(28i1,1x,i1,1x,28i1,1x,i1,1x,i1,1x,15i1,1x,3i1)') msgbits(1:77)
else
write(*,'(a14)') 'Message bits: '
write(*,'(50i1,1x,24i1)') msgbits
endif
write(*,*)
write(*,'(a17)') 'Channel symbols: '
write(*,'(10i1)') itone
write(*,*)
call sgran()
fsample=12000.0
icmplx=1
call gen_wspr4wave(itone,NN,NSPS,fsample,hmod,f0,c0,wave,icmplx,NMAX)
k=nint((xdt+1.0)/dt)-NSPS
c0=cshift(c0,-k)
if(k.gt.0) c0(0:k-1)=0.0
if(k.lt.0) c0(NMAX+k:NMAX-1)=0.0
do ifile=1,nfiles
c=c0
if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NMAX,NZ,fs,delay,fspread)
c=sig*c
wave=real(c)
if(snrdb.lt.90) then
do i=1,NMAX !Add gaussian noise at specified SNR
xnoise=gran()
wave(i)=wave(i) + xnoise
enddo
endif
gain=100.0
if(snrdb.lt.90.0) then
wave=gain*wave
else
datpk=maxval(abs(wave))
fac=32766.9/datpk
wave=fac*wave
endif
if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped."
iwave=nint(wave)
h=default_header(12000,NMAX)
write(fname,1102) ifile
1102 format('000000_',i6.6,'.wav')
open(10,file=fname,status='unknown',access='stream')
write(10) h,iwave !Save to *.wav file
close(10)
write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a17)
enddo
999 end program ft280sim

View File

@ -0,0 +1,16 @@
! FT4S280
! LDPC(280,101)/CRC24 code, six 4x4 Costas arrays for sync, ramp-up and ramp-down symbols
parameter (KK=77) !Information bits (77 + CRC24)
parameter (ND=140) !Data symbols
parameter (NS=24) !Sync symbols
parameter (NN=NS+ND) !Sync and data symbols (164)
parameter (NN2=NS+ND+2) !Total channel symbols (166)
parameter (NSPS=8400) !Samples per symbol at 12000 S/s
parameter (NZ=NSPS*NN) !Sync and Data samples (1,377,600)
parameter (NZ2=NSPS*NN2) !Total samples in shaped waveform (1,394,400)
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=28) !Downsample factor

View File

@ -30,10 +30,11 @@ program ft4sd
tt=NSPS*dt !Duration of "itone" symbols (s)
txt=NZ*dt !Transmission length (s)
hmod=1.0
Keff=91
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: ft4sd [-a <data_dir>] [-f fMHz] file1 [file2 ...]'
print*,'Usage: ft4sd [-a <data_dir>] [-f fMHz] [-h hmod] [-k Keff] file1 [file2 ...]'
go to 999
endif
iarg=1
@ -42,17 +43,24 @@ program ft4sd
if(arg(1:2).eq.'-a') then
call getarg(iarg+1,data_dir)
iarg=iarg+2
call getarg(iarg,arg)
endif
call getarg(iarg,arg)
if(arg(1:2).eq.'-f') then
call getarg(iarg+1,arg)
read(arg,*) fMHz
iarg=iarg+2
call getarg(iarg,arg)
endif
if(arg(1:2).eq.'-h') then
call getarg(iarg+1,arg)
read(arg,*) hmod
iarg=iarg+2
call getarg(iarg,arg)
endif
if(arg(1:2).eq.'-k') then
call getarg(iarg+1,arg)
read(arg,*) Keff
iarg=iarg+2
endif
ngood=0
@ -86,7 +94,6 @@ program ft4sd
del=1.5*hmod*fs/300.0
ndecodes=0
do icand=1,ncand
! do icand=1,1
fc0=candidates(icand,1)
xsnr=candidates(icand,2)
!write(*,*) 'candidates ',icand,fc0,xsnr
@ -124,7 +131,7 @@ program ft4sd
! if(smax .lt. 100.0 ) cycle
!isbest=375
!fc2=0
!fc2=-del
do ijitter=0,2
if(ijitter.eq.0) ioffset=0
if(ijitter.eq.1) ioffset=45
@ -185,10 +192,10 @@ program ft4sd
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(Keff.eq.77) ndeep=4
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)

95
lib/fsk4hf/genft280.f90 Normal file
View File

@ -0,0 +1,95 @@
subroutine genft280(msg0,ichk,msgsent,msgbits,i4tone)
! Encode an FT4 message
! Input:
! - msg0 requested message to be transmitted
! - ichk if ichk=1, return only msgsent
! - msgsent message as it will be decoded
! - i4tone array of audio tone values, {0,1,2,3}
! Frame structure:
! s4s4 d70 s4s4 d70 s4s4
! Message duration: TxT = 144*9600/12000 = 115.2 s
use packjt77
include 'ft4s280_params.f90'
character*37 msg0
character*37 message !Message to be generated
character*37 msgsent !Message as it will be received
character*77 c77
character*24 c24
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),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/
message=msg0
do i=1, 37
if(ichar(message(i:i)).eq.0) then
message(i:37)=' '
exit
endif
enddo
do i=1,37 !Strip leading blanks
if(message(1:1).ne.' ') exit
message=message(i+1:)
enddo
i3=-1
n3=-1
call pack77(message,i3,n3,c77)
call unpack77(c77,0,msgsent,unpk77_success) !Unpack to get msgsent
msgbits=0
read(c77,'(77i1)') msgbits(1:77)
call get_crc24(msgbits,101,ncrc24)
write(c24,'(b24.24)') ncrc24
read(c24,'(24i1)') msgbits(78:101)
if(ichk.eq.1) go to 999
if(unpk77_success) go to 2
1 msgbits=0
itone=0
msgsent='*** bad message *** '
go to 999
entry get_ft4s280_tones_from_101bits(msgbits,i4tone)
2 call encode280_101(msgbits,codeword)
! Grayscale mapping:
! bits tone
! 00 0
! 01 1
! 11 2
! 10 3
do i=1,ND
is=codeword(2*i)+2*codeword(2*i-1)
if(is.le.1) itmp(i)=is
if(is.eq.2) itmp(i)=3
if(is.eq.3) itmp(i)=2
enddo
i4tone(1:4)=icos4a
i4tone(5:8)=icos4b
i4tone(9:78)=itmp(1:70)
i4tone(79:82)=icos4a
i4tone(83:86)=icos4b
i4tone(87:156)=itmp(71:140)
i4tone(157:160)=icos4a
i4tone(161:164)=icos4b
999 return
end subroutine genft280

View File

@ -0,0 +1,117 @@
subroutine get_ft280_bitmetrics(cd,hmod,bitmetrics,badsync)
include 'ft4s280_params.f90'
parameter (NSS=20)
complex cd(0:NN*NSS-1)
complex cs(0:3,NN)
complex csymb(NSS)
complex c1(NSS,0:3) ! ideal waveforms, 20 samples per symbol, 4 tones
complex ccor(0:3,NN) ! correlations with each ideal waveform, for each symbol
complex cp(0:3) ! accumulated phase shift over symbol types 0:3
complex csum,cterm
integer icos8(0:7)
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 icos8/0,1,3,2,1,0,2,3/
data graymap/0,1,3,2/
data first/.true./
save first,one,c1,cp
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
twopi=8.0*atan(1.0)
dphi=twopi*hmod/NSS
do itone=0,3
dp=(itone-1.5)*dphi
phi=0.0
do j=1,NSS
c1(j,itone)=cmplx(cos(phi),sin(phi))
phi=mod(phi+dp,twopi)
enddo
cp(itone)=cmplx(cos(phi),sin(phi))
enddo
first=.false.
endif
do k=1,NN
i1=(k-1)*NSS
csymb=cd(i1:i1+NSS-1)
do itone=0,3
cs(itone,k)=sum(csymb*conjg(c1(:,itone)))
enddo
s4(0:3,k)=abs(cs(0:3,k))
enddo
! Sync quality check
is1=0
is2=0
is3=0
badsync=.false.
ibmax=0
do k=1,8
ip=maxloc(s4(:,k))
if(icos8(k-1).eq.(ip(1)-1)) is1=is1+1
ip=maxloc(s4(:,k+78))
if(icos8(k-1).eq.(ip(1)-1)) is2=is2+1
ip=maxloc(s4(:,k+156))
if(icos8(k-1).eq.(ip(1)-1)) is3=is3+1
enddo
nsync=is1+is2+is3 !Number of correct hard sync symbols, 0-24
badsync=.false.
! if(nsync .lt. 8) then
! badsync=.true.
! return
! endif
do nseq=4,1,-1 !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=4**nsym
do ks=1,NN-nsym+1,nsym
s2=0
do i=0,nt-1
csum=0
cterm=1
do j=0,nsym-1
ntone=mod(i/4**(nsym-1-j),4)
csum=csum+cs(graymap(ntone),ks+j)*cterm
cterm=cterm*conjg(cp(graymap(ntone)))
enddo
s2(i)=abs(csum)
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_ft280_bitmetrics

View File

@ -0,0 +1,182 @@
character*26 g(179)
data g/ &
"c919bcbfe4091279702a761e98", &
"51b952dddd36200cf73cc1ed30", &
"15871d32e8e888439180cf6fd8", &
"581f858f6c89ee5ccb91664358", &
"3515e85cedf905eda366a8fc20", &
"e9fcaa6aaa9bab21bc91174e80", &
"0ac73221d424e8747628b13968", &
"4999f7116446f1a7a7a1453a30", &
"0e92773bff2a6d4f09caa48898", &
"7dfaec97c17679f6c3b6a425f0", &
"00707d76a2a7d90297ee39f660", &
"8048cc93fc4ad84ccfc021e6e0", &
"0c13df64062fed419c9bf43400", &
"5523d84459c826b7bc3335d508", &
"828ee2552144d041ed44ada8e0", &
"3f1b89fbd93f674df4813f0898", &
"4e13df64062fed419c9bf43400", &
"5d8645307d3d442991d6efafd0", &
"e5cd9b98d73aab17ce04c4df10", &
"06d26e11e2d02e9cb4f191c2b0", &
"5630cebc5b3a09f7d4fe58fab0", &
"bbfa9591589229738ecbc19288", &
"c98654d1f1f16d507e9bb77cf0", &
"c2af2107bb2bdff49d909dc730", &
"51da7da0c9b1bd18a15f580068", &
"5bdfd83e7ca3097146a5359428", &
"34fc4d397d97ca3ceb272f49a0", &
"6716a6d027ade94010e9aa90b0", &
"62ac7bb089d1a13f6e89f92348", &
"737c3ab63210e195e92e8ad478", &
"db2da5b8a21d22a7122ad80e60", &
"1226525dba4221d4768a495878", &
"a99deb4c9b7d316917b1ece958", &
"8123fb46556f22a0b57bdc7eb0", &
"cc6a80e87a7a9bf8addb17a6a8", &
"3d42bb6ca1c8d30e6cee77aa10", &
"ad15a0c2f36d4409a458cc83c0", &
"766a04039736bd8be23513ae58", &
"257a3da77558d7c707170c30c8", &
"8e54a55fd9f00eb669ab787678", &
"4ef1a73cc9da8670d83bebc588", &
"be8bb82558d44fea1ab27376a0", &
"ea9db4f88c60edf410cb0128d8", &
"a84e19a5261818262ee7247278", &
"51f99e4ea17cf84038d4e00bd0", &
"610560db4095fc44d2465308a0", &
"7688745b59c3d6baa6950c4f50", &
"4b8794914d365b6802bd62a9c8", &
"f62c211d05ed28802b9d278298", &
"b9cd45b2ffa8c0dd688f8d2bc0", &
"68555e81f4227a48e76878bc98", &
"7ab58f11d41a2d38b80d2a7558", &
"aba2d33e69077b6acad393af68", &
"81e5f58fa3ab563e73706201a8", &
"7586aea816750c41671eaa7db8", &
"af37b0a97ba5334a3dd01948e8", &
"4fdc03c263a0c42dcc265f7dc8", &
"b23f2d7f28748944cdfffd5af0", &
"5c1e6f37dfba8feacaafdb0f78", &
"3a85b051f4f1c930d921f60828", &
"72319352bd8022ce2cae2e7858", &
"78b79f633ac6879e3ac3a005a0", &
"9f0c470609669953b23328de60", &
"86d3745d50142c82a066ab9490", &
"743e7bf411490f36a9799e37e8", &
"9b8378677870933ef360d7e418", &
"5f7adbf515b663a1434b0d47d8", &
"13249a96b14c6cdcfae5009eb0", &
"da9570e0a52125d0dc4dec4430", &
"ada13ce2dbcb57e2f5b31172f0", &
"84f5485886d4157e9d37efb4d0", &
"23f58c3200bab4ae5dee54edd0", &
"d4377aadf8acb19d4369613ac8", &
"17cefcf65c87885fb6c4d537a0", &
"59d70b8536488298930aaea7f8", &
"49e8dbb08c2ecdaa84bb6a5378", &
"e1694479ecc1f87e503f959e50", &
"dbb3fc94f0f70d4bd4dcf302d8", &
"4ccb3a56f80c236424683b1588", &
"f4f123c72596a00397d56fcdf8", &
"13f9cf266a6957b87cd2b576f0", &
"0904d341bc0878460cd8361ac0", &
"69fd534caf2cccf9c90659a038", &
"757d3d95089a5bc20a7b77c618", &
"30df1d7b8124415c73190b08d8", &
"d39319584987dce0c44176d5d8", &
"1a81df299eb7434a5b6b9322a0", &
"fe4acfab1c22c7bea222f1a6b0", &
"2f2fde37fa8f87a318f7bcda10", &
"fae712210c23665aa7a3f10620", &
"977b8407c7fd22d7715077ee78", &
"2ab2b355b3477df0b738c49d48", &
"93a2468cfd11a522b310069d88", &
"0e5ae6e789ded3c0d436359318", &
"9ece0b13a3c06d560a15d3c448", &
"838e8bbf5e671503ea72ba3118", &
"7c827de9a87d740763c69c6778", &
"1fe395e4e2e6d1373602243488", &
"f2c4efee3d0ce2e22749be9e20", &
"46405cca0e40f36ab83de4a998", &
"8b6e931355a79630ef2dbdbdb8", &
"10df1d3b8124415c72190b08d8", &
"cdff258b07a4f7cfe5c2210ba8", &
"1515e85cedf904eda366a8fc20", &
"a38276f2d077abc1da5e177868", &
"67a7b5ab66f21f391d306c3330", &
"29492cc630f9bad1fdedf0c990", &
"490a6dd38170eab178f7cebf78", &
"ca9db4e88c60edf410cf0128d8", &
"e3f1c23fa8531fb1e4c7768d88", &
"39d7d8fbbb689b2a9bedfd4dd0", &
"d1b952dd5d36200cf734c1ed30", &
"0820a5ccb970d1ab109d84d700", &
"58bc3c509fcd7874e9b1533ba8", &
"08ed7724ac66b7974499b12f40", &
"4738529b2fd04afd89184b64b8", &
"7155b496e3b9f687135b4c55b8", &
"b5d1d3cf38b1765dd730d8b960", &
"296de2c373773a869b9cf804c8", &
"1cdf18b99bcc47ae72bf59df68", &
"ad0888db89dd794be0b2660e98", &
"1f2a8db9db19cd4d69a735d930", &
"44b720007480382206fdbfbb18", &
"c63817aad3801fb993ea9032c0", &
"d44707db5a0b489fd48748cca8", &
"49f98a67c6e128a5300f7ccc50", &
"04849fa9da91d4514355406388", &
"dfad3a11788cf6f6517f987de8", &
"47078a19e38a0763cabd7c8d70", &
"aafa7f864f0da5bc78f8e57ba8", &
"8acb5a34e18e111023b3e7b1f8", &
"5acc41263d6aa1767e5e6acdc8", &
"27623a9f6c1174e35394191820", &
"1f2bde9c006b3b687964b1c5e0", &
"b01c6e357bce202244b4a88d08", &
"61c85d74d7e97576507c9b0e88", &
"bcad5a44d75ae40bc43559d268", &
"10584eaf319552194418563de0", &
"b29b011d717d10a22de0983980", &
"2f9b42d7d2299449491c612b20", &
"389ba33f5fec3bfb3a0ef86b50", &
"3df89f78c19fb27ae7ff19d360", &
"65ff6ba4e107aa919a6afb4ff0", &
"39b607c3f09679a62e134cd390", &
"94ad06f7b7414727d92f998930", &
"169200459898ae0bc7f06714a0", &
"c7a5a945adebb554cb4d86a830", &
"f37c3ab63230e195e92e8ad478", &
"559a51262e91aa9ba0fa96af48", &
"fb2998ca916a557463d00fb160", &
"aa32462ada57a76ae132fc8de8", &
"e6df6b19f58bfee0b96b731b90", &
"e984335d40a54fe914a6249110", &
"ea73d8f3f14bd9fe2374e39120", &
"3adab8e51c36f53584e3669c88", &
"74ef69f64dc4fef86c3b1fe640", &
"d01c6bc112d7ae3e4ba4820a78", &
"62923979fd3c3d1153bcaaf338", &
"038f72995b5072df8fe5f4dfa0", &
"9f07e7cea2f1476fb035978790", &
"2a5aad6a75d5c86cab38fd0070", &
"a254a09cc3180854688d2aa9c8", &
"0495639712a04820f7038ae7c0", &
"d99fc716ca825ad45cca8f4518", &
"01b8d558073c0377ce67344a50", &
"2fbd0f86a17c3f93713fbd09a0", &
"c29dc84bec7b4cd00dd1c17380", &
"5e6238b823f530ae017a03f0e0", &
"51203d329c68b061977d78d4c0", &
"1186729e08cf1dfbec30237968", &
"40363018b431224a1f559d2908", &
"e334e78442b614a0c9a377e1b8", &
"ff2eda86339f589f96382f52e0", &
"58a30e07fc7a37a4f858623778", &
"f5067fe407a4c3b94ce7b63e48", &
"1d09ced788a3642bc0ec640ec8", &
"17734ca67d53cd9d8595970668", &
"47953c2105bd94bff079672740", &
"3444682d1dc0ab486036c1b0d0"/

View File

@ -0,0 +1,476 @@
data Mn/ &
150, 151, 161, &
6, 164, 172, &
92, 128, 158, &
2, 63, 135, &
3, 14, 22, &
4, 18, 29, &
5, 17, 164, &
7, 99, 179, &
8, 88, 115, &
9, 62, 110, &
10, 107, 154, &
11, 50, 140, &
12, 28, 33, &
13, 31, 170, &
15, 69, 175, &
16, 77, 178, &
19, 70, 91, &
20, 95, 177, &
21, 96, 106, &
23, 129, 168, &
24, 49, 169, &
25, 65, 102, &
26, 82, 171, &
27, 45, 137, &
30, 89, 119, &
32, 148, 158, &
34, 94, 152, &
35, 44, 92, &
36, 39, 138, &
37, 55, 58, &
38, 121, 165, &
40, 81, 162, &
41, 139, 150, &
42, 43, 83, &
46, 80, 114, &
47, 52, 54, &
48, 166, 173, &
38, 53, 87, &
56, 64, 126, &
57, 67, 127, &
59, 156, 159, &
60, 97, 133, &
61, 118, 161, &
66, 100, 123, &
68, 124, 131, &
71, 101, 155, &
72, 74, 144, &
73, 112, 141, &
75, 136, 149, &
59, 78, 117, &
79, 130, 163, &
84, 93, 113, &
86, 108, 163, &
103, 146, 157, &
70, 104, 145, &
105, 128, 142, &
74, 109, 122, &
54, 111, 153, &
116, 154, 176, &
120, 132, 167, &
21, 125, 147, &
134, 143, 166, &
7, 81, 160, &
32, 99, 174, &
1, 93, 104, &
2, 69, 98, &
3, 33, 152, &
4, 46, 159, &
5, 126, 178, &
6, 127, 147, &
8, 101, 110, &
9, 73, 158, &
10, 120, 123, &
11, 122, 125, &
12, 58, 170, &
13, 88, 105, &
14, 133, 150, &
15, 92, 100, &
16, 90, 108, &
17, 44, 106, &
18, 35, 175, &
19, 94, 179, &
20, 97, 153, &
22, 109, 130, &
23, 63, 140, &
24, 37, 146, &
25, 141, 168, &
26, 95, 115, &
27, 107, 149, &
28, 91, 168, &
29, 134, 144, &
30, 31, 169, &
34, 40, 96, &
36, 156, 172, &
39, 61, 135, &
41, 42, 121, &
43, 57, 117, &
45, 62, 72, &
47, 137, 167, &
48, 83, 116, &
49, 65, 173, &
1, 50, 141, &
2, 8, 150, &
3, 62, 140, &
4, 104, 124, &
5, 128, 139, &
6, 64, 159, &
7, 103, 176, &
2, 11, 104, &
9, 71, 85, &
10, 80, 131, &
11, 17, 130, &
12, 148, 156, &
13, 39, 164, &
14, 15, 167, &
14, 32, 89, &
16, 114, 135, &
8, 164, 169, &
18, 107, 129, &
19, 53, 102, &
20, 134, 170, &
21, 43, 145, &
22, 24, 76, &
23, 44, 146, &
19, 22, 101, &
25, 41, 48, &
26, 46, 58, &
27, 82, 87, &
28, 78, 179, &
29, 73, 81, &
30, 116, 161, &
31, 96, 157, &
15, 58, 172, &
10, 33, 160, &
34, 110, 118, &
33, 35, 113, &
36, 166, 175, &
32, 37, 152, &
38, 57, 74, &
13, 82, 176, &
40, 42, 45, &
25, 57, 177, &
40, 120, 136, &
21, 92, 121, &
23, 34, 147, &
12, 45, 54, &
3, 46, 48, &
47, 91, 169, &
26, 61, 132, &
49, 123, 147, &
1, 79, 88, &
51, 97, 101, &
52, 155, 177, &
24, 72, 105, &
54, 84, 106, &
55, 63, 126, &
56, 72, 163, &
38, 63, 170, &
37, 71, 178, &
20, 49, 59, &
30, 60, 117, &
61, 65, 137, &
41, 98, 119, &
47, 51, 62, &
6, 76, 131, &
55, 70, 81, &
66, 111, 119, &
60, 67, 94, &
68, 112, 132, &
9, 69, 157, &
70, 75, 89, &
69, 108, 153, &
44, 53, 77, &
29, 130, 149, &
65, 103, 125, &
74, 85, 156, &
56, 67, 68, &
77, 138, 144, &
28, 95, 138, &
79, 133, 142, &
35, 50, 86, &
73, 78, 137, &
27, 126, 175, &
83, 100, 143, &
42, 142, 168, &
40, 48, 158, &
86, 95, 174, &
39, 109, 129, &
59, 88, 125, &
6, 89, 155, &
36, 90, 102, &
75, 97, 141, &
43, 146, 148, &
93, 149, 168, &
52, 83, 94, &
80, 87, 106, &
91, 96, 143, &
3, 43, 126, &
98, 154, 162, &
99, 115, 173, &
5, 84, 100, &
64, 133, 154, &
90, 117, 158, &
7, 108, 151, &
4, 128, 167, &
105, 127, 136, &
1, 83, 114, &
107, 127, 134, &
4, 108, 170, &
92, 109, 171, &
110, 113, 122, &
111, 124, 166, &
12, 112, 150, &
2, 95, 105, &
17, 114, 118, &
99, 139, 144, &
116, 165, 178, &
5, 22, 73, &
16, 115, 162, &
13, 34, 41, &
120, 122, 151, &
121, 160, 172, &
8, 37, 102, &
123, 140, 165, &
7, 53, 93, &
9, 10, 130, &
11, 30, 58, &
31, 66, 179, &
14, 31, 45, &
15, 88, 129, &
18, 101, 148, &
16, 62, 127, &
17, 20, 68, &
19, 86, 98, &
25, 106, 163, &
135, 152, 163, &
23, 124, 137, &
21, 28, 71, &
24, 26, 153, &
29, 90, 123, &
32, 113, 134, &
35, 57, 169, &
27, 50, 139, &
33, 60, 65, &
38, 61, 142, &
145, 153, 154, &
39, 67, 81, &
36, 84, 133, &
18, 161, 173, &
93, 155, 171, &
42, 99, 131, &
49, 87, 162, &
51, 56, 168, &
47, 125, 144, &
44, 143, 159, &
46, 75, 138, &
52, 78, 107, &
54, 109, 174, &
64, 110, 179, &
159, 165, 174, &
66, 135, 171, &
63, 76, 117, &
59, 111, 120, &
72, 160, 166, &
70, 118, 156, &
55, 157, 173, &
74, 100, 176, &
77, 112, 145, &
69, 141, 147, &
94, 140, 151, &
51, 82, 104, &
85, 98, 167, &
80, 119, 146, &
97, 122, 172, &
90, 96, 132, &
79, 91, 178, &
103, 136, 152, &
1, 76, 85, &
115, 121, 149, &
116, 175, 177/
data Nm/ &
65, 102, 151, 207, 278, 0, &
4, 66, 103, 109, 214, 0, &
5, 67, 104, 147, 198, 0, &
6, 68, 105, 205, 209, 0, &
7, 69, 106, 201, 218, 0, &
2, 70, 107, 165, 190, 0, &
8, 63, 108, 204, 225, 0, &
9, 71, 103, 118, 223, 0, &
10, 72, 110, 170, 226, 0, &
11, 73, 111, 134, 226, 0, &
12, 74, 109, 112, 227, 0, &
13, 75, 113, 146, 213, 0, &
14, 76, 114, 140, 220, 0, &
5, 77, 115, 116, 229, 0, &
15, 78, 115, 133, 230, 0, &
16, 79, 117, 219, 232, 0, &
7, 80, 112, 215, 233, 0, &
6, 81, 119, 231, 249, 0, &
17, 82, 120, 125, 234, 0, &
18, 83, 121, 160, 233, 0, &
19, 61, 122, 144, 238, 0, &
5, 84, 123, 125, 218, 0, &
20, 85, 124, 145, 237, 0, &
21, 86, 123, 154, 239, 0, &
22, 87, 126, 142, 235, 0, &
23, 88, 127, 149, 239, 0, &
24, 89, 128, 183, 243, 0, &
13, 90, 129, 179, 238, 0, &
6, 91, 130, 174, 240, 0, &
25, 92, 131, 161, 227, 0, &
14, 92, 132, 228, 229, 0, &
26, 64, 116, 138, 241, 0, &
13, 67, 134, 136, 244, 0, &
27, 93, 135, 145, 220, 0, &
28, 81, 136, 181, 242, 0, &
29, 94, 137, 191, 248, 0, &
30, 86, 138, 159, 223, 0, &
31, 38, 139, 158, 245, 0, &
29, 95, 114, 188, 247, 0, &
32, 93, 141, 143, 186, 0, &
33, 96, 126, 163, 220, 0, &
34, 96, 141, 185, 251, 0, &
34, 97, 122, 193, 198, 0, &
28, 80, 124, 173, 255, 0, &
24, 98, 141, 146, 229, 0, &
35, 68, 127, 147, 256, 0, &
36, 99, 148, 164, 254, 0, &
37, 100, 126, 147, 186, 0, &
21, 101, 150, 160, 252, 0, &
12, 102, 181, 243, 0, 0, &
152, 164, 253, 271, 0, 0, &
36, 153, 195, 257, 0, 0, &
38, 120, 173, 225, 0, 0, &
36, 58, 146, 155, 258, 0, &
30, 156, 166, 266, 0, 0, &
39, 157, 177, 253, 0, 0, &
40, 97, 139, 142, 242, 0, &
30, 75, 127, 133, 227, 0, &
41, 50, 160, 189, 263, 0, &
42, 161, 168, 244, 0, 0, &
43, 95, 149, 162, 245, 0, &
10, 98, 104, 164, 232, 0, &
4, 85, 156, 158, 262, 0, &
39, 107, 202, 259, 0, 0, &
22, 101, 162, 175, 244, 0, &
44, 167, 228, 261, 0, 0, &
40, 168, 177, 247, 0, 0, &
45, 169, 177, 233, 0, 0, &
15, 66, 170, 172, 269, 0, &
17, 55, 166, 171, 265, 0, &
46, 110, 159, 238, 0, 0, &
47, 98, 154, 157, 264, 0, &
48, 72, 130, 182, 218, 0, &
47, 57, 139, 176, 267, 0, &
49, 171, 192, 256, 0, 0, &
123, 165, 262, 278, 0, 0, &
16, 173, 178, 268, 0, 0, &
50, 129, 182, 257, 0, 0, &
51, 151, 180, 276, 0, 0, &
35, 111, 196, 273, 0, 0, &
32, 63, 130, 166, 247, 0, &
23, 128, 140, 271, 0, 0, &
34, 100, 184, 195, 207, 0, &
52, 155, 201, 248, 0, 0, &
110, 176, 272, 278, 0, 0, &
53, 181, 187, 234, 0, 0, &
38, 128, 196, 252, 0, 0, &
9, 76, 151, 189, 230, 0, &
25, 116, 171, 190, 0, 0, &
79, 191, 203, 240, 275, 0, &
17, 90, 148, 197, 276, 0, &
3, 28, 78, 144, 210, 0, &
52, 65, 194, 225, 250, 0, &
27, 82, 168, 195, 270, 0, &
18, 88, 179, 187, 214, 0, &
19, 93, 132, 197, 275, 0, &
42, 83, 152, 192, 274, 0, &
66, 163, 199, 234, 272, 0, &
8, 64, 200, 216, 251, 0, &
44, 78, 184, 201, 267, 0, &
46, 71, 125, 152, 231, 0, &
22, 120, 191, 223, 0, 0, &
54, 108, 175, 277, 0, 0, &
55, 65, 105, 109, 271, 0, &
56, 76, 154, 206, 214, 0, &
19, 80, 155, 196, 235, 0, &
11, 89, 119, 208, 257, 0, &
53, 79, 172, 204, 209, 0, &
57, 84, 188, 210, 258, 0, &
10, 71, 135, 211, 259, 0, &
58, 167, 212, 263, 0, 0, &
48, 169, 213, 268, 0, 0, &
52, 136, 211, 241, 0, 0, &
35, 117, 207, 215, 0, 0, &
9, 88, 200, 219, 279, 0, &
59, 100, 131, 217, 280, 0, &
50, 97, 161, 203, 262, 0, &
43, 135, 215, 265, 0, 0, &
25, 163, 167, 273, 0, 0, &
60, 73, 143, 221, 263, 0, &
31, 96, 144, 222, 279, 0, &
57, 74, 211, 221, 274, 0, &
44, 73, 150, 224, 240, 0, &
45, 105, 212, 237, 0, 0, &
61, 74, 175, 189, 254, 0, &
39, 69, 156, 183, 198, 0, &
40, 70, 206, 208, 232, 0, &
3, 56, 106, 205, 0, 0, &
20, 119, 188, 230, 0, 0, &
51, 84, 112, 174, 226, 0, &
45, 111, 165, 251, 0, 0, &
60, 149, 169, 275, 0, 0, &
42, 77, 180, 202, 248, 0, &
62, 91, 121, 208, 241, 0, &
4, 95, 117, 236, 261, 0, &
49, 143, 206, 277, 0, 0, &
24, 99, 162, 182, 237, 0, &
29, 178, 179, 256, 0, 0, &
33, 106, 216, 243, 0, 0, &
12, 85, 104, 224, 270, 0, &
48, 87, 102, 192, 269, 0, &
56, 180, 185, 245, 0, 0, &
62, 184, 197, 255, 0, 0, &
47, 91, 178, 216, 254, 0, &
55, 122, 246, 268, 0, 0, &
54, 86, 124, 193, 273, 0, &
61, 70, 145, 150, 269, 0, &
26, 113, 193, 231, 0, 0, &
49, 89, 174, 194, 279, 0, &
1, 33, 77, 103, 213, 0, &
1, 204, 221, 270, 0, 0, &
27, 67, 138, 236, 277, 0, &
58, 83, 172, 239, 246, 0, &
11, 59, 199, 202, 246, 0, &
46, 153, 190, 250, 0, 0, &
41, 94, 113, 176, 265, 0, &
54, 132, 170, 266, 0, 0, &
3, 26, 72, 186, 203, 0, &
41, 68, 107, 255, 260, 0, &
63, 134, 222, 264, 0, 0, &
1, 43, 131, 249, 0, 0, &
32, 199, 219, 252, 0, 0, &
51, 53, 157, 235, 236, 0, &
2, 7, 114, 118, 0, 0, &
31, 217, 224, 260, 0, 0, &
37, 62, 137, 212, 264, 0, &
60, 99, 115, 205, 272, 0, &
20, 87, 90, 185, 194, 253, &
21, 92, 118, 148, 242, 0, &
14, 75, 121, 158, 209, 0, &
23, 210, 250, 261, 0, 0, &
2, 94, 133, 222, 274, 0, &
37, 101, 200, 249, 266, 0, &
64, 187, 258, 260, 0, 0, &
15, 81, 137, 183, 280, 0, &
59, 108, 140, 267, 0, 0, &
18, 142, 153, 280, 0, 0, &
16, 69, 159, 217, 276, 0, &
8, 82, 129, 228, 259, 0/
data nrw/ &
5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, &
5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, &
5,5,5,5,5,5,5,5,5,4,4,4,4,5,4,4,5,5,5,4, &
5,5,5,4,5,4,4,4,5,5,4,5,5,5,4,4,4,4,4,4, &
5,4,5,4,4,4,4,5,4,5,5,5,5,5,5,5,5,5,5,5, &
5,4,4,5,5,5,5,5,5,5,4,4,4,4,5,5,5,4,4,5, &
5,5,5,4,5,5,5,4,4,5,4,4,5,5,5,4,5,4,4,5, &
5,4,4,5,4,5,5,4,5,5,4,5,5,5,4,5,4,5,5,4, &
4,4,5,4,4,5,5,6,5,5,4,5,5,4,5,4,4,5,5/
ncw=3

View File

@ -0,0 +1,144 @@
program ldpcsim280_101
! End-to-end test of the (280,101)/crc24 encoder and decoders.
use packjt77
parameter(N=280, 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(280)
integer*1 cw(280)
integer*1 codeword(N),message101(101)
integer ncrc24
real rxdata(N),llr(N)
real llrd(280)
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. ldpcsim280_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 encode280_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 bpdecode280_101(llr,apmask,max_iterations,message101,cw,nharderror,niterations,nchecks)
dmin=0.0
if( (nharderror .lt. 0) .and. (ndeep .ge. 0) ) then
! call osd280_101(llr, Keff, apmask, ndeep, message101, cw, nharderror, dmin)
maxsuper=2
call decode280_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 ldpcsim280_101

403
lib/fsk4hf/osd280_101.f90 Normal file
View File

@ -0,0 +1,403 @@
subroutine osd280_101(llr,k,apmask,ndeep,message101,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (280,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=280
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 encode280_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=17
elseif(ndeep.eq.5) then
nord=3
npre1=1
npre2=1
nt=40
ntheta=12
ntau=15
elseif(ndeep.eq.6) then
nord=4
npre1=1
npre2=1
nt=95
ntheta=12
ntau=15
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 osd280_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