mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2026-06-06 16:05:09 -04:00
Rename fst240 to fst4 in Fortran routines.
This commit is contained in:
Steven Franke
@@ -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
|
||||
@@ -0,0 +1,154 @@
|
||||
subroutine decode240_101(llr,Keff,maxosd,norder,apmask,message101,cw,ntype,nharderror,dmin)
|
||||
!
|
||||
! A hybrid bp/osd decoder for the (240,101) code.
|
||||
!
|
||||
! maxosd<0: do bp only
|
||||
! maxosd=0: do bp and then call osd once with channel llrs
|
||||
! maxosd>1: do bp and then call osd maxosd times with saved bp outputs
|
||||
! norder : osd decoding depth
|
||||
!
|
||||
integer, parameter:: N=240, K=101, M=N-K
|
||||
integer*1 cw(N),apmask(N)
|
||||
integer*1 nxor(N),hdec(N)
|
||||
integer*1 message101(101),m101(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),zsave(N,3)
|
||||
real llr(N)
|
||||
real Tmn
|
||||
|
||||
include "ldpc_240_101_parity.f90"
|
||||
|
||||
maxiterations=30
|
||||
nosd=0
|
||||
if(maxosd.gt.3) maxosd=3
|
||||
if(maxosd.eq.0) then ! osd with channel llrs
|
||||
nosd=1
|
||||
zsave(:,1)=llr
|
||||
elseif(maxosd.gt.0) then !
|
||||
nosd=maxosd
|
||||
elseif(maxosd.lt.0) then ! just bp
|
||||
nosd=0
|
||||
endif
|
||||
|
||||
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
|
||||
zsum=0.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
|
||||
zsum=zsum+zn
|
||||
if(iter.gt.0 .and. iter.le.maxosd) then
|
||||
zsave(:,iter)=zsum
|
||||
endif
|
||||
|
||||
! 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
|
||||
enddo
|
||||
if( ncheck .eq. 0 ) then ! we have a codeword - if crc is good, return it
|
||||
m101=0
|
||||
m101(1:101)=cw(1:101)
|
||||
call get_crc24(m101,101,nbadcrc)
|
||||
if(nbadcrc.eq.0) then
|
||||
message101=cw(1:101)
|
||||
hdec=0
|
||||
where(llr .ge. 0) hdec=1
|
||||
nxor=ieor(hdec,cw)
|
||||
nharderror=sum(nxor)
|
||||
dmin=sum(nxor*abs(llr))
|
||||
ntype=1
|
||||
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
|
||||
exit
|
||||
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
|
||||
|
||||
do i=1,nosd
|
||||
zn=zsave(:,i)
|
||||
call osd240_101(zn,Keff,apmask,norder,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))
|
||||
ntype=2
|
||||
return
|
||||
endif
|
||||
enddo
|
||||
|
||||
ntype=0
|
||||
nharderror=-1
|
||||
dminosd=0.0
|
||||
|
||||
return
|
||||
end subroutine decode240_101
|
||||
@@ -0,0 +1,154 @@
|
||||
subroutine decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw,ntype,nharderror,dmin)
|
||||
!
|
||||
! A hybrid bp/osd decoder for the (240,74) code.
|
||||
!
|
||||
! maxosd<0: do bp only
|
||||
! maxosd=0: do bp and then call osd once with channel llrs
|
||||
! maxosd>1: do bp and then call osd maxosd times with saved bp outputs
|
||||
! norder : osd decoding depth
|
||||
!
|
||||
integer, parameter:: N=240, K=74, M=N-K
|
||||
integer*1 cw(N),apmask(N)
|
||||
integer*1 nxor(N),hdec(N)
|
||||
integer*1 message74(74),m74(74)
|
||||
integer nrw(M),ncw
|
||||
integer Nm(5,M)
|
||||
integer Mn(3,N) ! 3 checks per bit
|
||||
integer synd(M)
|
||||
real tov(3,N)
|
||||
real toc(5,M)
|
||||
real tanhtoc(5,M)
|
||||
real zn(N),zsum(N),zsave(N,3)
|
||||
real llr(N)
|
||||
real Tmn
|
||||
|
||||
include "ldpc_240_74_parity.f90"
|
||||
|
||||
maxiterations=30
|
||||
nosd=0
|
||||
if(maxosd.gt.3) maxosd=3
|
||||
if(maxosd.eq.0) then ! osd with channel llrs
|
||||
nosd=1
|
||||
zsave(:,1)=llr
|
||||
elseif(maxosd.gt.0) then !
|
||||
nosd=maxosd
|
||||
elseif(maxosd.lt.0) then ! just bp
|
||||
nosd=0
|
||||
endif
|
||||
|
||||
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
|
||||
zsum=0.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
|
||||
zsum=zsum+zn
|
||||
if(iter.gt.0 .and. iter.le.maxosd) then
|
||||
zsave(:,iter)=zsum
|
||||
endif
|
||||
|
||||
! 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
|
||||
enddo
|
||||
if( ncheck .eq. 0 ) then ! we have a codeword - if crc is good, return it
|
||||
m74=0
|
||||
m74(1:74)=cw(1:74)
|
||||
call get_crc24(m74,74,nbadcrc)
|
||||
if(nbadcrc.eq.0) then
|
||||
message74=cw(1:74)
|
||||
hdec=0
|
||||
where(llr .ge. 0) hdec=1
|
||||
nxor=ieor(hdec,cw)
|
||||
nharderror=sum(nxor)
|
||||
dmin=sum(nxor*abs(llr))
|
||||
ntype=1
|
||||
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
|
||||
exit
|
||||
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:5,i)=tanh(-toc(1:5,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
|
||||
|
||||
do i=1,nosd
|
||||
zn=zsave(:,i)
|
||||
call osd240_74(zn,Keff,apmask,norder,message74,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))
|
||||
ntype=2
|
||||
return
|
||||
endif
|
||||
enddo
|
||||
|
||||
ntype=0
|
||||
nharderror=-1
|
||||
dminosd=0.0
|
||||
|
||||
return
|
||||
end subroutine decode240_74
|
||||
@@ -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
|
||||
@@ -0,0 +1,46 @@
|
||||
subroutine encode240_74(message,codeword)
|
||||
use, intrinsic :: iso_c_binding
|
||||
use iso_c_binding, only: c_loc,c_size_t
|
||||
use crc
|
||||
|
||||
integer, parameter:: N=240, K=74, M=N-K
|
||||
character*24 c24
|
||||
integer*1 codeword(N)
|
||||
integer*1 gen(M,K)
|
||||
integer*1 message(K)
|
||||
integer*1 pchecks(M)
|
||||
integer*4 ncrc24
|
||||
include "ldpc_240_74_generator.f90"
|
||||
logical first
|
||||
data first/.true./
|
||||
save first,gen
|
||||
|
||||
if( first ) then ! fill the generator matrix
|
||||
gen=0
|
||||
do i=1,M
|
||||
do j=1,19
|
||||
read(g(i)(j:j),"(Z1)") istr
|
||||
ibmax=4
|
||||
if(j.eq.19) ibmax=2
|
||||
do jj=1, ibmax
|
||||
icol=(j-1)*4+jj
|
||||
if( btest(istr,4-jj) ) gen(i,icol)=1
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
first=.false.
|
||||
endif
|
||||
|
||||
do i=1,M
|
||||
nsum=0
|
||||
do j=1,K
|
||||
nsum=nsum+message(j)*gen(i,j)
|
||||
enddo
|
||||
pchecks(i)=mod(nsum,2)
|
||||
enddo
|
||||
|
||||
codeword(1:K)=message
|
||||
codeword(K+1:N)=pchecks
|
||||
|
||||
return
|
||||
end subroutine encode240_74
|
||||
@@ -0,0 +1,7 @@
|
||||
! FST4
|
||||
! LDPC(240,101)/CRC24 code, five 8x4 sync
|
||||
|
||||
parameter (KK=77) !Information bits (77 + CRC24)
|
||||
parameter (ND=120) !Data symbols
|
||||
parameter (NS=40) !Sync symbols
|
||||
parameter (NN=NS+ND) !Sync and data symbols (160)
|
||||
@@ -0,0 +1,155 @@
|
||||
program fst4sim
|
||||
|
||||
! Generate simulated signals for experimental slow FT4 mode
|
||||
|
||||
use wavhdr
|
||||
use packjt77
|
||||
include 'fst4_params.f90' !Set various constants
|
||||
type(hdr) h !Header for .wav file
|
||||
logical*1 wspr_hint
|
||||
character arg*12,fname*17
|
||||
character msg37*37,msgsent37*37,c77*77
|
||||
complex, allocatable :: c0(:)
|
||||
complex, allocatable :: c(:)
|
||||
real, allocatable :: wave(:)
|
||||
integer hmod
|
||||
integer itone(NN)
|
||||
integer*1 msgbits(101)
|
||||
integer*2, allocatable :: iwave(:) !Generated full-length waveform
|
||||
|
||||
! Get command-line argument(s)
|
||||
nargs=iargc()
|
||||
if(nargs.ne.10) then
|
||||
print*,'Need 10 arguments, got ',nargs
|
||||
print*,'Usage: fst4sim "message" TRsec f0 DT h fdop del nfiles snr W'
|
||||
print*,'Examples: fst4sim "K1JT K9AN EN50" 60 1500 0.0 1 0.1 1.0 10 -15 F'
|
||||
print*,'W (T or F) argument is hint to encoder to use WSPR message when there is abiguity'
|
||||
go to 999
|
||||
endif
|
||||
call getarg(1,msg37) !Message to be transmitted
|
||||
call getarg(2,arg)
|
||||
read(arg,*) nsec !TR sequence length, seconds
|
||||
call getarg(3,arg)
|
||||
read(arg,*) f00 !Frequency (only used for single-signal)
|
||||
call getarg(4,arg)
|
||||
read(arg,*) xdt !Time offset from nominal (s)
|
||||
call getarg(5,arg)
|
||||
read(arg,*) hmod !Modulation index, h
|
||||
call getarg(6,arg)
|
||||
read(arg,*) fspread !Watterson frequency spread (Hz)
|
||||
call getarg(7,arg)
|
||||
read(arg,*) delay !Watterson delay (ms)
|
||||
call getarg(8,arg)
|
||||
read(arg,*) nfiles !Number of files
|
||||
call getarg(9,arg)
|
||||
read(arg,*) snrdb !SNR_2500
|
||||
call getarg(10,arg)
|
||||
read(arg,*) wspr_hint !0:break ties as 77-bit 1:break ties as 50-bit
|
||||
|
||||
nfiles=abs(nfiles)
|
||||
twopi=8.0*atan(1.0)
|
||||
fs=12000.0 !Sample rate (Hz)
|
||||
dt=1.0/fs !Sample interval (s)
|
||||
nsps=0
|
||||
if(nsec.eq.15) nsps=720
|
||||
if(nsec.eq.30) nsps=1680
|
||||
if(nsec.eq.60) nsps=3888
|
||||
if(nsec.eq.120) nsps=8200
|
||||
if(nsec.eq.300) nsps=21504
|
||||
if(nsec.eq.900) nsps=66560
|
||||
if(nsec.eq.1800) nsps=134400
|
||||
if(nsps.eq.0) then
|
||||
print*,'Invalid TR sequence length.'
|
||||
go to 999
|
||||
endif
|
||||
baud=12000.0/nsps !Keying rate (baud)
|
||||
nmax=nsec*12000
|
||||
nz=nsps*NN
|
||||
txt=nz*dt !Transmission length (s)
|
||||
tt=nsps*dt !Duration of symbols (s)
|
||||
nwave=max(nmax,(NN+2)*nsps)
|
||||
allocate( c0(0:nwave-1) )
|
||||
allocate( c(0:nwave-1) )
|
||||
allocate( wave(nwave) )
|
||||
allocate( iwave(nmax) )
|
||||
|
||||
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
|
||||
|
||||
if(wspr_hint) then
|
||||
i3=0
|
||||
n3=6
|
||||
else
|
||||
i3=-1
|
||||
n3=-1
|
||||
endif
|
||||
call pack77(msg37,i3,n3,c77)
|
||||
if(i3.eq.0.and.n3.eq.6) iwspr=1
|
||||
call genfst4(msg37,0,msgsent37,msgbits,itone,iwspr)
|
||||
write(*,*)
|
||||
write(*,'(a9,a37,a3,L2,a7,i2)') 'Message: ',msgsent37,'W:',wspr_hint,' iwspr:',iwspr
|
||||
write(*,1000) f00,xdt,hmod,txt,snrdb
|
||||
1000 format('f0:',f9.3,' DT:',f6.2,' hmod:',i6,' 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(*,'(77i1,1x,24i1)') msgbits
|
||||
endif
|
||||
write(*,*)
|
||||
write(*,'(a17)') 'Channel symbols: '
|
||||
write(*,'(10i1)') itone
|
||||
write(*,*)
|
||||
|
||||
! call sgran()
|
||||
|
||||
fsample=12000.0
|
||||
icmplx=1
|
||||
f0=f00+1.5*hmod*baud
|
||||
call gen_fst4wave(itone,NN,nsps,nwave,fsample,hmod,f0,icmplx,c0,wave)
|
||||
k=nint((xdt+1.0)/dt)
|
||||
if(nsec.eq.15) k=nint((xdt+0.5)/dt)
|
||||
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,nwave,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(:size(iwave)))
|
||||
h=default_header(12000,nmax)
|
||||
if(nmax/12000.le.30) then
|
||||
write(fname,1102) ifile
|
||||
1102 format('000000_',i6.6,'.wav')
|
||||
else
|
||||
write(fname,1104) ifile
|
||||
1104 format('000000_',i4.4,'.wav')
|
||||
endif
|
||||
open(10,file=trim(fname),status='unknown',access='stream')
|
||||
write(10) h,iwave !Save to *.wav file
|
||||
close(10)
|
||||
write(*,1110) ifile,xdt,f00,snrdb,fname
|
||||
1110 format(i4,f7.2,f8.2,f7.1,2x,a17)
|
||||
enddo
|
||||
|
||||
999 end program fst4sim
|
||||
@@ -0,0 +1,91 @@
|
||||
subroutine gen_fst4wave(itone,nsym,nsps,nwave,fsample,hmod,f0, &
|
||||
icmplx,cwave,wave)
|
||||
|
||||
parameter(NTAB=65536)
|
||||
real wave(nwave)
|
||||
complex cwave(nwave),ctab(0:NTAB-1)
|
||||
real, allocatable, save :: pulse(:)
|
||||
real, allocatable :: dphi(:)
|
||||
integer hmod
|
||||
integer itone(nsym)
|
||||
logical first
|
||||
data first/.true./
|
||||
data nsps0/-99/
|
||||
save first,twopi,dt,tsym,nsps0,ctab
|
||||
|
||||
if(first) then
|
||||
twopi=8.0*atan(1.0)
|
||||
do i=0,NTAB-1
|
||||
phi=i*twopi/NTAB
|
||||
ctab(i)=cmplx(cos(phi),sin(phi))
|
||||
enddo
|
||||
endif
|
||||
|
||||
if(first.or.nsps.ne.nsps0) then
|
||||
if(allocated(pulse)) deallocate(pulse)
|
||||
allocate(pulse(1:3*nsps))
|
||||
dt=1.0/fsample
|
||||
tsym=nsps/fsample
|
||||
! Compute the smoothed frequency-deviation pulse
|
||||
do i=1,3*nsps
|
||||
tt=(i-1.5*nsps)/real(nsps)
|
||||
pulse(i)=gfsk_pulse(2.0,tt)
|
||||
enddo
|
||||
first=.false.
|
||||
nsps0=nsps
|
||||
endif
|
||||
|
||||
! Compute the smoothed frequency waveform.
|
||||
! Length = (nsym+2)*nsps samples, zero-padded
|
||||
allocate( dphi(0:(nsym+2)*nsps-1) )
|
||||
dphi_peak=twopi*hmod/real(nsps)
|
||||
dphi=0.0
|
||||
do j=1,nsym
|
||||
ib=(j-1)*nsps
|
||||
ie=ib+3*nsps-1
|
||||
dphi(ib:ie) = dphi(ib:ie) + dphi_peak*pulse(1:3*nsps)*itone(j)
|
||||
enddo
|
||||
|
||||
! Calculate and insert the audio waveform
|
||||
phi=0.0
|
||||
dphi = dphi + twopi*(f0-1.5*hmod/tsym)*dt !Shift frequency up by f0
|
||||
if(icmplx.eq.0) wave=0.
|
||||
if(icmplx.eq.1) cwave=0.
|
||||
k=0
|
||||
do j=0,(nsym+2)*nsps-1
|
||||
k=k+1
|
||||
i=phi*float(NTAB)/twopi
|
||||
i=iand(i,NTAB-1)
|
||||
if(icmplx.eq.0) then
|
||||
wave(k)=real(ctab(i))
|
||||
else
|
||||
cwave(k)=ctab(i)
|
||||
endif
|
||||
phi=phi+dphi(j)
|
||||
if(phi.gt.twopi) phi=phi-twopi
|
||||
enddo
|
||||
|
||||
! Compute the ramp-up and ramp-down symbols
|
||||
kshift=nsps
|
||||
if(icmplx.eq.0) then
|
||||
wave(1:nsps)=0.0
|
||||
wave(nsps+1:nsps+nsps/4)=wave(nsps+1:nsps+nsps/4) * &
|
||||
(1.0-cos(twopi*(/(i,i=0,nsps/4-1)/)/real(nsps/2)))/2.0
|
||||
k1=nsym*nsps+3*nsps/4+1
|
||||
wave((nsym+1)*nsps+1:)=0.0
|
||||
wave(k1:k1+nsps/4)=wave(k1:k1+nsps/4) * &
|
||||
(1.0+cos(twopi*(/(i,i=0,nsps/4)/)/real(nsps/2)))/2.0
|
||||
wave=cshift(wave,kshift)
|
||||
else
|
||||
cwave(1:nsps)=0.0
|
||||
cwave(nsps+1:nsps+nsps/4)=cwave(nsps+1:nsps+nsps/4) * &
|
||||
(1.0-cos(twopi*(/(i,i=0,nsps/4-1)/)/real(nsps/2)))/2.0
|
||||
k1=nsym*nsps+3*nsps/4+1
|
||||
cwave((nsym+1)*nsps+1:)=0.0
|
||||
cwave(k1:k1+nsps/4)=cwave(k1:k1+nsps/4) * &
|
||||
(1.0+cos(twopi*(/(i,i=0,nsps/4)/)/real(nsps/2)))/2.0
|
||||
cwave=cshift(cwave,kshift)
|
||||
endif
|
||||
|
||||
return
|
||||
end subroutine gen_fst4wave
|
||||
@@ -0,0 +1,108 @@
|
||||
subroutine genfst240_64(msg0,ichk,msgsent,msgbits,i4tone,iwspr)
|
||||
|
||||
! 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}
|
||||
! - iwspr 0: (240,101)/crc24, 1: (240,74)/crc24
|
||||
!
|
||||
! Frame structure:
|
||||
! s8 d30 s8 d30 s8 d30 s8 d30 s8
|
||||
|
||||
use packjt77
|
||||
include 'fst240_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 isyncword1(8),isyncword2(8)
|
||||
integer ncrc24
|
||||
integer graymap64(64)
|
||||
logical unpk77_success
|
||||
data isyncword1/3,1,4,0,6,5,2/
|
||||
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/
|
||||
data graymap64/ 0, 1, 3, 2, 6, 7, 5, 4,12,13,15,14,10,11, 9, 8, &
|
||||
24,25,27,26,30,31,29,28,20,21,23,22,18,19,17,16, &
|
||||
48,49,51,50,54,55,53,52,60,61,63,62,58,59,57,56, &
|
||||
40,41,43,42,46,47,45,44,36,37,39,38,34,35,33,32/
|
||||
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
|
||||
iwspr=0
|
||||
if(i3.eq.0.and.n3.eq.6) then
|
||||
iwspr=1
|
||||
read(c77,'(50i1)') msgbits(1:50)
|
||||
call get_crc24(msgbits,74,ncrc24)
|
||||
write(c24,'(b24.24)') ncrc24
|
||||
read(c24,'(24i1)') msgbits(51:74)
|
||||
else
|
||||
read(c77,'(77i1)') msgbits(1:77)
|
||||
msgbits(1:77)=mod(msgbits(1:77)+rvec,2)
|
||||
call get_crc24(msgbits,101,ncrc24)
|
||||
write(c24,'(b24.24)') ncrc24
|
||||
read(c24,'(24i1)') msgbits(78:101)
|
||||
endif
|
||||
|
||||
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_fst240_tones_from_bits(msgbits,i4tone,iwspr)
|
||||
|
||||
2 continue
|
||||
if(iwspr.eq.0) then
|
||||
call encode240_101(msgbits,codeword)
|
||||
else
|
||||
call encode240_74(msgbits(1:74),codeword)
|
||||
endif
|
||||
|
||||
! Grayscale mapping:
|
||||
! bits tone
|
||||
|
||||
do i=1,40
|
||||
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: 8)=isyncword1
|
||||
i4tone( 9: 38)=itmp( 1: 30)
|
||||
i4tone( 39: 46)=isyncword2
|
||||
i4tone( 47: 76)=itmp( 31: 60)
|
||||
i4tone( 77: 84)=isyncword1
|
||||
i4tone( 85:114)=itmp( 61: 90)
|
||||
i4tone(115:122)=isyncword2
|
||||
i4tone(123:152)=itmp( 91:120)
|
||||
i4tone(153:160)=isyncword1
|
||||
|
||||
999 return
|
||||
|
||||
end subroutine genfst240_64
|
||||
|
||||
subroutine graycode(in
|
||||
@@ -0,0 +1,111 @@
|
||||
subroutine genfst4(msg0,ichk,msgsent,msgbits,i4tone,iwspr)
|
||||
|
||||
! 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}
|
||||
! - iwspr in: 0: FST240 1: FST240W
|
||||
! out 0: (240,101)/crc24, 1: (240,74)/crc24
|
||||
!
|
||||
! Frame structure:
|
||||
! s8 d30 s8 d30 s8 d30 s8 d30 s8
|
||||
|
||||
use packjt77
|
||||
include 'fst4_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 isyncword1(8),isyncword2(8)
|
||||
integer ncrc24
|
||||
logical unpk77_success
|
||||
data isyncword1/0,1,3,2,1,0,2,3/
|
||||
data isyncword2/2,3,1,0,3,2,0,1/
|
||||
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
|
||||
if(iwspr.eq.1) then
|
||||
i3=0
|
||||
n3=6
|
||||
endif
|
||||
call pack77(message,i3,n3,c77)
|
||||
call unpack77(c77,0,msgsent,unpk77_success) !Unpack to get msgsent
|
||||
msgbits=0
|
||||
iwspr=0
|
||||
if(i3.eq.0.and.n3.eq.6) then
|
||||
iwspr=1
|
||||
read(c77,'(50i1)') msgbits(1:50)
|
||||
call get_crc24(msgbits,74,ncrc24)
|
||||
write(c24,'(b24.24)') ncrc24
|
||||
read(c24,'(24i1)') msgbits(51:74)
|
||||
else
|
||||
read(c77,'(77i1)') msgbits(1:77)
|
||||
msgbits(1:77)=mod(msgbits(1:77)+rvec,2)
|
||||
call get_crc24(msgbits,101,ncrc24)
|
||||
write(c24,'(b24.24)') ncrc24
|
||||
read(c24,'(24i1)') msgbits(78:101)
|
||||
endif
|
||||
|
||||
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_fst4_tones_from_bits(msgbits,i4tone,iwspr)
|
||||
|
||||
2 continue
|
||||
if(iwspr.eq.0) then
|
||||
call encode240_101(msgbits,codeword)
|
||||
else
|
||||
call encode240_74(msgbits(1:74),codeword)
|
||||
endif
|
||||
|
||||
! 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: 8)=isyncword1
|
||||
i4tone( 9: 38)=itmp( 1: 30)
|
||||
i4tone( 39: 46)=isyncword2
|
||||
i4tone( 47: 76)=itmp( 31: 60)
|
||||
i4tone( 77: 84)=isyncword1
|
||||
i4tone( 85:114)=itmp( 61: 90)
|
||||
i4tone(115:122)=isyncword2
|
||||
i4tone(123:152)=itmp( 91:120)
|
||||
i4tone(153:160)=isyncword1
|
||||
|
||||
999 return
|
||||
|
||||
end subroutine genfst4
|
||||
@@ -0,0 +1,25 @@
|
||||
subroutine get_crc24(mc,len,ncrc)
|
||||
!
|
||||
! 1. To calculate 24-bit CRC, mc(1:len-24) is the message and mc(len-23:len) are zero.
|
||||
! 2. To check a received CRC, mc(1:len) is the received message plus CRC.
|
||||
! ncrc will be zero if the received message/CRC are consistent.
|
||||
!
|
||||
character c24*24
|
||||
integer*1 mc(len)
|
||||
integer*1 r(25),p(25)
|
||||
integer ncrc
|
||||
! polynomial for 24-bit CRC 0x100065b
|
||||
data p/1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,1,1,0,1,1/
|
||||
|
||||
! divide by polynomial
|
||||
r=mc(1:25)
|
||||
do i=0,len-25
|
||||
r(25)=mc(i+25)
|
||||
r=mod(r+r(1)*p,2)
|
||||
r=cshift(r,1)
|
||||
enddo
|
||||
|
||||
write(c24,'(24b1)') r(1:24)
|
||||
read(c24,'(b24.24)') ncrc
|
||||
|
||||
end subroutine get_crc24
|
||||
@@ -0,0 +1,131 @@
|
||||
subroutine get_fst4_bitmetrics(cd,nss,hmod,nmax,nhicoh,bitmetrics,s4,badsync)
|
||||
|
||||
include 'fst4_params.f90'
|
||||
complex cd(0:NN*nss-1)
|
||||
complex cs(0:3,NN)
|
||||
complex csymb(nss)
|
||||
complex, allocatable, save :: c1(:,:) ! ideal waveforms, 20 samples per symbol, 4 tones
|
||||
complex cp(0:3) ! accumulated phase shift over symbol types 0:3
|
||||
complex csum,cterm
|
||||
integer isyncword1(0:7),isyncword2(0:7)
|
||||
integer graymap(0:3)
|
||||
integer ip(1)
|
||||
integer hmod
|
||||
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 isyncword1/0,1,3,2,1,0,2,3/
|
||||
data isyncword2/2,3,1,0,3,2,0,1/
|
||||
data graymap/0,1,3,2/
|
||||
data first/.true./,nss0/-1/
|
||||
save first,one,cp,nss0
|
||||
|
||||
if(nss.ne.nss0 .and. allocated(c1)) deallocate(c1)
|
||||
if(first .or. nss.ne.nss0) then
|
||||
allocate(c1(nss,0:3))
|
||||
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
|
||||
is4=0
|
||||
is5=0
|
||||
badsync=.false.
|
||||
ibmax=0
|
||||
|
||||
do k=1,8
|
||||
ip=maxloc(s4(:,k))
|
||||
if(isyncword1(k-1).eq.(ip(1)-1)) is1=is1+1
|
||||
ip=maxloc(s4(:,k+38))
|
||||
if(isyncword2(k-1).eq.(ip(1)-1)) is2=is2+1
|
||||
ip=maxloc(s4(:,k+76))
|
||||
if(isyncword1(k-1).eq.(ip(1)-1)) is3=is3+1
|
||||
ip=maxloc(s4(:,k+114))
|
||||
if(isyncword2(k-1).eq.(ip(1)-1)) is4=is4+1
|
||||
ip=maxloc(s4(:,k+152))
|
||||
if(isyncword1(k-1).eq.(ip(1)-1)) is5=is5+1
|
||||
enddo
|
||||
nsync=is1+is2+is3+is4+is5 !Number of correct hard sync symbols, 0-40
|
||||
badsync=.false.
|
||||
if(nsync .lt. 16) then
|
||||
badsync=.true.
|
||||
return
|
||||
endif
|
||||
|
||||
bitmetrics=0.0
|
||||
do nseq=1,nmax !Try coherent sequences of 1, 2, and 4 symbols
|
||||
if(nseq.eq.1) nsym=1
|
||||
if(nseq.eq.2) nsym=2
|
||||
if(nhicoh.eq.0) then
|
||||
if(nseq.eq.3) nsym=3
|
||||
if(nseq.eq.4) nsym=4
|
||||
else
|
||||
if(nseq.eq.3) nsym=4
|
||||
if(nseq.eq.4) nsym=8
|
||||
endif
|
||||
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.3) ibmax=5
|
||||
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_fst4_bitmetrics
|
||||
@@ -0,0 +1,131 @@
|
||||
subroutine get_fst4_bitmetrics2(cd,nss,hmod,nsizes,bitmetrics,s4hmod,badsync)
|
||||
|
||||
include 'fst4_params.f90'
|
||||
complex cd(0:NN*nss-1)
|
||||
complex csymb(nss)
|
||||
complex, allocatable, save :: c1(:,:) ! ideal waveforms, 4 tones
|
||||
complex cp(0:3) ! accumulated phase shift over symbol types 0:3
|
||||
complex csum,cterm
|
||||
integer isyncword1(0:7),isyncword2(0:7)
|
||||
integer graymap(0:3)
|
||||
integer ip(1)
|
||||
integer hmod
|
||||
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,4),s4hmod(0:3,NN)
|
||||
data isyncword1/0,1,3,2,1,0,2,3/
|
||||
data isyncword2/2,3,1,0,3,2,0,1/
|
||||
data graymap/0,1,3,2/
|
||||
data first/.true./,nss0/-1/
|
||||
save first,one,cp,nss0
|
||||
|
||||
if(nss.ne.nss0 .and. allocated(c1)) deallocate(c1)
|
||||
if(first .or. nss.ne.nss0) then
|
||||
allocate(c1(nss,0:3))
|
||||
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
|
||||
s4(itone,k,1)=abs(sum(csymb*conjg(c1(:,itone))))
|
||||
s4(itone,k,2)=abs(sum(csymb( 1:nss/2)*conjg(c1( 1:nss/2,itone)))) + &
|
||||
abs(sum(csymb(nss/2+1: nss)*conjg(c1(nss/2+1: nss,itone))))
|
||||
s4(itone,k,3)=abs(sum(csymb( 1: nss/4)*conjg(c1( 1: nss/4,itone)))) + &
|
||||
abs(sum(csymb( nss/4+1: nss/2)*conjg(c1( nss/4+1: nss/2,itone)))) + &
|
||||
abs(sum(csymb( nss/2+1:3*nss/4)*conjg(c1( nss/2+1:3*nss/4,itone)))) + &
|
||||
abs(sum(csymb(3*nss/4+1: nss)*conjg(c1(3*nss/4+1: nss,itone))))
|
||||
s4(itone,k,4)=abs(sum(csymb( 1: nss/8)*conjg(c1( 1: nss/8,itone)))) + &
|
||||
abs(sum(csymb( nss/8+1: nss/4)*conjg(c1( nss/8+1: nss/4,itone)))) + &
|
||||
abs(sum(csymb( nss/4+1:3*nss/8)*conjg(c1( nss/4+1:3*nss/8,itone)))) + &
|
||||
abs(sum(csymb(3*nss/8+1: nss/2)*conjg(c1(3*nss/8+1: nss/2,itone)))) + &
|
||||
abs(sum(csymb( nss/2+1:5*nss/8)*conjg(c1( nss/2+1:5*nss/8,itone)))) + &
|
||||
abs(sum(csymb(5*nss/8+1:3*nss/4)*conjg(c1(5*nss/8+1:3*nss/4,itone)))) + &
|
||||
abs(sum(csymb(3*nss/4+1:7*nss/8)*conjg(c1(3*nss/4+1:7*nss/8,itone)))) + &
|
||||
abs(sum(csymb(7*nss/8+1: nss)*conjg(c1(7*nss/8+1: nss,itone))))
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Sync quality check
|
||||
is1=0
|
||||
is2=0
|
||||
is3=0
|
||||
is4=0
|
||||
is5=0
|
||||
badsync=.false.
|
||||
ibmax=0
|
||||
|
||||
is1=0; is2=0; is3=0; is4=0; is5=0
|
||||
do k=1,8
|
||||
ip=maxloc(s4(:,k,1))
|
||||
if(isyncword1(k-1).eq.(ip(1)-1)) is1=is1+1
|
||||
ip=maxloc(s4(:,k+38,1))
|
||||
if(isyncword2(k-1).eq.(ip(1)-1)) is2=is2+1
|
||||
ip=maxloc(s4(:,k+76,1))
|
||||
if(isyncword1(k-1).eq.(ip(1)-1)) is3=is3+1
|
||||
ip=maxloc(s4(:,k+114,1))
|
||||
if(isyncword2(k-1).eq.(ip(1)-1)) is4=is4+1
|
||||
ip=maxloc(s4(:,k+152,1))
|
||||
if(isyncword1(k-1).eq.(ip(1)-1)) is5=is5+1
|
||||
enddo
|
||||
nsync=is1+is2+is3+is4+is5 !Number of correct hard sync symbols, 0-40
|
||||
badsync=.false.
|
||||
|
||||
if(nsync .lt. 16) then
|
||||
badsync=.true.
|
||||
return
|
||||
endif
|
||||
|
||||
bitmetrics=0.0
|
||||
do nsub=1,nsizes
|
||||
do ks=1,NN
|
||||
s2=0
|
||||
do i=0,3
|
||||
s2(i)=s4(graymap(i),ks,nsub)
|
||||
enddo
|
||||
ipt=1+(ks-1)*2
|
||||
ibmax=1
|
||||
do ib=0,ibmax
|
||||
bm=maxval(s2(0:3),one(0:3,ibmax-ib)) - &
|
||||
maxval(s2(0:3),.not.one(0:3,ibmax-ib))
|
||||
if(ipt+ib.gt.2*NN) cycle
|
||||
bitmetrics(ipt+ib,nsub)=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 the s4 array corresponding to N=1/hmod. Will be used for SNR calculation
|
||||
if(hmod.eq.1) s4hmod(:,:)=s4(:,:,1)
|
||||
if(hmod.eq.2) s4hmod(:,:)=s4(:,:,2)
|
||||
if(hmod.eq.4) s4hmod(:,:)=s4(:,:,3)
|
||||
if(hmod.eq.8) s4hmod(:,:)=s4(:,:,4)
|
||||
return
|
||||
|
||||
end subroutine get_fst4_bitmetrics2
|
||||
@@ -0,0 +1,142 @@
|
||||
character*26 g(139)
|
||||
|
||||
data g/ &
|
||||
"e28df133efbc554bcd30eb1828", &
|
||||
"b1adf97787f81b4ac02e0caff8", &
|
||||
"e70c43adce5036f847af367560", &
|
||||
"c26663f7f7acafdf5abacb6f30", &
|
||||
"eba93204ddfa3bcf994aea8998", &
|
||||
"126b51e33c6a740afa0d5ce990", &
|
||||
"b41a1569e6fede1f2f5395cb68", &
|
||||
"1d3af0bb43fddbc670a291cc70", &
|
||||
"e0aebd9921e2c9e1d453ffccb0", &
|
||||
"897d1370f0df94b8b27a5e4fb8", &
|
||||
"5e97539338003b13fa8198ad38", &
|
||||
"7276b87da4a4d777e2752fdd48", &
|
||||
"989888bd3a85835e2bc6a560f8", &
|
||||
"7ec4f4a56199ab0a8d6e102478", &
|
||||
"207007665090258782d1b38a98", &
|
||||
"1ea1f61cd7f0b7eed7dd346ab8", &
|
||||
"08f150b27c7f18a027783de0e8", &
|
||||
"d42324a4e21b62d548d7865858", &
|
||||
"2e029656269d4fe46e167d21d0", &
|
||||
"7d84acb7737b0ca6b6f2ef5eb0", &
|
||||
"6674ca04528ad4782bf5e15248", &
|
||||
"118ce9825f563ae4963af7a0b0", &
|
||||
"fb06248cc985e314b1b36ccd38", &
|
||||
"1c478b7a5aec7e1cfc9c24eb70", &
|
||||
"185a0f06a84f7f4f484c455020", &
|
||||
"98b840a3a70688cd58588e3e30", &
|
||||
"cfb7719de83a3baf582e5b2aa0", &
|
||||
"9d8cc6b5a01fdbfa307a769048", &
|
||||
"ed776a728ca162d6fcc8996760", &
|
||||
"8d2b068128dfb2f8d22c79db50", &
|
||||
"bd2ba50007789ffb7324aa9190", &
|
||||
"fd95008fe88812025e78065610", &
|
||||
"3027849be8e99f9ef68eac1020", &
|
||||
"88574e1ea39d87414b15e803a8", &
|
||||
"89365b330e76e6dde740dced08", &
|
||||
"c83f37b913ed0f6b802aaf21d8", &
|
||||
"bdca7c1959caa7488b7eb13030", &
|
||||
"794e0b4888e1ef42992287dd98", &
|
||||
"526ac87fbaa790c6cd58864e08", &
|
||||
"940518ba1a51c1da55bc8b2d70", &
|
||||
"59c5e51ebfbd02ab30ff822378", &
|
||||
"c81fff87866e04f8f3948c7f10", &
|
||||
"7913513f3e2a3c0f76b69f6d68", &
|
||||
"e43cc04da189c44803c4f740a0", &
|
||||
"fdca7c1959ca85488b7eb13030", &
|
||||
"95b07fce9b7b1bf4f057ca61b8", &
|
||||
"d7db48a86691a0c0c9305aac90", &
|
||||
"0d50bf79a59464597c43ba8058", &
|
||||
"4a9c34b23fd5eaff8c9dc215e0", &
|
||||
"3d5305a6f0427938eeb9d1c118", &
|
||||
"55d8b6b58039f7a3a2d592a900", &
|
||||
"784f349ecb74c4abbdbb073b90", &
|
||||
"5973bbb2205f9d6a5c9a55c238", &
|
||||
"5d2ee61006fec94f69f6b0f460", &
|
||||
"9e1f52ef1e6589990dd0ce0cc8", &
|
||||
"85b7b48f4b45775c9f8a36cc90", &
|
||||
"ae1d6a0171168f6d70804b79f8", &
|
||||
"a467aa9aa6cdc7094677c730d8", &
|
||||
"dcf2f56c9ae20fb57e89b916d0", &
|
||||
"3ae98d26ae96ea714c1a5146d0", &
|
||||
"103c89581446805b8c71b2e638", &
|
||||
"6783f3dfec835dd4e92131cc20", &
|
||||
"52f88428c50f12c55876f7d8a8", &
|
||||
"51fcb0e56a22fa3b7140aeaa80", &
|
||||
"07c54871155603e65325f66cd8", &
|
||||
"a8dd4fac47a113ee5706eef180", &
|
||||
"f6cdc6f4cc1fa7e4db15bf86f8", &
|
||||
"2e1c6a0171168f6d70c04a79f8", &
|
||||
"2a90ab82bef6424db981752dc8", &
|
||||
"845a1db59c193249d937e889d0", &
|
||||
"a929d379f1769cb4baa4e41e90", &
|
||||
"0c2a5829548d82223d6f566d48", &
|
||||
"420087bc5c4e2f5bc139ad0220", &
|
||||
"6df8d880ae7209fe52c69ede00", &
|
||||
"dfbdcef29a985fd40d052d1a88", &
|
||||
"8567fc332342b1ed8408f5fa00", &
|
||||
"c908feb4e1866a24ca0c702a08", &
|
||||
"645f5ee59f9f64fd43a5f2ec30", &
|
||||
"bee56991e877baf3e9cf11b770", &
|
||||
"649ea2e4194ca51be28abf3430", &
|
||||
"90e7394c551bd58d00686d5420", &
|
||||
"4e3cf731f8f89e8414214afaf0", &
|
||||
"dcbf16aa8180a7712571e94f98", &
|
||||
"9b456c015999c52b7fbd1ab390", &
|
||||
"397ab76924659c4b8b3be4ac58", &
|
||||
"4f5038c4f9da4b02bdfa178278", &
|
||||
"4892fada978c98dd4fd363c450", &
|
||||
"6c8af64b426bc474431c110c98", &
|
||||
"84a553be5ef0e57390a5af05b0", &
|
||||
"bed4a9347c9a2064f6d63ac0f8", &
|
||||
"d973bbb2605f9d6a5c9a57c238", &
|
||||
"1e3bee9a99fe10d3864ee669d8", &
|
||||
"a590771ff185d807cb32f46000", &
|
||||
"9a498fc4b549d81c625f80fc90", &
|
||||
"28b3e72878aadee7e0e2617950", &
|
||||
"96ce025d621a91396aa8f3ec20", &
|
||||
"4f5a77becf838a590d6d406ea8", &
|
||||
"52d3856dfb9fe78012f10e25c0", &
|
||||
"b45323c2b28b4752ca0675d2e0", &
|
||||
"3bae5a8452a785beb35851ad18", &
|
||||
"65098832d20d915e75bea336e8", &
|
||||
"5eb6f3c331098e8c0fbfa3aee0", &
|
||||
"ef19d974a25540c8998fbf1df0", &
|
||||
"403ea58feff08cf92d5cacc780", &
|
||||
"6ba93204ddfa7bcb994aea8998", &
|
||||
"653909166aa7bead4bd9c90020", &
|
||||
"089cb20e639bc5a44da66f17c0", &
|
||||
"10f803949961359e994f5ade88", &
|
||||
"15b7ec1e6106cd55ef7d996590", &
|
||||
"c99e99de9d85d2b999a17a95d8", &
|
||||
"ca3e161b97148bac6dd28a6178", &
|
||||
"e1ab199c992cb4c22aee115358", &
|
||||
"ea8a4d0e96d3d9f827899b6d88", &
|
||||
"8af4992d60223f021569a8ab60", &
|
||||
"5087771abceb87a6d872291fe8", &
|
||||
"d045e0812e217bb7bbdac92f30", &
|
||||
"ccccd78ae5fa6e191f21c06908", &
|
||||
"54545f37df6fed4734ef6509b0", &
|
||||
"b0780327d899cbc03d95a81a48", &
|
||||
"a4229c31f2b85e44a322273d50", &
|
||||
"d182ab001c2085ea7be26a20d0", &
|
||||
"1a82c30b4fba7dfaafb8d287a8", &
|
||||
"d974fba598e7fb0630c1587db0", &
|
||||
"b5c078a8cbab3e73728659ea20", &
|
||||
"626bbf9eed1a8715c3a7d38f60", &
|
||||
"c1efe9aa67130865fda93d8be8", &
|
||||
"d39796dbce155df6306e7b77c0", &
|
||||
"c7e7c1f032d7209b4549e84aa8", &
|
||||
"d5799b30a1605baf6b9cd04960", &
|
||||
"0baf2d21051a926dfd87046d70", &
|
||||
"da8bf7d1e305c499b573c02cc8", &
|
||||
"0ccaa7fffb9ae3e42dd0688328", &
|
||||
"b951b62e18f5290ac13c195130", &
|
||||
"79b006f001961fb233be80d0e8", &
|
||||
"56637b6dedfd6e050f06404a48", &
|
||||
"e0c4bf71a15597523bbd57bde0", &
|
||||
"1312231ffa04426a34a8fab038", &
|
||||
"db5f6f0455d24b8358d1cbc3d8", &
|
||||
"d559e31b34d21f48e1f501af30"/
|
||||
@@ -0,0 +1,393 @@
|
||||
data Mn/ &
|
||||
57, 100, 134, &
|
||||
56, 99, 136, &
|
||||
1, 12, 15, &
|
||||
2, 23, 72, &
|
||||
3, 133, 137, &
|
||||
4, 93, 125, &
|
||||
5, 68, 139, &
|
||||
6, 38, 55, &
|
||||
7, 40, 78, &
|
||||
8, 30, 84, &
|
||||
9, 17, 122, &
|
||||
10, 34, 95, &
|
||||
11, 36, 138, &
|
||||
13, 90, 132, &
|
||||
14, 50, 117, &
|
||||
16, 57, 83, &
|
||||
18, 22, 121, &
|
||||
19, 60, 89, &
|
||||
20, 98, 107, &
|
||||
21, 37, 61, &
|
||||
24, 26, 75, &
|
||||
25, 88, 115, &
|
||||
27, 49, 127, &
|
||||
28, 74, 119, &
|
||||
29, 111, 114, &
|
||||
31, 91, 129, &
|
||||
32, 96, 104, &
|
||||
30, 33, 130, &
|
||||
35, 65, 135, &
|
||||
41, 42, 87, &
|
||||
44, 108, 131, &
|
||||
45, 94, 101, &
|
||||
45, 46, 97, &
|
||||
47, 102, 134, &
|
||||
48, 64, 104, &
|
||||
19, 51, 116, &
|
||||
20, 52, 67, &
|
||||
53, 104, 113, &
|
||||
12, 54, 103, &
|
||||
58, 66, 88, &
|
||||
62, 80, 124, &
|
||||
63, 70, 71, &
|
||||
73, 114, 123, &
|
||||
76, 85, 128, &
|
||||
77, 106, 109, &
|
||||
46, 79, 126, &
|
||||
61, 81, 110, &
|
||||
82, 92, 120, &
|
||||
86, 105, 112, &
|
||||
66, 100, 118, &
|
||||
23, 51, 136, &
|
||||
1, 40, 53, &
|
||||
2, 73, 81, &
|
||||
3, 63, 130, &
|
||||
4, 68, 136, &
|
||||
5, 60, 78, &
|
||||
6, 72, 131, &
|
||||
7, 115, 124, &
|
||||
8, 89, 120, &
|
||||
9, 15, 44, &
|
||||
10, 22, 93, &
|
||||
11, 49, 100, &
|
||||
13, 55, 80, &
|
||||
14, 76, 95, &
|
||||
16, 54, 111, &
|
||||
17, 41, 110, &
|
||||
18, 69, 139, &
|
||||
21, 24, 116, &
|
||||
25, 39, 71, &
|
||||
26, 69, 90, &
|
||||
27, 101, 133, &
|
||||
28, 64, 126, &
|
||||
29, 94, 103, &
|
||||
31, 56, 57, &
|
||||
32, 91, 102, &
|
||||
33, 35, 129, &
|
||||
34, 47, 128, &
|
||||
36, 86, 117, &
|
||||
37, 74, 75, &
|
||||
38, 79, 106, &
|
||||
42, 82, 123, &
|
||||
43, 77, 99, &
|
||||
48, 70, 92, &
|
||||
50, 109, 118, &
|
||||
52, 112, 119, &
|
||||
58, 62, 108, &
|
||||
59, 84, 134, &
|
||||
57, 65, 122, &
|
||||
67, 97, 113, &
|
||||
83, 127, 135, &
|
||||
85, 121, 125, &
|
||||
87, 132, 137, &
|
||||
96, 98, 105, &
|
||||
73, 107, 138, &
|
||||
1, 83, 89, &
|
||||
2, 41, 70, &
|
||||
3, 35, 131, &
|
||||
4, 111, 128, &
|
||||
5, 29, 99, &
|
||||
6, 25, 31, &
|
||||
7, 19, 96, &
|
||||
1, 39, 110, &
|
||||
2, 7, 117, &
|
||||
3, 49, 109, &
|
||||
4, 81, 96, &
|
||||
5, 100, 108, &
|
||||
6, 51, 124, &
|
||||
2, 20, 132, &
|
||||
8, 80, 137, &
|
||||
9, 56, 67, &
|
||||
10, 63, 102, &
|
||||
11, 16, 101, &
|
||||
12, 115, 122, &
|
||||
13, 32, 128, &
|
||||
14, 15, 130, &
|
||||
14, 70, 99, &
|
||||
11, 51, 69, &
|
||||
17, 89, 105, &
|
||||
18, 83, 99, &
|
||||
19, 44, 79, &
|
||||
20, 106, 133, &
|
||||
10, 21, 123, &
|
||||
22, 23, 61, &
|
||||
16, 22, 60, &
|
||||
24, 38, 114, &
|
||||
25, 37, 42, &
|
||||
26, 43, 52, &
|
||||
27, 68, 71, &
|
||||
28, 65, 139, &
|
||||
29, 62, 69, &
|
||||
30, 92, 126, &
|
||||
31, 78, 123, &
|
||||
13, 44, 78, &
|
||||
33, 40, 120, &
|
||||
7, 34, 119, &
|
||||
4, 35, 77, &
|
||||
12, 36, 52, &
|
||||
25, 98, 136, &
|
||||
5, 24, 133, &
|
||||
1, 80, 91, &
|
||||
33, 96, 97, &
|
||||
34, 41, 91, &
|
||||
32, 37, 117, &
|
||||
26, 72, 125, &
|
||||
19, 65, 75, &
|
||||
45, 131, 136, &
|
||||
46, 55, 70, &
|
||||
47, 48, 50, &
|
||||
6, 48, 94, &
|
||||
3, 74, 79, &
|
||||
39, 50, 126, &
|
||||
23, 118, 127, &
|
||||
21, 36, 113, &
|
||||
53, 77, 134, &
|
||||
30, 54, 55, &
|
||||
17, 46, 135, &
|
||||
9, 92, 102, &
|
||||
57, 85, 87, &
|
||||
58, 125, 138, &
|
||||
59, 76, 93, &
|
||||
60, 66, 107, &
|
||||
47, 132, 138, &
|
||||
29, 85, 131, &
|
||||
43, 73, 108, &
|
||||
64, 75, 129, &
|
||||
28, 38, 53, &
|
||||
61, 106, 122, &
|
||||
56, 71, 114, &
|
||||
27, 57, 120, &
|
||||
62, 67, 130, &
|
||||
54, 104, 118, &
|
||||
8, 68, 115, &
|
||||
72, 86, 111, &
|
||||
73, 74, 94, &
|
||||
49, 105, 113, &
|
||||
42, 86, 121, &
|
||||
40, 59, 109, &
|
||||
35, 88, 95, &
|
||||
31, 107, 112, &
|
||||
58, 64, 87, &
|
||||
68, 79, 104, &
|
||||
1, 5, 121, &
|
||||
15, 82, 93, &
|
||||
18, 88, 116, &
|
||||
82, 84, 119, &
|
||||
7, 71, 103, &
|
||||
4, 80, 94, &
|
||||
63, 81, 84, &
|
||||
66, 76, 137, &
|
||||
83, 124, 129, &
|
||||
90, 112, 116, &
|
||||
89, 111, 134, &
|
||||
6, 21, 120, &
|
||||
3, 16, 25, &
|
||||
12, 28, 131, &
|
||||
45, 95, 110, &
|
||||
17, 93, 124, &
|
||||
97, 121, 127, &
|
||||
98, 103, 135, &
|
||||
8, 99, 138, &
|
||||
41, 101, 139, &
|
||||
13, 24, 105, &
|
||||
14, 53, 107, &
|
||||
10, 64, 98, &
|
||||
11, 35, 78, &
|
||||
90, 100, 103, &
|
||||
9, 72, 101, &
|
||||
18, 74, 92, &
|
||||
15, 73, 87, &
|
||||
2, 88, 113, &
|
||||
20, 55, 85, &
|
||||
19, 67, 110, &
|
||||
26, 27, 95, &
|
||||
22, 50, 114, &
|
||||
29, 49, 81, &
|
||||
32, 52, 83, &
|
||||
30, 37, 77, &
|
||||
39, 128, 135, &
|
||||
23, 128, 130, &
|
||||
36, 76, 126, &
|
||||
33, 132, 139, &
|
||||
34, 89, 118, &
|
||||
38, 58, 127, &
|
||||
31, 54, 125, &
|
||||
40, 70, 75, &
|
||||
41, 109, 116, &
|
||||
43, 60, 63, &
|
||||
44, 84, 86, &
|
||||
42, 47, 62, &
|
||||
45, 82, 90, &
|
||||
43, 46, 91, &
|
||||
48, 112, 122, &
|
||||
51, 102, 133, &
|
||||
59, 61, 108, &
|
||||
65, 117, 137, &
|
||||
56, 66, 96, &
|
||||
59, 69, 104, &
|
||||
39, 69, 119, &
|
||||
97, 115, 123, &
|
||||
106, 111, 129/
|
||||
|
||||
data Nm/ &
|
||||
3, 52, 95, 102, 140, 182, &
|
||||
4, 53, 96, 103, 108, 210, &
|
||||
5, 54, 97, 104, 150, 194, &
|
||||
6, 55, 98, 105, 136, 187, &
|
||||
7, 56, 99, 106, 139, 182, &
|
||||
8, 57, 100, 107, 149, 193, &
|
||||
9, 58, 101, 103, 135, 186, &
|
||||
10, 59, 109, 172, 200, 0, &
|
||||
11, 60, 110, 157, 207, 0, &
|
||||
12, 61, 111, 122, 204, 0, &
|
||||
13, 62, 112, 117, 205, 0, &
|
||||
3, 39, 113, 137, 195, 0, &
|
||||
14, 63, 114, 133, 202, 0, &
|
||||
15, 64, 115, 116, 203, 0, &
|
||||
3, 60, 115, 183, 209, 0, &
|
||||
16, 65, 112, 124, 194, 0, &
|
||||
11, 66, 118, 156, 197, 0, &
|
||||
17, 67, 119, 184, 208, 0, &
|
||||
18, 36, 101, 120, 145, 212, &
|
||||
19, 37, 108, 121, 211, 0, &
|
||||
20, 68, 122, 153, 193, 0, &
|
||||
17, 61, 123, 124, 214, 0, &
|
||||
4, 51, 123, 152, 219, 0, &
|
||||
21, 68, 125, 139, 202, 0, &
|
||||
22, 69, 100, 126, 138, 194, &
|
||||
21, 70, 127, 144, 213, 0, &
|
||||
23, 71, 128, 169, 213, 0, &
|
||||
24, 72, 129, 166, 195, 0, &
|
||||
25, 73, 99, 130, 163, 215, &
|
||||
10, 28, 131, 155, 217, 0, &
|
||||
26, 74, 100, 132, 179, 224, &
|
||||
27, 75, 114, 143, 216, 0, &
|
||||
28, 76, 134, 141, 221, 0, &
|
||||
12, 77, 135, 142, 222, 0, &
|
||||
29, 76, 97, 136, 178, 205, &
|
||||
13, 78, 137, 153, 220, 0, &
|
||||
20, 79, 126, 143, 217, 0, &
|
||||
8, 80, 125, 166, 223, 0, &
|
||||
69, 102, 151, 218, 238, 0, &
|
||||
9, 52, 134, 177, 225, 0, &
|
||||
30, 66, 96, 142, 201, 226, &
|
||||
30, 81, 126, 176, 229, 0, &
|
||||
82, 127, 164, 227, 231, 0, &
|
||||
31, 60, 120, 133, 228, 0, &
|
||||
32, 33, 146, 196, 230, 0, &
|
||||
33, 46, 147, 156, 231, 0, &
|
||||
34, 77, 148, 162, 229, 0, &
|
||||
35, 83, 148, 149, 232, 0, &
|
||||
23, 62, 104, 175, 215, 0, &
|
||||
15, 84, 148, 151, 214, 0, &
|
||||
36, 51, 107, 117, 233, 0, &
|
||||
37, 85, 127, 137, 216, 0, &
|
||||
38, 52, 154, 166, 203, 0, &
|
||||
39, 65, 155, 171, 224, 0, &
|
||||
8, 63, 147, 155, 211, 0, &
|
||||
2, 74, 110, 168, 236, 0, &
|
||||
1, 16, 74, 88, 158, 169, &
|
||||
40, 86, 159, 180, 223, 0, &
|
||||
87, 160, 177, 234, 237, 0, &
|
||||
18, 56, 124, 161, 227, 0, &
|
||||
20, 47, 123, 167, 234, 0, &
|
||||
41, 86, 130, 170, 229, 0, &
|
||||
42, 54, 111, 188, 227, 0, &
|
||||
35, 72, 165, 180, 204, 0, &
|
||||
29, 88, 129, 145, 235, 0, &
|
||||
40, 50, 161, 189, 236, 0, &
|
||||
37, 89, 110, 170, 212, 0, &
|
||||
7, 55, 128, 172, 181, 0, &
|
||||
67, 70, 117, 130, 237, 238, &
|
||||
42, 83, 96, 116, 147, 225, &
|
||||
42, 69, 128, 168, 186, 0, &
|
||||
4, 57, 144, 173, 207, 0, &
|
||||
43, 53, 94, 164, 174, 209, &
|
||||
24, 79, 150, 174, 208, 0, &
|
||||
21, 79, 145, 165, 225, 0, &
|
||||
44, 64, 160, 189, 220, 0, &
|
||||
45, 82, 136, 154, 217, 0, &
|
||||
9, 56, 132, 133, 205, 0, &
|
||||
46, 80, 120, 150, 181, 0, &
|
||||
41, 63, 109, 140, 187, 0, &
|
||||
47, 53, 105, 188, 215, 0, &
|
||||
48, 81, 183, 185, 230, 0, &
|
||||
16, 90, 95, 119, 190, 216, &
|
||||
10, 87, 185, 188, 228, 0, &
|
||||
44, 91, 158, 163, 211, 0, &
|
||||
49, 78, 173, 176, 228, 0, &
|
||||
30, 92, 158, 180, 209, 0, &
|
||||
22, 40, 178, 184, 210, 0, &
|
||||
18, 59, 95, 118, 192, 222, &
|
||||
14, 70, 191, 206, 230, 0, &
|
||||
26, 75, 140, 142, 231, 0, &
|
||||
48, 83, 131, 157, 208, 0, &
|
||||
6, 61, 160, 183, 197, 0, &
|
||||
32, 73, 149, 174, 187, 0, &
|
||||
12, 64, 178, 196, 213, 0, &
|
||||
27, 93, 101, 105, 141, 236, &
|
||||
33, 89, 141, 198, 239, 0, &
|
||||
19, 93, 138, 199, 204, 0, &
|
||||
2, 82, 99, 116, 119, 200, &
|
||||
1, 50, 62, 106, 206, 0, &
|
||||
32, 71, 112, 201, 207, 0, &
|
||||
34, 75, 111, 157, 233, 0, &
|
||||
39, 73, 186, 199, 206, 0, &
|
||||
27, 35, 38, 171, 181, 237, &
|
||||
49, 93, 118, 175, 202, 0, &
|
||||
45, 80, 121, 167, 240, 0, &
|
||||
19, 94, 161, 179, 203, 0, &
|
||||
31, 86, 106, 164, 234, 0, &
|
||||
45, 84, 104, 177, 226, 0, &
|
||||
47, 66, 102, 196, 212, 0, &
|
||||
25, 65, 98, 173, 192, 240, &
|
||||
49, 85, 179, 191, 232, 0, &
|
||||
38, 89, 153, 175, 210, 0, &
|
||||
25, 43, 125, 168, 214, 0, &
|
||||
22, 58, 113, 172, 239, 0, &
|
||||
36, 68, 184, 191, 226, 0, &
|
||||
15, 78, 103, 143, 235, 0, &
|
||||
50, 84, 152, 171, 222, 0, &
|
||||
24, 85, 135, 185, 238, 0, &
|
||||
48, 59, 134, 169, 193, 0, &
|
||||
17, 91, 176, 182, 198, 0, &
|
||||
11, 88, 113, 167, 232, 0, &
|
||||
43, 81, 122, 132, 239, 0, &
|
||||
41, 58, 107, 190, 197, 0, &
|
||||
6, 91, 144, 159, 224, 0, &
|
||||
46, 72, 131, 151, 220, 0, &
|
||||
23, 90, 152, 198, 223, 0, &
|
||||
44, 77, 98, 114, 218, 219, &
|
||||
26, 76, 165, 190, 240, 0, &
|
||||
28, 54, 115, 170, 219, 0, &
|
||||
31, 57, 97, 146, 163, 195, &
|
||||
14, 92, 108, 162, 221, 0, &
|
||||
5, 71, 121, 139, 233, 0, &
|
||||
1, 34, 87, 154, 192, 0, &
|
||||
29, 90, 156, 199, 218, 0, &
|
||||
2, 51, 55, 138, 146, 0, &
|
||||
5, 92, 109, 189, 235, 0, &
|
||||
13, 94, 159, 162, 200, 0, &
|
||||
7, 67, 129, 201, 221, 0/
|
||||
|
||||
data nrw/ &
|
||||
6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,6,5, &
|
||||
5,5,5,5,6,5,5,5,6,5,6,5,5,5,6,5,5,5,5,5, &
|
||||
6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,5,5,5, &
|
||||
5,5,5,5,5,5,5,5,6,6,5,5,6,5,5,5,5,5,5,5, &
|
||||
5,5,6,5,5,5,5,5,6,5,5,5,5,5,5,6,5,5,6,5, &
|
||||
5,5,5,6,5,5,5,5,5,5,6,5,5,5,5,5,5,5,5,5, &
|
||||
5,5,5,5,5,5,5,6,5,5,6,5,5,5,5,5,5,5,5/
|
||||
|
||||
ncw=3
|
||||
@@ -0,0 +1,170 @@
|
||||
character*19 g(166)
|
||||
|
||||
data g/ &
|
||||
"de8b3201e3c59f55a14", &
|
||||
"2e06d352ebc5b74c4fc", &
|
||||
"2e16d6cf5a725c3244c", &
|
||||
"84f5587edca6d777de4", &
|
||||
"e152b1e2b5965093ecc", &
|
||||
"244b4828a2ccf2b5f58", &
|
||||
"5fbbaade810e123c730", &
|
||||
"6b7e92a99a918df3d44", &
|
||||
"bbcec6a63ab757a7278", &
|
||||
"f5f3f0b89a21ceccdb0", &
|
||||
"a248c5f1ec2bc816290", &
|
||||
"c84bbad839a5fe76d0c", &
|
||||
"ad724129bbf4c7f4570", &
|
||||
"91adb56e7623a2575cc", &
|
||||
"cbe995bdf156df2c9e4", &
|
||||
"92ff6ea492c08c150e0", &
|
||||
"c4ddbe5a02f6a933384", &
|
||||
"d2e9befc131dc483858", &
|
||||
"68567543d1eebcb080c", &
|
||||
"21fa61d559f9baf6abc", &
|
||||
"911c4fbbafc72e3db28", &
|
||||
"7c0b534af4b7d583d50", &
|
||||
"12ce371b90ee9dfe72c", &
|
||||
"15a604148872e251ec4", &
|
||||
"3a3c9f3eb0e0f96edc0", &
|
||||
"705919ffb636f96b390", &
|
||||
"43daaaa8163d6bc2bd4", &
|
||||
"96e11ea798b74b10e98", &
|
||||
"811150609c9dee8230c", &
|
||||
"be713f85ab34380f4b0", &
|
||||
"5a02c4abaaccb8f24c4", &
|
||||
"67bdebb8863d04768cc", &
|
||||
"5a449cd90c3dbdfe844", &
|
||||
"9c7a54d1c4ef7418b84", &
|
||||
"cd82fefaaf9cd28cd8c", &
|
||||
"ca47e847fabb0054a38", &
|
||||
"f0b30cef6aab9e37f98", &
|
||||
"d948d912fbcc1708710", &
|
||||
"cce1a7b355053d98270", &
|
||||
"4cf227c225a9063dd48", &
|
||||
"2db92612e9ba1418e24", &
|
||||
"3d215c04c762c3d6a28", &
|
||||
"77de65500b5624ceb0c", &
|
||||
"fd1a1df99ded2fb9d88", &
|
||||
"2a19392c71438410fb8", &
|
||||
"a9b486a9d26ed579754", &
|
||||
"b698d244ac78d97a498", &
|
||||
"3d7975b74d727a5e704", &
|
||||
"38094225a2bce0e1940", &
|
||||
"3d3e58fae40fac342b0", &
|
||||
"7732e839a066e337714", &
|
||||
"69356c082b7753a47b0", &
|
||||
"3e868a55dc403a802ac", &
|
||||
"a0157a14a6bf7fdbbcc", &
|
||||
"1ab628e11a7ab4a7c44", &
|
||||
"9da3a2247d7449052f4", &
|
||||
"199a8a7b114816b97f4", &
|
||||
"b1c5cde2542061704cc", &
|
||||
"432fa8d3a153eafbdc8", &
|
||||
"c4ece7e400d8a89c448", &
|
||||
"316ecf74e4b983f007c", &
|
||||
"6a14fa8e713bb5e8adc", &
|
||||
"da4b957ded8374e3640", &
|
||||
"0a804dba7c7e4533300", &
|
||||
"52c342ed033f86580e0", &
|
||||
"1667da8d6fcf4272470", &
|
||||
"da2f7038d550fa88d8c", &
|
||||
"685bcbab1d9dd2c2a44", &
|
||||
"4c93008b3156b3636bc", &
|
||||
"726998d6327ac797c3c", &
|
||||
"44ece7e400d8a8dc448", &
|
||||
"01f9add00dfe823a948", &
|
||||
"dbb95f5ce9e371ad720", &
|
||||
"fc746ee5c76827a8728", &
|
||||
"b25408029506467f4b4", &
|
||||
"9b5c9219e21126b7cf8", &
|
||||
"39ae9f48ba9d1a24f04", &
|
||||
"7de2699623eb507f938", &
|
||||
"b9c6e903ee91dd32934", &
|
||||
"397510d2c6cb5e81de8", &
|
||||
"20157a14aebf7fdbbec", &
|
||||
"067f76ea5817a465980", &
|
||||
"9248f3cea0869feb994", &
|
||||
"23cde2678004ebe5f80", &
|
||||
"5b81fe6848f58e3cfa8", &
|
||||
"a9099ace96bff092904", &
|
||||
"4afa4b0802b33215438", &
|
||||
"f4f740396b030360858", &
|
||||
"fc613f77a35ee1163b8", &
|
||||
"1a4dc27d7e8cc835ff4", &
|
||||
"e9b056f153b39def7ec", &
|
||||
"b62eb777a2f953c7efc", &
|
||||
"388ae4de514b62d238c", &
|
||||
"891529af40e85317160", &
|
||||
"474f1afeb724dbd2ba8", &
|
||||
"11d70880fd88fdd307c", &
|
||||
"29f26a3acb76e6a517c", &
|
||||
"df3e902ff9cadcf776c", &
|
||||
"e3c42da8445965c09f0", &
|
||||
"ce277a6aeccc316dc58", &
|
||||
"4d7841fb71543abd9b8", &
|
||||
"e63230d2d465fb44750", &
|
||||
"b6e11fa798b74b14e98", &
|
||||
"05f189d37c5616547b4", &
|
||||
"ebdb51a81d1e883baa8", &
|
||||
"bf5bc736663bcd53ae0", &
|
||||
"2f8d1cc0936142c08fc", &
|
||||
"436b22fc36d917b6928", &
|
||||
"044b482822ccf2b5f58", &
|
||||
"37b2e839a066e3b7714", &
|
||||
"2a9b4b765c581f0c51c", &
|
||||
"10a7d44cecf8e6628dc", &
|
||||
"ad95f02df6d5502dd4c", &
|
||||
"bbd34f8afd63deaf564", &
|
||||
"cabddfeb01fce632788", &
|
||||
"66b57babeedd6124114", &
|
||||
"7813e0454fbd462be8c", &
|
||||
"b6105ed6f01ea621d04", &
|
||||
"9f68bbcec679d1c088c", &
|
||||
"673da96e414fc7a0f40", &
|
||||
"5568adb935e11084abc", &
|
||||
"f6dd308de5e5c4f6fb0", &
|
||||
"3b49e80d40ae596c7b4", &
|
||||
"a3cde2478004ebe5f80", &
|
||||
"dd8e4f309e919d5ed94", &
|
||||
"5a4020d387757d7bc28", &
|
||||
"64f9e02ae32362a255c", &
|
||||
"630d5942d392334b0dc", &
|
||||
"0bd7e9f4229b2dee210", &
|
||||
"bca549a9467d3a2550c", &
|
||||
"2fef7b1f578c5e28d04", &
|
||||
"f35e0fdda1be4b3b35c", &
|
||||
"69ed575e7cc537d2394", &
|
||||
"7dfdcfbfd5ef3093680", &
|
||||
"b3b2921af97f251d328", &
|
||||
"5622d0fe90363522364", &
|
||||
"fcd4fc7fa04a69d2ac4", &
|
||||
"1119ea451502ed9ab34", &
|
||||
"970ee777ec969a41754", &
|
||||
"688d14f8afec76783dc", &
|
||||
"4d0b8a1028578407420", &
|
||||
"d3d2138d9fa268da3e8", &
|
||||
"df1bdbff898e006394c", &
|
||||
"8ac478a916bb0b77684", &
|
||||
"93881997428e2c17a94", &
|
||||
"4aa510e746245e90c08", &
|
||||
"e00cb8543f85a5d58b8", &
|
||||
"9100d8eb74031073044", &
|
||||
"38710e4235bd1e4003c", &
|
||||
"6aef311cac4c4dccfd4", &
|
||||
"58430f577f51c36b3e0", &
|
||||
"12082fa5d4268a95b4c", &
|
||||
"7a7435a0aca071e64d0", &
|
||||
"cd8250ebadc95de15b0", &
|
||||
"debad40c852e99d64dc", &
|
||||
"4e6caa5e7c86efef748", &
|
||||
"a5d4cbb97e726e3c580", &
|
||||
"7e3a0a2c73ef8553640", &
|
||||
"b60bfc2fd2bd8f530dc", &
|
||||
"32dbef097a5f84b0318", &
|
||||
"4cc7c1cf434300be380", &
|
||||
"896840945be8eabf7f0", &
|
||||
"36c9b10ec694819a0a0", &
|
||||
"349f46a799ef95a47c8", &
|
||||
"9bdcd4ce2563e560b74", &
|
||||
"b19fcd7111a335c52ec"/
|
||||
|
||||
@@ -0,0 +1,423 @@
|
||||
data Mn/ &
|
||||
84, 101, 144, &
|
||||
10, 14, 138, &
|
||||
87, 148, 166, &
|
||||
1, 50, 67, &
|
||||
2, 53, 74, &
|
||||
3, 83, 113, &
|
||||
4, 90, 121, &
|
||||
5, 63, 128, &
|
||||
6, 124, 138, &
|
||||
8, 22, 108, &
|
||||
11, 28, 159, &
|
||||
12, 18, 142, &
|
||||
13, 24, 145, &
|
||||
15, 131, 149, &
|
||||
16, 44, 93, &
|
||||
17, 41, 47, &
|
||||
19, 37, 129, &
|
||||
20, 33, 94, &
|
||||
21, 100, 154, &
|
||||
23, 71, 141, &
|
||||
25, 89, 95, &
|
||||
26, 105, 153, &
|
||||
27, 36, 58, &
|
||||
29, 59, 166, &
|
||||
30, 52, 126, &
|
||||
31, 61, 77, &
|
||||
32, 84, 111, &
|
||||
34, 97, 155, &
|
||||
38, 98, 127, &
|
||||
39, 76, 143, &
|
||||
40, 55, 92, &
|
||||
42, 147, 158, &
|
||||
43, 82, 148, &
|
||||
45, 49, 109, &
|
||||
46, 70, 86, &
|
||||
48, 78, 139, &
|
||||
51, 101, 104, &
|
||||
54, 63, 96, &
|
||||
56, 81, 125, &
|
||||
57, 117, 164, &
|
||||
60, 75, 107, &
|
||||
39, 62, 132, &
|
||||
64, 110, 118, &
|
||||
24, 65, 146, &
|
||||
66, 80, 134, &
|
||||
68, 91, 114, &
|
||||
69, 123, 162, &
|
||||
72, 88, 152, &
|
||||
79, 99, 130, &
|
||||
85, 112, 124, &
|
||||
99, 103, 157, &
|
||||
106, 115, 133, &
|
||||
116, 120, 140, &
|
||||
119, 161, 165, &
|
||||
64, 122, 137, &
|
||||
34, 89, 135, &
|
||||
136, 138, 163, &
|
||||
93, 144, 159, &
|
||||
35, 130, 150, &
|
||||
62, 151, 164, &
|
||||
104, 153, 160, &
|
||||
1, 106, 166, &
|
||||
2, 132, 152, &
|
||||
3, 11, 105, &
|
||||
4, 18, 160, &
|
||||
5, 53, 91, &
|
||||
6, 109, 141, &
|
||||
7, 111, 113, &
|
||||
8, 54, 136, &
|
||||
9, 61, 92, &
|
||||
10, 40, 101, &
|
||||
12, 30, 146, &
|
||||
13, 37, 82, &
|
||||
14, 29, 95, &
|
||||
1, 47, 131, &
|
||||
2, 8, 139, &
|
||||
3, 58, 130, &
|
||||
4, 96, 115, &
|
||||
5, 119, 129, &
|
||||
6, 60, 148, &
|
||||
7, 95, 163, &
|
||||
2, 35, 56, &
|
||||
9, 67, 79, &
|
||||
10, 75, 122, &
|
||||
11, 17, 121, &
|
||||
12, 137, 145, &
|
||||
13, 36, 152, &
|
||||
14, 15, 155, &
|
||||
15, 134, 143, &
|
||||
16, 106, 125, &
|
||||
11, 106, 157, &
|
||||
18, 99, 118, &
|
||||
19, 50, 94, &
|
||||
20, 126, 158, &
|
||||
21, 41, 135, &
|
||||
22, 24, 71, &
|
||||
23, 42, 136, &
|
||||
22, 109, 161, &
|
||||
25, 39, 46, &
|
||||
26, 45, 55, &
|
||||
27, 77, 82, &
|
||||
28, 73, 166, &
|
||||
29, 69, 76, &
|
||||
30, 108, 150, &
|
||||
31, 91, 146, &
|
||||
14, 32, 147, &
|
||||
33, 35, 107, &
|
||||
34, 103, 111, &
|
||||
8, 94, 122, &
|
||||
13, 70, 151, &
|
||||
32, 37, 142, &
|
||||
3, 38, 87, &
|
||||
25, 51, 92, &
|
||||
40, 57, 72, &
|
||||
21, 108, 153, &
|
||||
23, 26, 142, &
|
||||
43, 44, 48, &
|
||||
30, 43, 62, &
|
||||
7, 45, 154, &
|
||||
16, 46, 149, &
|
||||
1, 53, 75, &
|
||||
33, 44, 160, &
|
||||
49, 86, 157, &
|
||||
19, 80, 159, &
|
||||
51, 116, 138, &
|
||||
52, 92, 98, &
|
||||
6, 12, 47, &
|
||||
54, 83, 101, &
|
||||
24, 55, 102, &
|
||||
56, 63, 120, &
|
||||
17, 57, 82, &
|
||||
38, 154, 162, &
|
||||
59, 74, 151, &
|
||||
53, 144, 164, &
|
||||
61, 85, 117, &
|
||||
62, 66, 90, &
|
||||
48, 113, 145, &
|
||||
64, 65, 128, &
|
||||
27, 29, 65, &
|
||||
58, 63, 134, &
|
||||
9, 74, 83, &
|
||||
68, 109, 113, &
|
||||
41, 61, 69, &
|
||||
36, 60, 155, &
|
||||
42, 64, 144, &
|
||||
40, 90, 130, &
|
||||
28, 110, 135, &
|
||||
20, 59, 112, &
|
||||
70, 110, 124, &
|
||||
54, 76, 105, &
|
||||
4, 77, 111, &
|
||||
78, 104, 143, &
|
||||
66, 67, 91, &
|
||||
80, 81, 88, &
|
||||
50, 101, 132, &
|
||||
71, 97, 120, &
|
||||
72, 131, 158, &
|
||||
84, 133, 141, &
|
||||
5, 85, 99, &
|
||||
49, 89, 133, &
|
||||
87, 132, 140, &
|
||||
34, 88, 104, &
|
||||
89, 105, 147, &
|
||||
6, 76, 102, &
|
||||
18, 31, 163, &
|
||||
52, 96, 140, &
|
||||
93, 102, 165, &
|
||||
79, 104, 165, &
|
||||
81, 100, 126, &
|
||||
95, 121, 152, &
|
||||
97, 123, 153, &
|
||||
37, 98, 114, &
|
||||
8, 91, 155, &
|
||||
100, 114, 160, &
|
||||
2, 26, 28, &
|
||||
93, 116, 150, &
|
||||
68, 103, 166, &
|
||||
78, 117, 125, &
|
||||
86, 107, 127, &
|
||||
4, 59, 136, &
|
||||
9, 37, 97, &
|
||||
7, 30, 75, &
|
||||
80, 148, 153, &
|
||||
73, 138, 164, &
|
||||
10, 39, 103, &
|
||||
39, 146, 156, &
|
||||
48, 129, 136, &
|
||||
5, 17, 51, &
|
||||
112, 149, 161, &
|
||||
11, 24, 126, &
|
||||
1, 70, 78, &
|
||||
14, 113, 118, &
|
||||
10, 119, 141, &
|
||||
13, 33, 105, &
|
||||
19, 57, 89, &
|
||||
12, 25, 56, &
|
||||
16, 18, 54, &
|
||||
84, 124, 162, &
|
||||
20, 41, 134, &
|
||||
15, 45, 82, &
|
||||
115, 118, 123, &
|
||||
128, 139, 149, &
|
||||
127, 156, 159, &
|
||||
21, 141, 152, &
|
||||
23, 130, 156, &
|
||||
3, 160, 164, &
|
||||
22, 90, 110, &
|
||||
35, 61, 109, &
|
||||
31, 87, 158, &
|
||||
42, 60, 106, &
|
||||
137, 140, 157, &
|
||||
27, 114, 124, &
|
||||
32, 62, 125, &
|
||||
34, 38, 128, &
|
||||
40, 123, 139, &
|
||||
29, 66, 86, &
|
||||
36, 52, 161, &
|
||||
43, 63, 133, &
|
||||
46, 73, 108, &
|
||||
44, 135, 146, &
|
||||
47, 115, 127, &
|
||||
49, 74, 116, &
|
||||
58, 102, 122, &
|
||||
55, 85, 132, &
|
||||
50, 65, 150, &
|
||||
67, 145, 162, &
|
||||
53, 71, 77, &
|
||||
69, 88, 142, &
|
||||
68, 72, 93, &
|
||||
9, 64, 95, &
|
||||
92, 94, 111, &
|
||||
81, 83, 119, &
|
||||
98, 143, 163, &
|
||||
73, 79, 96, &
|
||||
35, 129, 131, &
|
||||
99, 100, 151, &
|
||||
7, 112, 159, &
|
||||
117, 137, 156, &
|
||||
120, 147, 154, &
|
||||
107, 121, 165/
|
||||
|
||||
data Nm/ &
|
||||
4, 62, 75, 121, 191, &
|
||||
5, 63, 76, 82, 175, &
|
||||
6, 64, 77, 112, 206, &
|
||||
7, 65, 78, 151, 180, &
|
||||
8, 66, 79, 159, 188, &
|
||||
9, 67, 80, 127, 164, &
|
||||
68, 81, 119, 182, 237, &
|
||||
10, 69, 76, 109, 173, &
|
||||
70, 83, 141, 181, 230, &
|
||||
2, 71, 84, 185, 193, &
|
||||
11, 64, 85, 91, 190, &
|
||||
12, 72, 86, 127, 196, &
|
||||
13, 73, 87, 110, 194, &
|
||||
2, 74, 88, 106, 192, &
|
||||
14, 88, 89, 200, 0, &
|
||||
15, 90, 120, 197, 0, &
|
||||
16, 85, 131, 188, 0, &
|
||||
12, 65, 92, 165, 197, &
|
||||
17, 93, 124, 195, 0, &
|
||||
18, 94, 148, 199, 0, &
|
||||
19, 95, 115, 204, 0, &
|
||||
10, 96, 98, 207, 0, &
|
||||
20, 97, 116, 205, 0, &
|
||||
13, 44, 96, 129, 190, &
|
||||
21, 99, 113, 196, 0, &
|
||||
22, 100, 116, 175, 0, &
|
||||
23, 101, 139, 212, 0, &
|
||||
11, 102, 147, 175, 0, &
|
||||
24, 74, 103, 139, 216, &
|
||||
25, 72, 104, 118, 182, &
|
||||
26, 105, 165, 209, 0, &
|
||||
27, 106, 111, 213, 0, &
|
||||
18, 107, 122, 194, 0, &
|
||||
28, 56, 108, 162, 214, &
|
||||
59, 82, 107, 208, 235, &
|
||||
23, 87, 144, 217, 0, &
|
||||
17, 73, 111, 172, 181, &
|
||||
29, 112, 132, 214, 0, &
|
||||
30, 42, 99, 185, 186, &
|
||||
31, 71, 114, 146, 215, &
|
||||
16, 95, 143, 199, 0, &
|
||||
32, 97, 145, 210, 0, &
|
||||
33, 117, 118, 218, 0, &
|
||||
15, 117, 122, 220, 0, &
|
||||
34, 100, 119, 200, 0, &
|
||||
35, 99, 120, 219, 0, &
|
||||
16, 75, 127, 221, 0, &
|
||||
36, 117, 137, 187, 0, &
|
||||
34, 123, 160, 222, 0, &
|
||||
4, 93, 155, 225, 0, &
|
||||
37, 113, 125, 188, 0, &
|
||||
25, 126, 166, 217, 0, &
|
||||
5, 66, 121, 134, 227, &
|
||||
38, 69, 128, 150, 197, &
|
||||
31, 100, 129, 224, 0, &
|
||||
39, 82, 130, 196, 0, &
|
||||
40, 114, 131, 195, 0, &
|
||||
23, 77, 140, 223, 0, &
|
||||
24, 133, 148, 180, 0, &
|
||||
41, 80, 144, 210, 0, &
|
||||
26, 70, 135, 143, 208, &
|
||||
42, 60, 118, 136, 213, &
|
||||
8, 38, 130, 140, 218, &
|
||||
43, 55, 138, 145, 230, &
|
||||
44, 138, 139, 225, 0, &
|
||||
45, 136, 153, 216, 0, &
|
||||
4, 83, 153, 226, 0, &
|
||||
46, 142, 177, 229, 0, &
|
||||
47, 103, 143, 228, 0, &
|
||||
35, 110, 149, 191, 0, &
|
||||
20, 96, 156, 227, 0, &
|
||||
48, 114, 157, 229, 0, &
|
||||
102, 184, 219, 234, 0, &
|
||||
5, 133, 141, 222, 0, &
|
||||
41, 84, 121, 182, 0, &
|
||||
30, 103, 150, 164, 0, &
|
||||
26, 101, 151, 227, 0, &
|
||||
36, 152, 178, 191, 0, &
|
||||
49, 83, 168, 234, 0, &
|
||||
45, 124, 154, 183, 0, &
|
||||
39, 154, 169, 232, 0, &
|
||||
33, 73, 101, 131, 200, &
|
||||
6, 128, 141, 232, 0, &
|
||||
1, 27, 158, 198, 0, &
|
||||
50, 135, 159, 224, 0, &
|
||||
35, 123, 179, 216, 0, &
|
||||
3, 112, 161, 209, 0, &
|
||||
48, 154, 162, 228, 0, &
|
||||
21, 56, 160, 163, 195, &
|
||||
7, 136, 146, 207, 0, &
|
||||
46, 66, 105, 153, 173, &
|
||||
31, 70, 113, 126, 231, &
|
||||
15, 58, 167, 176, 229, &
|
||||
18, 93, 109, 231, 0, &
|
||||
21, 74, 81, 170, 230, &
|
||||
38, 78, 166, 234, 0, &
|
||||
28, 156, 171, 181, 0, &
|
||||
29, 126, 172, 233, 0, &
|
||||
49, 51, 92, 159, 236, &
|
||||
19, 169, 174, 236, 0, &
|
||||
1, 37, 71, 128, 155, &
|
||||
129, 164, 167, 223, 0, &
|
||||
51, 108, 177, 185, 0, &
|
||||
37, 61, 152, 162, 168, &
|
||||
22, 64, 150, 163, 194, &
|
||||
52, 62, 90, 91, 210, &
|
||||
41, 107, 179, 240, 0, &
|
||||
10, 104, 115, 219, 0, &
|
||||
34, 67, 98, 142, 208, &
|
||||
43, 147, 149, 207, 0, &
|
||||
27, 68, 108, 151, 231, &
|
||||
50, 148, 189, 237, 0, &
|
||||
6, 68, 137, 142, 192, &
|
||||
46, 172, 174, 212, 0, &
|
||||
52, 78, 201, 221, 0, &
|
||||
53, 125, 176, 222, 0, &
|
||||
40, 135, 178, 238, 0, &
|
||||
43, 92, 192, 201, 0, &
|
||||
54, 79, 193, 232, 0, &
|
||||
53, 130, 156, 239, 0, &
|
||||
7, 85, 170, 240, 0, &
|
||||
55, 84, 109, 223, 0, &
|
||||
47, 171, 201, 215, 0, &
|
||||
9, 50, 149, 198, 212, &
|
||||
39, 90, 178, 213, 0, &
|
||||
25, 94, 169, 190, 0, &
|
||||
29, 179, 203, 221, 0, &
|
||||
8, 138, 202, 214, 0, &
|
||||
17, 79, 187, 235, 0, &
|
||||
49, 59, 77, 146, 205, &
|
||||
14, 75, 157, 235, 0, &
|
||||
42, 63, 155, 161, 224, &
|
||||
52, 158, 160, 218, 0, &
|
||||
45, 89, 140, 199, 0, &
|
||||
56, 95, 147, 220, 0, &
|
||||
57, 69, 97, 180, 187, &
|
||||
55, 86, 211, 238, 0, &
|
||||
2, 9, 57, 125, 184, &
|
||||
36, 76, 202, 215, 0, &
|
||||
53, 161, 166, 211, 0, &
|
||||
20, 67, 158, 193, 204, &
|
||||
12, 111, 116, 228, 0, &
|
||||
30, 89, 152, 233, 0, &
|
||||
1, 58, 134, 145, 0, &
|
||||
13, 86, 137, 226, 0, &
|
||||
44, 72, 105, 186, 220, &
|
||||
32, 106, 163, 239, 0, &
|
||||
3, 33, 80, 183, 0, &
|
||||
14, 120, 189, 202, 0, &
|
||||
59, 104, 176, 225, 0, &
|
||||
60, 110, 133, 236, 0, &
|
||||
48, 63, 87, 170, 204, &
|
||||
22, 61, 115, 171, 183, &
|
||||
19, 119, 132, 239, 0, &
|
||||
28, 88, 144, 173, 0, &
|
||||
186, 203, 205, 238, 0, &
|
||||
51, 91, 123, 211, 0, &
|
||||
32, 94, 157, 209, 0, &
|
||||
11, 58, 124, 203, 237, &
|
||||
61, 65, 122, 174, 206, &
|
||||
54, 98, 189, 217, 0, &
|
||||
47, 132, 198, 226, 0, &
|
||||
57, 81, 165, 233, 0, &
|
||||
40, 60, 134, 184, 206, &
|
||||
54, 167, 168, 240, 0, &
|
||||
3, 24, 62, 102, 177/
|
||||
|
||||
data nrw/ &
|
||||
5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,5,4,4, &
|
||||
4,4,4,5,4,4,4,4,5,5,4,4,4,5,5,4,5,4,5,5, &
|
||||
4,4,4,4,4,4,4,4,4,4,4,4,5,5,4,4,4,4,4,4, &
|
||||
5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, &
|
||||
4,5,4,4,4,4,4,4,5,4,5,5,5,4,5,4,4,4,5,4, &
|
||||
5,4,4,5,5,5,4,4,5,4,5,4,5,4,4,4,4,4,4,4, &
|
||||
4,4,4,5,4,4,4,4,4,5,4,5,4,4,4,5,4,5,4,4, &
|
||||
5,4,4,4,4,5,4,4,4,4,4,5,5,4,4,4,4,4,5,5, &
|
||||
4,4,4,5,4,5/
|
||||
|
||||
ncw=3
|
||||
|
||||
@@ -0,0 +1,137 @@
|
||||
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)
|
||||
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,*) norder
|
||||
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(*,*) "norder : ",norder
|
||||
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
|
||||
dmin=0.0
|
||||
maxosd=2
|
||||
call decode240_101(llr, Keff, maxosd, norder, apmask, message101, cw, ntype, nharderror, dmin)
|
||||
if(nharderror.ge.0) then
|
||||
n2err=0
|
||||
do i=1,N
|
||||
if( cw(i).ne.codeword(i) ) 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
|
||||
@@ -0,0 +1,125 @@
|
||||
program ldpcsim240_74
|
||||
|
||||
! End-to-end test of the (240,74)/crc24 encoder and decoders.
|
||||
|
||||
use packjt77
|
||||
|
||||
parameter(N=240, K=74, M=N-K)
|
||||
character*8 arg
|
||||
character*37 msg0,msg
|
||||
character*77 c77
|
||||
character*24 c24
|
||||
integer*1 msgbits(74)
|
||||
integer*1 apmask(240)
|
||||
integer*1 cw(240)
|
||||
integer*1 codeword(N),message74(74)
|
||||
integer ncrc24
|
||||
real rxdata(N),llr(N)
|
||||
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_74 20 5 1000 0.85 64 "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 [50,74]'
|
||||
print*,'WSPR-format message is optional'
|
||||
return
|
||||
endif
|
||||
call getarg(1,arg)
|
||||
read(arg,*) max_iterations
|
||||
call getarg(2,arg)
|
||||
read(arg,*) norder
|
||||
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(*,*) "norder : ",norder
|
||||
write(*,*) "s : ",s
|
||||
write(*,*) "K : ",Keff
|
||||
|
||||
msgbits=0
|
||||
read(c77,'(50i1)') msgbits(1:50)
|
||||
write(*,*) 'message'
|
||||
write(*,'(50i1)') msgbits(1:50)
|
||||
|
||||
call get_crc24(msgbits,74,ncrc24)
|
||||
write(c24,'(b24.24)') ncrc24
|
||||
read(c24,'(24i1)') msgbits(51:74)
|
||||
write(*,'(24i1)') msgbits(51:74)
|
||||
write(*,*) 'message with crc24'
|
||||
write(*,'(74i1)') msgbits(1:74)
|
||||
call encode240_74(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
|
||||
dmin=0.0
|
||||
maxosd=0
|
||||
call decode240_74(llr, Keff, maxosd, norder, apmask, message74, cw, ntype, nharderror, dmin)
|
||||
if(nharderror.ge.0) then
|
||||
n2err=0
|
||||
do i=1,N
|
||||
if( cw(i).ne.codeword(i) ) 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,e10.3)") db,esn0,ngood,nue,pberr
|
||||
|
||||
enddo
|
||||
|
||||
end program ldpcsim240_74
|
||||
@@ -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=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 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
|
||||
|
||||
@@ -0,0 +1,403 @@
|
||||
subroutine osd240_74(llr,k,apmask,ndeep,message74,cw,nhardmin,dmin)
|
||||
!
|
||||
! An ordered-statistics decoder for the (240,74) code.
|
||||
! Message payload is 50 bits. Any or all of a 24-bit CRC can be
|
||||
! used for detecting incorrect codewords. The remaining CRC bits are
|
||||
! cascaded with the LDPC code for the purpose of improving the
|
||||
! distance spectrum of the code.
|
||||
!
|
||||
! If p1 (0.le.p1.le.24) is the number of CRC24 bits that are
|
||||
! to be used for bad codeword detection, then the argument k should
|
||||
! be set to 77+p1.
|
||||
!
|
||||
! Valid values for k are in the range [50,74].
|
||||
!
|
||||
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 message74(74)
|
||||
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=74-k and p1+p2=24.
|
||||
!
|
||||
! The last p2 bits of the CRC24 are cascaded with the LDPC code.
|
||||
!
|
||||
! The first p1=k-50 CRC24 bits will be used for error detection.
|
||||
!
|
||||
allocate( gen(k,N) )
|
||||
gen=0
|
||||
do i=1,k
|
||||
message74=0
|
||||
message74(i)=1
|
||||
if(i.le.50) then
|
||||
call get_crc24(message74,74,ncrc24)
|
||||
write(c24,'(b24.24)') ncrc24
|
||||
read(c24,'(24i1)') message74(51:74)
|
||||
message74(51:k)=0
|
||||
endif
|
||||
call encode240_74(message74,cw)
|
||||
gen(i,:)=cw
|
||||
enddo
|
||||
|
||||
first=.false.
|
||||
endif
|
||||
|
||||
rx=llr
|
||||
apmaskr=apmask
|
||||
|
||||
! Hard decisions on the received word.
|
||||
hdec=0
|
||||
where(rx .ge. 0) hdec=1
|
||||
|
||||
! Use magnitude of received symbols as a measure of reliability.
|
||||
absrx=abs(rx)
|
||||
call indexx(absrx,N,indx)
|
||||
|
||||
! Re-order the columns of the generator matrix in order of decreasing reliability.
|
||||
do i=1,N
|
||||
genmrb(1:k,i)=gen(1:k,indx(N+1-i))
|
||||
indices(i)=indx(N+1-i)
|
||||
enddo
|
||||
|
||||
! Do gaussian elimination to create a generator matrix with the most reliable
|
||||
! received bits in positions 1:k in order of decreasing reliability (more or less).
|
||||
do id=1,k ! diagonal element indices
|
||||
do icol=id,k+20 ! The 20 is ad hoc - beware
|
||||
iflag=0
|
||||
if( genmrb(id,icol) .eq. 1 ) then
|
||||
iflag=1
|
||||
if( icol .ne. id ) then ! reorder column
|
||||
temp(1:k)=genmrb(1:k,id)
|
||||
genmrb(1:k,id)=genmrb(1:k,icol)
|
||||
genmrb(1:k,icol)=temp(1:k)
|
||||
itmp=indices(id)
|
||||
indices(id)=indices(icol)
|
||||
indices(icol)=itmp
|
||||
endif
|
||||
do ii=1,k
|
||||
if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then
|
||||
genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N))
|
||||
endif
|
||||
enddo
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
g2=transpose(genmrb)
|
||||
|
||||
! The hard decisions for the k MRB bits define the order 0 message, m0.
|
||||
! Encode m0 using the modified generator matrix to find the "order 0" codeword.
|
||||
! Flip various combinations of bits in m0 and re-encode to generate a list of
|
||||
! codewords. Return the member of the list that has the smallest Euclidean
|
||||
! distance to the received word.
|
||||
|
||||
hdec=hdec(indices) ! hard decisions from received symbols
|
||||
m0=hdec(1:k) ! zero'th order message
|
||||
absrx=absrx(indices)
|
||||
rx=rx(indices)
|
||||
apmaskr=apmaskr(indices)
|
||||
|
||||
call mrbencode74(m0,c0,g2,N,k)
|
||||
nxor=ieor(c0,hdec)
|
||||
nhardmin=sum(nxor)
|
||||
dmin=sum(nxor*absrx)
|
||||
|
||||
cw=c0
|
||||
ntotal=0
|
||||
nrejected=0
|
||||
npre1=0
|
||||
npre2=0
|
||||
|
||||
if(ndeep.eq.0) goto 998 ! norder=0
|
||||
if(ndeep.gt.6) ndeep=6
|
||||
if( ndeep.eq. 1) then
|
||||
nord=1
|
||||
npre1=0
|
||||
npre2=0
|
||||
nt=40
|
||||
ntheta=12
|
||||
elseif(ndeep.eq.2) then
|
||||
nord=1
|
||||
npre1=1
|
||||
npre2=0
|
||||
nt=40
|
||||
ntheta=12
|
||||
elseif(ndeep.eq.3) then
|
||||
nord=1
|
||||
npre1=1
|
||||
npre2=1
|
||||
nt=40
|
||||
ntheta=12
|
||||
ntau=14
|
||||
elseif(ndeep.eq.4) then
|
||||
nord=2
|
||||
npre1=1
|
||||
npre2=1
|
||||
nt=40
|
||||
ntheta=12
|
||||
ntau=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 mrbencode74(me,ce,g2,N,k)
|
||||
e2sub=ieor(ce(k+1:N),hdec(k+1:N))
|
||||
e2=e2sub
|
||||
nd1kpt=sum(e2sub(1:nt))+1
|
||||
d1=sum(ieor(me(1:k),hdec(1:k))*absrx(1:k))
|
||||
else
|
||||
e2=ieor(e2sub,g2(k+1:N,n1))
|
||||
nd1kpt=sum(e2(1:nt))+2
|
||||
endif
|
||||
if(nd1kpt .le. ntheta) then
|
||||
call mrbencode74(me,ce,g2,N,k)
|
||||
nxor=ieor(ce,hdec)
|
||||
if(n1.eq.iflag) then
|
||||
dd=d1+sum(e2sub*absrx(k+1:N))
|
||||
else
|
||||
dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(k+1:N))
|
||||
endif
|
||||
if( dd .lt. dmin ) then
|
||||
dmin=dd
|
||||
cw=ce
|
||||
nhardmin=sum(nxor)
|
||||
nd1kptbest=nd1kpt
|
||||
endif
|
||||
else
|
||||
nrejected=nrejected+1
|
||||
endif
|
||||
enddo
|
||||
! Get the next test error pattern, iflag will go negative
|
||||
! when the last pattern with weight iorder has been generated.
|
||||
call nextpat74(misub,k,iorder,iflag)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if(npre2.eq.1) then
|
||||
reset=.true.
|
||||
ntotal=0
|
||||
do i1=k,1,-1
|
||||
do i2=i1-1,1,-1
|
||||
ntotal=ntotal+1
|
||||
mi(1:ntau)=ieor(g2(k+1:k+ntau,i1),g2(k+1:k+ntau,i2))
|
||||
call boxit74(reset,mi(1:ntau),ntau,ntotal,i1,i2)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
ncount2=0
|
||||
ntotal2=0
|
||||
reset=.true.
|
||||
! Now run through again and do the second pre-processing rule
|
||||
misub(1:k-nord)=0
|
||||
misub(k-nord+1:k)=1
|
||||
iflag=k-nord+1
|
||||
do while(iflag .ge.0)
|
||||
me=ieor(m0,misub)
|
||||
call mrbencode74(me,ce,g2,N,k)
|
||||
e2sub=ieor(ce(k+1:N),hdec(k+1:N))
|
||||
do i2=0,ntau
|
||||
ntotal2=ntotal2+1
|
||||
ui=0
|
||||
if(i2.gt.0) ui(i2)=1
|
||||
r2pat=ieor(e2sub,ui)
|
||||
778 continue
|
||||
call fetchit74(reset,r2pat(1:ntau),ntau,in1,in2)
|
||||
if(in1.gt.0.and.in2.gt.0) then
|
||||
ncount2=ncount2+1
|
||||
mi=misub
|
||||
mi(in1)=1
|
||||
mi(in2)=1
|
||||
if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:k),mi).eq.1)) cycle
|
||||
me=ieor(m0,mi)
|
||||
call mrbencode74(me,ce,g2,N,k)
|
||||
nxor=ieor(ce,hdec)
|
||||
dd=sum(nxor*absrx)
|
||||
if( dd .lt. dmin ) then
|
||||
dmin=dd
|
||||
cw=ce
|
||||
nhardmin=sum(nxor)
|
||||
endif
|
||||
goto 778
|
||||
endif
|
||||
enddo
|
||||
call nextpat74(misub,k,nord,iflag)
|
||||
enddo
|
||||
endif
|
||||
|
||||
998 continue
|
||||
! Re-order the codeword to [message bits][parity bits] format.
|
||||
cw(indices)=cw
|
||||
hdec(indices)=hdec
|
||||
message74=cw(1:74)
|
||||
call get_crc24(message74,74,nbadcrc)
|
||||
if(nbadcrc.ne.0) nhardmin=-nhardmin
|
||||
|
||||
return
|
||||
end subroutine osd240_74
|
||||
|
||||
subroutine mrbencode74(me,codeword,g2,N,K)
|
||||
integer*1 me(K),codeword(N),g2(N,K)
|
||||
! fast encoding for low-weight test patterns
|
||||
codeword=0
|
||||
do i=1,K
|
||||
if( me(i) .eq. 1 ) then
|
||||
codeword=ieor(codeword,g2(1:N,i))
|
||||
endif
|
||||
enddo
|
||||
return
|
||||
end subroutine mrbencode74
|
||||
|
||||
subroutine nextpat74(mi,k,iorder,iflag)
|
||||
integer*1 mi(k),ms(k)
|
||||
! generate the next test error pattern
|
||||
ind=-1
|
||||
do i=1,k-1
|
||||
if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i
|
||||
enddo
|
||||
if( ind .lt. 0 ) then ! no more patterns of this order
|
||||
iflag=ind
|
||||
return
|
||||
endif
|
||||
ms=0
|
||||
ms(1:ind-1)=mi(1:ind-1)
|
||||
ms(ind)=1
|
||||
ms(ind+1)=0
|
||||
if( ind+1 .lt. k ) then
|
||||
nz=iorder-sum(ms)
|
||||
ms(k-nz+1:k)=1
|
||||
endif
|
||||
mi=ms
|
||||
do i=1,k ! iflag will point to the lowest-index 1 in mi
|
||||
if(mi(i).eq.1) then
|
||||
iflag=i
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
return
|
||||
end subroutine nextpat74
|
||||
|
||||
subroutine boxit74(reset,e2,ntau,npindex,i1,i2)
|
||||
integer*1 e2(1:ntau)
|
||||
integer indexes(5000,2),fp(0:525000),np(5000)
|
||||
logical reset
|
||||
common/boxes/indexes,fp,np
|
||||
|
||||
if(reset) then
|
||||
patterns=-1
|
||||
fp=-1
|
||||
np=-1
|
||||
sc=-1
|
||||
indexes=-1
|
||||
reset=.false.
|
||||
endif
|
||||
|
||||
indexes(npindex,1)=i1
|
||||
indexes(npindex,2)=i2
|
||||
ipat=0
|
||||
do i=1,ntau
|
||||
if(e2(i).eq.1) then
|
||||
ipat=ipat+ishft(1,ntau-i)
|
||||
endif
|
||||
enddo
|
||||
|
||||
ip=fp(ipat) ! see what's currently stored in fp(ipat)
|
||||
if(ip.eq.-1) then
|
||||
fp(ipat)=npindex
|
||||
else
|
||||
do while (np(ip).ne.-1)
|
||||
ip=np(ip)
|
||||
enddo
|
||||
np(ip)=npindex
|
||||
endif
|
||||
return
|
||||
end subroutine boxit74
|
||||
|
||||
subroutine fetchit74(reset,e2,ntau,i1,i2)
|
||||
integer indexes(5000,2),fp(0:525000),np(5000)
|
||||
integer lastpat
|
||||
integer*1 e2(ntau)
|
||||
logical reset
|
||||
common/boxes/indexes,fp,np
|
||||
save lastpat,inext
|
||||
|
||||
if(reset) then
|
||||
lastpat=-1
|
||||
reset=.false.
|
||||
endif
|
||||
|
||||
ipat=0
|
||||
do i=1,ntau
|
||||
if(e2(i).eq.1) then
|
||||
ipat=ipat+ishft(1,ntau-i)
|
||||
endif
|
||||
enddo
|
||||
index=fp(ipat)
|
||||
|
||||
if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices
|
||||
i1=indexes(index,1)
|
||||
i2=indexes(index,2)
|
||||
inext=np(index)
|
||||
elseif(lastpat.eq.ipat .and. inext.gt.0) then
|
||||
i1=indexes(inext,1)
|
||||
i2=indexes(inext,2)
|
||||
inext=np(inext)
|
||||
else
|
||||
i1=-1
|
||||
i2=-1
|
||||
inext=-1
|
||||
endif
|
||||
lastpat=ipat
|
||||
return
|
||||
end subroutine fetchit74
|
||||
|
||||
Reference in New Issue
Block a user