WSJT-X/lib/fst240_decode.f90

800 lines
28 KiB
Fortran
Raw Normal View History

2020-06-27 09:53:11 -04:00
module fst240_decode
type :: fst240_decoder
procedure(fst240_decode_callback), pointer :: callback
2020-06-27 09:53:11 -04:00
contains
procedure :: decode
end type fst240_decoder
abstract interface
subroutine fst240_decode_callback (this,nutc,sync,nsnr,dt,freq, &
decoded,nap,qual,ntrperiod,lwspr)
import fst240_decoder
implicit none
class(fst240_decoder), intent(inout) :: this
integer, intent(in) :: nutc
real, intent(in) :: sync
integer, intent(in) :: nsnr
real, intent(in) :: dt
real, intent(in) :: freq
character(len=37), intent(in) :: decoded
integer, intent(in) :: nap
real, intent(in) :: qual
integer, intent(in) :: ntrperiod
logical, intent(in) :: lwspr
end subroutine fst240_decode_callback
end interface
2020-06-27 09:53:11 -04:00
contains
subroutine decode(this,callback,iwave,nutc,nQSOProgress,nfqso, &
nfa,nfb,nsubmode,ndepth,ntrperiod,nexp_decode,ntol,nzhsym, &
emedelay,lapcqonly,napwid,mycall,hiscall,nfsplit,iwspr)
2020-06-27 09:53:11 -04:00
use timer_module, only: timer
use packjt77
2020-07-07 21:25:25 -04:00
use, intrinsic :: iso_c_binding
include 'fst240/fst240_params.f90'
parameter (MAXCAND=100)
class(fst240_decoder), intent(inout) :: this
procedure(fst240_decode_callback) :: callback
character*37 decodes(100)
character*37 msg,msgsent
character*77 c77
character*12 mycall,hiscall
character*12 mycall0,hiscall0
complex, allocatable :: c2(:)
complex, allocatable :: cframe(:)
complex, allocatable :: c_bigfft(:) !Complex waveform
2020-07-07 21:25:25 -04:00
real, allocatable, target :: r_data(:)
complex, pointer, dimension(:) :: c_data_ptr
real llr(240),llra(240),llrb(240),llrc(240),llrd(240)
real candidates(100,4)
real bitmetrics(320,4)
real s4(0:3,NN)
real minsync
logical lapcqonly
integer itone(NN)
integer hmod
integer*1 apmask(240),cw(240)
integer*1 hbits(320)
integer*1 message101(101),message74(74),message77(77)
integer*1 rvec(77)
integer apbits(240)
integer nappasses(0:5) ! # of decoding passes for QSO states 0-5
integer naptypes(0:5,4) ! (nQSOProgress,decoding pass)
integer mcq(29),mrrr(19),m73(19),mrr73(19)
logical badsync,unpk77_success,single_decode
logical first,nohiscall,lwspr
integer*2 iwave(300*12000)
data mcq/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0/
data mrrr/0,1,1,1,1,1,1,0,1,0,0,1,0,0,1,0,0,0,1/
data m73/0,1,1,1,1,1,1,0,1,0,0,1,0,1,0,0,0,0,1/
data mrr73/0,1,1,1,1,1,1,0,0,1,1,1,0,1,0,1,0,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/
data first/.true./
save first,apbits,nappasses,naptypes,mycall0,hiscall0
this%callback => callback
dxcall13=hiscall ! initialize for use in packjt77
mycall13=mycall
fMHz=1.0
if(first) then
mcq=2*mod(mcq+rvec(1:29),2)-1
mrrr=2*mod(mrrr+rvec(59:77),2)-1
m73=2*mod(m73+rvec(59:77),2)-1
mrr73=2*mod(mrr73+rvec(59:77),2)-1
nappasses(0)=2
nappasses(1)=2
nappasses(2)=2
nappasses(3)=2
nappasses(4)=2
nappasses(5)=3
! iaptype
!------------------------
! 1 CQ ??? ??? (29 ap bits)
! 2 MyCall ??? ??? (29 ap bits)
! 3 MyCall DxCall ??? (58 ap bits)
! 4 MyCall DxCall RRR (77 ap bits)
! 5 MyCall DxCall 73 (77 ap bits)
! 6 MyCall DxCall RR73 (77 ap bits)
!********
naptypes(0,1:4)=(/1,2,0,0/) ! Tx6 selected (CQ)
naptypes(1,1:4)=(/2,3,0,0/) ! Tx1
naptypes(2,1:4)=(/2,3,0,0/) ! Tx2
naptypes(3,1:4)=(/3,6,0,0/) ! Tx3
naptypes(4,1:4)=(/3,6,0,0/) ! Tx4
naptypes(5,1:4)=(/3,1,2,0/) ! Tx5
mycall0=''
hiscall0=''
first=.false.
endif
l1=index(mycall,char(0))
if(l1.ne.0) mycall(l1:)=" "
l1=index(hiscall,char(0))
if(l1.ne.0) hiscall(l1:)=" "
if(mycall.ne.mycall0 .or. hiscall.ne.hiscall0) then
apbits=0
apbits(1)=99
apbits(30)=99
if(len(trim(mycall)) .lt. 3) go to 10
nohiscall=.false.
hiscall0=hiscall
if(len(trim(hiscall0)).lt.3) then
hiscall0=mycall ! use mycall for dummy hiscall - mycall won't be hashed.
nohiscall=.true.
endif
msg=trim(mycall)//' '//trim(hiscall0)//' RR73'
i3=-1
n3=-1
call pack77(msg,i3,n3,c77)
call unpack77(c77,1,msgsent,unpk77_success)
if(i3.ne.1 .or. (msg.ne.msgsent) .or. .not.unpk77_success) go to 10
read(c77,'(77i1)') message77
message77=mod(message77+rvec,2)
call encode174_91(message77,cw)
apbits=2*cw-1
if(nohiscall) apbits(30)=99
10 continue
mycall0=mycall
hiscall0=hiscall
endif
!************************************
hmod=2**nsubmode
if(nfqso+nqsoprogress.eq.-999) return
Keff=91
2020-06-27 09:53:11 -04:00
nmax=15*12000
single_decode=iand(nexp_decode,32).eq.32
if(ntrperiod.eq.15) then
nsps=720
nmax=15*12000
ndown=18/hmod !nss=40,80,160,400
if(hmod.eq.4) ndown=4
if(hmod.eq.8) ndown=2
nfft1=int(nmax/ndown)*ndown
else if(ntrperiod.eq.30) then
nsps=1680
nmax=30*12000
ndown=42/hmod !nss=40,80,168,336
nfft1=359856
if(hmod.eq.4) then
ndown=10
nfft1=nmax
endif
if(hmod.eq.8) then
ndown=5
nfft1=nmax
endif
else if(ntrperiod.eq.60) then
nsps=3888
nmax=60*12000
ndown=96/hmod !nss=36,81,162,324
if(hmod.eq.1) ndown=108
nfft1=int(719808/ndown)*ndown
else if(ntrperiod.eq.120) then
nsps=8200
nmax=120*12000
ndown=200/hmod !nss=40,82,164,328
if(hmod.eq.1) ndown=205
nfft1=int(nmax/ndown)*ndown
else if(ntrperiod.eq.300) then
nsps=21504
nmax=300*12000
ndown=512/hmod !nss=42,84,168,336
nfft1=int((nmax-200)/ndown)*ndown
end if
nss=nsps/ndown
fs=12000.0 !Sample rate
fs2=fs/ndown
nspsec=nint(fs2)
dt=1.0/fs !Sample interval (s)
dt2=1.0/fs2
tt=nsps*dt !Duration of "itone" symbols (s)
baud=1.0/tt
sigbw=4.0*hmod*baud
nfft2=nfft1/ndown !make sure that nfft1 is exactly nfft2*ndown
nfft1=nfft2*ndown
nh1=nfft1/2
allocate( r_data(1:nfft1+2) )
2020-07-07 21:25:25 -04:00
call c_f_pointer (c_loc (r_data), c_data_ptr, [(nfft1+2)/2]) ! c_data_ptr shares memory with r_data
allocate( c_bigfft(0:nfft1/2) )
allocate( c2(0:nfft2-1) )
allocate( cframe(0:160*nss-1) )
if(ndepth.eq.3) then
nblock=4
if(hmod.eq.1) nblock=4 ! number of block sizes to try
jittermax=2
norder=3
elseif(ndepth.eq.2) then
nblock=1
if(hmod.eq.1) nblock=3
jittermax=0
norder=3
elseif(ndepth.eq.1) then
nblock=1
jittermax=0
norder=3
endif
! The big fft is done once and is used for calculating the smoothed spectrum
2020-06-27 09:53:11 -04:00
! and also for downconverting/downsampling each candidate.
r_data(1:nfft1)=iwave(1:nfft1)
r_data(nfft1+1:nfft1+2)=0.0
2020-07-07 21:25:25 -04:00
call four2a(c_data_ptr,nfft1,1,-1,0)
c_bigfft=cmplx(r_data(1:nfft1+2:2),r_data(2:nfft1+2:2))
! write(*,3001) iwspr,nfa,nfb,nfsplit,ndepth
!3001 format('a',5i5)
! iwspr=1 !### For hardwired tests ###
if(iwspr.eq.0) then
itype1=1
itype2=1
elseif( iwspr.eq.1 ) then
itype1=2
itype2=2
elseif( iwspr.eq.2 ) then
itype1=1
itype2=2
endif
if(hmod.eq.1) then
if(fMHz.lt.2.0) then
nsyncoh=8 ! Use N=8 for sync
nhicoh=1 ! Use N=1,2,4,8 for symbol estimation
else
nsyncoh=4 ! Use N=4 for sync
nhicoh=0 ! Use N=1,2,3,4 for symbol estimation
endif
else
if(hmod.eq.2) nsyncoh=1
if(hmod.eq.4) nsyncoh=1
if(hmod.eq.8) nsyncoh=1
endif
do iqorw=itype1,itype2 ! iqorw=1 for QSO mode and iqorw=2 for wspr-type messages
if( iwspr.lt.2 ) then
if( single_decode ) then
fa=max(100,nint(nfqso+1.5*hmod*baud-ntol))
fb=min(4800,nint(nfqso+1.5*hmod*baud+ntol))
else
fa=max(100,nfa)
fb=min(4800,nfb)
endif
elseif( iwspr.eq.2 .and. iqorw.eq.1 ) then
fa=max(100,nfa)
fb=nfsplit
elseif( iwspr.eq.2 .and. iqorw.eq.2 ) then
fa=nfsplit
fb=min(4800,nfb)
endif
if(hmod.eq.1) then
if(ntrperiod.eq.15) minsync=1.15
if(ntrperiod.gt.15) minsync=1.20
elseif(hmod.gt.1) then
minsync=1.5
endif
2020-06-27 09:53:11 -04:00
! Get first approximation of candidate frequencies
call get_candidates_fst240(c_bigfft,nfft1,nsps,hmod,fs,fa,fb, &
minsync,ncand,candidates,base)
2020-06-27 09:53:11 -04:00
ndecodes=0
decodes=' '
2020-06-27 09:53:11 -04:00
isbest=0
fc2=0.
do icand=1,ncand
fc0=candidates(icand,1)
detmet=candidates(icand,2)
2020-06-27 09:53:11 -04:00
! Downconvert and downsample a slice of the spectrum centered on the
! rough estimate of the candidates frequency.
! Output array c2 is complex baseband sampled at 12000/ndown Sa/sec.
! The size of the downsampled c2 array is nfft2=nfft1/ndown
call fst240_downsample(c_bigfft,nfft1,ndown,fc0,sigbw,c2)
call timer('sync240 ',0)
do isync=0,1
if(isync.eq.0) then
fc1=0.0
if(emedelay.lt.0.1) then ! search offsets from 0 s to 2 s
is0=1.5*nspsec
ishw=1.5*nspsec
else ! search plus or minus 1.5 s centered on emedelay
is0=nint(emedelay*nspsec)
ishw=1.5*nspsec
endif
isst=4*hmod
ifhw=12
df=.1*baud
else if(isync.eq.1) then
fc1=fc2
is0=isbest
ishw=4*hmod
isst=1*hmod
ifhw=7
df=.02*baud
2020-06-27 09:53:11 -04:00
endif
smax=0.0
do if=-ifhw,ifhw
fc=fc1+df*if
do istart=max(1,is0-ishw),is0+ishw,isst
call sync_fst240(c2,istart,fc,hmod,nsyncoh,nfft2,nss,fs2,sync)
if(sync.gt.smax) then
fc2=fc
isbest=istart
smax=sync
endif
enddo
enddo
2020-06-27 09:53:11 -04:00
enddo
call timer('sync240 ',1)
2020-06-27 09:53:11 -04:00
fc_synced = fc0 + fc2
dt_synced = (isbest-fs2)*dt2 !nominal dt is 1 second so frame starts at sample fs2
candidates(icand,3)=fc_synced
candidates(icand,4)=isbest
enddo
2020-06-27 09:53:11 -04:00
! remove duplicate candidates
do icand=1,ncand
fc=candidates(icand,3)
isbest=nint(candidates(icand,4))
do ic2=1,ncand
fc2=candidates(ic2,3)
isbest2=nint(candidates(ic2,4))
if(ic2.ne.icand .and. fc2.gt.0.0) then
if(abs(fc2-fc).lt.0.10*baud) then ! same frequency
if(abs(isbest2-isbest).le.2) then
candidates(ic2,3)=-1
endif
endif
2020-06-27 09:53:11 -04:00
endif
enddo
enddo
ic=0
do icand=1,ncand
if(candidates(icand,3).gt.0) then
ic=ic+1
candidates(ic,:)=candidates(icand,:)
endif
enddo
ncand=ic
xsnr=0.
do icand=1,ncand
sync=candidates(icand,2)
fc_synced=candidates(icand,3)
isbest=nint(candidates(icand,4))
xdt=(isbest-nspsec)/fs2
if(ntrperiod.eq.15) xdt=(isbest-real(nspsec)/2.0)/fs2
call fst240_downsample(c_bigfft,nfft1,ndown,fc_synced,sigbw,c2)
do ijitter=0,jittermax
if(ijitter.eq.0) ioffset=0
if(ijitter.eq.1) ioffset=1
if(ijitter.eq.2) ioffset=-1
is0=isbest+ioffset
if(is0.lt.0) cycle
cframe=c2(is0:is0+160*nss-1)
bitmetrics=0
if(hmod.eq.1) then
call get_fst240_bitmetrics(cframe,nss,hmod,nblock,nhicoh,bitmetrics,s4,badsync)
else
call get_fst240_bitmetrics2(cframe,nss,hmod,nblock,bitmetrics,s4,badsync)
endif
if(badsync) cycle
hbits=0
where(bitmetrics(:,1).ge.0) hbits=1
ns1=count(hbits( 1: 16).eq.(/0,0,0,1,1,0,1,1,0,1,0,0,1,1,1,0/))
ns2=count(hbits( 77: 92).eq.(/1,1,1,0,0,1,0,0,1,0,1,1,0,0,0,1/))
ns3=count(hbits(153:168).eq.(/0,0,0,1,1,0,1,1,0,1,0,0,1,1,1,0/))
ns4=count(hbits(229:244).eq.(/1,1,1,0,0,1,0,0,1,0,1,1,0,0,0,1/))
ns5=count(hbits(305:320).eq.(/0,0,0,1,1,0,1,1,0,1,0,0,1,1,1,0/))
nsync_qual=ns1+ns2+ns3+ns4+ns5
! if(nsync_qual.lt. 46) cycle !### Value ?? ###
scalefac=2.83
llra( 1: 60)=bitmetrics( 17: 76, 1)
llra( 61:120)=bitmetrics( 93:152, 1)
llra(121:180)=bitmetrics(169:228, 1)
llra(181:240)=bitmetrics(245:304, 1)
llra=scalefac*llra
llrb( 1: 60)=bitmetrics( 17: 76, 2)
llrb( 61:120)=bitmetrics( 93:152, 2)
llrb(121:180)=bitmetrics(169:228, 2)
llrb(181:240)=bitmetrics(245:304, 2)
llrb=scalefac*llrb
llrc( 1: 60)=bitmetrics( 17: 76, 3)
llrc( 61:120)=bitmetrics( 93:152, 3)
llrc(121:180)=bitmetrics(169:228, 3)
llrc(181:240)=bitmetrics(245:304, 3)
llrc=scalefac*llrc
llrd( 1: 60)=bitmetrics( 17: 76, 4)
llrd( 61:120)=bitmetrics( 93:152, 4)
llrd(121:180)=bitmetrics(169:228, 4)
llrd(181:240)=bitmetrics(245:304, 4)
llrd=scalefac*llrd
apmag=maxval(abs(llra))*1.1
ntmax=nblock+nappasses(nQSOProgress)
if(lapcqonly) ntmax=nblock+1
if(ndepth.eq.1) ntmax=nblock
apmask=0
if(iqorw.eq.2) then ! 50-bit msgs, no ap decoding
nblock=4
ntmax=nblock
endif
do itry=1,ntmax
if(itry.eq.1) llr=llra
if(itry.eq.2.and.itry.le.nblock) llr=llrb
if(itry.eq.3.and.itry.le.nblock) llr=llrc
if(itry.eq.4.and.itry.le.nblock) llr=llrd
if(itry.le.nblock) then
apmask=0
iaptype=0
endif
napwid=1.2*(4.0*baud*hmod)
if(itry.gt.nblock) then
if(nblock.eq.1) llr=llra
if(nblock.gt.1) llr=llrc
iaptype=naptypes(nQSOProgress,itry-nblock)
if(lapcqonly) iaptype=1
if(iaptype.ge.2 .and. apbits(1).gt.1) cycle ! No, or nonstandard, mycall
if(iaptype.ge.3 .and. apbits(30).gt.1) cycle ! No, or nonstandard, dxcall
if(iaptype.eq.1) then ! CQ
apmask=0
apmask(1:29)=1
llr(1:29)=apmag*mcq(1:29)
endif
if(iaptype.eq.2) then ! MyCall ??? ???
apmask=0
apmask(1:29)=1
llr(1:29)=apmag*apbits(1:29)
2020-06-27 09:53:11 -04:00
endif
if(iaptype.eq.3) then ! MyCall DxCall ???
apmask=0
apmask(1:58)=1
llr(1:58)=apmag*apbits(1:58)
2020-06-27 09:53:11 -04:00
endif
if(iaptype.eq.4 .or. iaptype.eq.5 .or. iaptype .eq.6) then
apmask=0
apmask(1:77)=1
llr(1:58)=apmag*apbits(1:58)
if(iaptype.eq.4) llr(59:77)=apmag*mrrr(1:19)
if(iaptype.eq.5) llr(59:77)=apmag*m73(1:19)
if(iaptype.eq.6) llr(59:77)=apmag*mrr73(1:19)
endif
endif
dmin=0.0
nharderrors=-1
unpk77_success=.false.
if(iqorw.eq.1) then
maxosd=2
2020-06-30 16:50:50 -04:00
Keff=91
norder=3
call timer('d240_101',0)
call decode240_101(llr,Keff,maxosd,norder,apmask,message101, &
cw,ntype,nharderrors,dmin)
call timer('d240_101',1)
elseif(iqorw.eq.2) then
maxosd=2
call timer('d240_74 ',0)
Keff=64
norder=4
call decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw, &
ntype,nharderrors,dmin)
call timer('d240_74 ',1)
2020-06-27 09:53:11 -04:00
endif
2020-06-30 16:50:50 -04:00
if(nharderrors .ge.0) then
2020-07-01 17:42:22 -04:00
if(count(cw.eq.1).eq.0) then
nharderrors=-nharderrors
cycle
endif
2020-06-30 16:50:50 -04:00
if(iqorw.eq.1) then
write(c77,'(77i1)') mod(message101(1:77)+rvec,2)
call unpack77(c77,0,msg,unpk77_success)
else
write(c77,'(50i1)') message74(1:50)
c77(51:77)='000000000000000000000110000'
call unpack77(c77,0,msg,unpk77_success)
endif
if(unpk77_success) then
idupe=0
do i=1,ndecodes
if(decodes(i).eq.msg) idupe=1
enddo
if(idupe.eq.1) goto 2002
ndecodes=ndecodes+1
decodes(ndecodes)=msg
2020-06-30 17:05:52 -04:00
if(iqorw.eq.1) then
call get_fst240_tones_from_bits(message101,itone,0)
else
call get_fst240_tones_from_bits(message74,itone,1)
endif
xsig=0
do i=1,NN
xsig=xsig+s4(itone(i),i)**2
enddo
arg=400.0*(xsig/base)-1.0
if(arg.gt.0.0) then
xsnr=10*log10(arg)-21.0-11.7*log10(nsps/800.0)
else
xsnr=-99.9
endif
else
cycle
endif
nsnr=nint(xsnr)
qual=0.
fsig=fc_synced - 1.5*hmod*baud
2020-07-04 14:06:50 -04:00
! write(21,'(i6,8i6,f7.1,f9.2,f7.1,1x,f7.2,1x,f7.1,1x,a37)') &
! nutc,icand,itry,nsyncoh,iaptype,ijitter,ntype,nsync_qual,nharderrors,dmin,sync,xsnr,xdt,fsig,msg
! flush(21)
call this%callback(nutc,smax1,nsnr,xdt,fsig,msg, &
iaptype,qual,ntrperiod,lwspr)
goto 2002
endif
enddo ! metrics
enddo ! istart jitter
2020-06-30 16:50:50 -04:00
2002 enddo !candidate list
enddo ! iqorw
2020-06-27 09:53:11 -04:00
return
end subroutine decode
2020-06-27 09:53:11 -04:00
subroutine sync_fst240(cd0,i0,f0,hmod,ncoh,np,nss,fs,sync)
2020-06-27 09:53:11 -04:00
! Compute sync power for a complex, downsampled FST240 signal.
include 'fst240/fst240_params.f90'
complex cd0(0:np-1)
complex, allocatable, save :: csync1(:),csync2(:)
complex, allocatable, save :: csynct1(:),csynct2(:)
complex ctwk(8*nss)
complex z1,z2,z3,z4,z5
logical first
integer hmod,isyncword1(0:7),isyncword2(0:7)
real f0save
data isyncword1/0,1,3,2,1,0,2,3/
data isyncword2/2,3,1,0,3,2,0,1/
data first/.true./,f0save/-99.9/,nss0/-1/
save first,twopi,dt,fac,f0save,nss0
p(z1)=(real(z1*fac)**2 + aimag(z1*fac)**2)**0.5 !Compute power
if(nss.ne.nss0 .and. allocated(csync1)) deallocate(csync1,csync2,csynct1,csynct2)
if(first .or. nss.ne.nss0) then
allocate( csync1(8*nss), csync2(8*nss) )
allocate( csynct1(8*nss), csynct2(8*nss) )
twopi=8.0*atan(1.0)
dt=1/fs
k=1
phi1=0.0
phi2=0.0
do i=0,7
dphi1=twopi*hmod*(isyncword1(i)-1.5)/real(nss)
dphi2=twopi*hmod*(isyncword2(i)-1.5)/real(nss)
do j=1,nss
csync1(k)=cmplx(cos(phi1),sin(phi1))
csync2(k)=cmplx(cos(phi2),sin(phi2))
phi1=mod(phi1+dphi1,twopi)
phi2=mod(phi2+dphi2,twopi)
k=k+1
enddo
2020-06-27 09:53:11 -04:00
enddo
first=.false.
nss0=nss
fac=1.0/(8.0*nss)
2020-06-27 09:53:11 -04:00
endif
if(f0.ne.f0save) then
dphi=twopi*f0*dt
phi=0.0
do i=1,8*nss
ctwk(i)=cmplx(cos(phi),sin(phi))
phi=mod(phi+dphi,twopi)
enddo
csynct1=ctwk*csync1
csynct2=ctwk*csync2
f0save=f0
2020-06-27 09:53:11 -04:00
endif
i1=i0 !Costas arrays
i2=i0+38*nss
i3=i0+76*nss
i4=i0+114*nss
i5=i0+152*nss
s1=0.0
s2=0.0
s3=0.0
s4=0.0
s5=0.0
if(ncoh.gt.0) then
nsec=8/ncoh
do i=1,nsec
is=(i-1)*ncoh*nss
z1=0
if(i1+is.ge.1) then
z1=sum(cd0(i1+is:i1+is+ncoh*nss-1)*conjg(csynct1(is+1:is+ncoh*nss)))
endif
z2=sum(cd0(i2+is:i2+is+ncoh*nss-1)*conjg(csynct2(is+1:is+ncoh*nss)))
z3=sum(cd0(i3+is:i3+is+ncoh*nss-1)*conjg(csynct1(is+1:is+ncoh*nss)))
z4=sum(cd0(i4+is:i4+is+ncoh*nss-1)*conjg(csynct2(is+1:is+ncoh*nss)))
z5=0
if(i5+is+ncoh*nss-1.le.np) then
z5=sum(cd0(i5+is:i5+is+ncoh*nss-1)*conjg(csynct1(is+1:is+ncoh*nss)))
endif
s1=s1+abs(z1)/(8*nss)
s2=s2+abs(z2)/(8*nss)
s3=s3+abs(z3)/(8*nss)
s4=s4+abs(z4)/(8*nss)
s5=s5+abs(z5)/(8*nss)
enddo
else
nsub=-ncoh
nps=nss/nsub
do i=1,8
do isub=1,nsub
is=(i-1)*nss+(isub-1)*nps
if(i1+is.ge.1) then
s1=s1+abs(sum(cd0(i1+is:i1+is+nps-1)*conjg(csynct1(is+1:is+nps))))
endif
s2=s2+abs(sum(cd0(i2+is:i2+is+nps-1)*conjg(csynct1(is+1:is+nps))))
s3=s3+abs(sum(cd0(i3+is:i3+is+nps-1)*conjg(csynct1(is+1:is+nps))))
s4=s4+abs(sum(cd0(i4+is:i4+is+nps-1)*conjg(csynct1(is+1:is+nps))))
s5=0
if(i5+is+ncoh*nss-1.le.np) then
s5=s5+abs(sum(cd0(i5+is:i5+is+nps-1)*conjg(csynct1(is+1:is+nps))))
endif
enddo
enddo
endif
sync = s1+s2+s3+s4+s5
return
end subroutine sync_fst240
2020-06-27 09:53:11 -04:00
subroutine fst240_downsample(c_bigfft,nfft1,ndown,f0,sigbw,c1)
2020-06-27 09:53:11 -04:00
! Output: Complex data in c(), sampled at 12000/ndown Hz
complex c_bigfft(0:nfft1/2)
complex c1(0:nfft1/ndown-1)
df=12000.0/nfft1
i0=nint(f0/df)
ih=nint( ( f0 + 1.3*sigbw/2.0 )/df)
nbw=ih-i0+1
c1=0.
c1(0)=c_bigfft(i0)
nfft2=nfft1/ndown
do i=1,nbw
if(i0+i.le.nfft1/2) c1(i)=c_bigfft(i0+i)
if(i0-i.ge.0) c1(nfft2-i)=c_bigfft(i0-i)
2020-06-27 09:53:11 -04:00
enddo
c1=c1/nfft2
call four2a(c1,nfft2,1,1,1) !c2c FFT back to time domain
return
end subroutine fst240_downsample
subroutine get_candidates_fst240(c_bigfft,nfft1,nsps,hmod,fs,fa,fb, &
minsync,ncand,candidates,base)
complex c_bigfft(0:nfft1/2) !Full length FFT of raw data
integer hmod !Modulation index (submode)
integer im(1) !For maxloc
real candidates(100,4) !Candidate list
real s(18000) !Low resolution power spectrum
real s2(18000) !CCF of s() with 4 tones
real xdb(-3:3) !Model 4-tone CCF peaks
real minsync
data xdb/0.25,0.50,0.75,1.0,0.75,0.50,0.25/
nh1=nfft1/2
df1=fs/nfft1
baud=fs/nsps !Keying rate
df2=baud/2.0
nd=df2/df1 !s() sums this many bins of big FFT
ndh=nd/2
ia=nint(max(100.0,fa)/df2) !Low frequency search limit
ib=nint(min(4800.0,fb)/df2) !High frequency limit
signal_bw=4*(12000.0/nsps)*hmod
analysis_bw=min(4800.0,fb)-max(100.0,fa)
xnoise_bw=10.0*signal_bw !Is this a good compromise?
if(analysis_bw.gt.xnoise_bw) then
ina=ia
inb=ib
else
fcenter=(fa+fb)/2.0 !If noise_bw > analysis_bw,
fl = max(100.0,fcenter-xnoise_bw/2.)/df2 !we'll search over noise_bw
fh = min(4800.0,fcenter+xnoise_bw/2.)/df2
ina=nint(fl)
inb=nint(fh)
endif
s=0. !Compute low-resloution power spectrum
do i=ina,inb ! noise analysis window includes signal analysis window
j0=nint(i*df2/df1)
do j=j0-ndh,j0+ndh
s(i)=s(i) + real(c_bigfft(j))**2 + aimag(c_bigfft(j))**2
enddo
enddo
ina=max(ina,1+3*hmod) !Don't run off the ends
inb=min(inb,18000-3*hmod)
s2=0.
do i=ina,inb !Compute CCF of s() and 4 tones
s2(i)=s(i-hmod*3) + s(i-hmod) +s(i+hmod) +s(i+hmod*3)
enddo
call pctile(s2(ina+hmod*3:inb-hmod*3),inb-ina+1-hmod*6,30,base)
s2=s2/base !Normalize wrt noise level
ncand=0
candidates=0
if(ia.lt.3) ia=3
if(ib.gt.18000-2) ib=18000-2
! Find candidates, using the CLEAN algorithm to remove a model of each one
! from s2() after it has been found.
pval=99.99
do while(ncand.lt.100)
im=maxloc(s2(ia:ib))
iploc=ia+im(1)-1 !Index of CCF peak
pval=s2(iploc) !Peak value
if(pval.lt.minsync) exit
if(s2(iploc).gt.thresh) then !Is this a possible candidate?
do i=-3,+3 !Remove 0.9 of a model CCF at
k=iploc+2*hmod*i !this frequency from s2()
if(k.ge.ia .and. k.le.ib) then
s2(k)=max(0.,s2(k)-0.9*pval*xdb(i))
endif
enddo
ncand=ncand+1
candidates(ncand,1)=df2*iploc !Candidate frequency
candidates(ncand,2)=pval !Rough estimate of SNR
endif
enddo
2020-06-27 09:53:11 -04:00
return
end subroutine get_candidates_fst240
2020-06-27 09:53:11 -04:00
end module fst240_decode