WSJT-X/lib/superfox/qpc_sync.f90

102 lines
2.1 KiB
Fortran
Raw Normal View History

subroutine qpc_sync(crcvd0,fsample,isync,fsync,ftol,f2,t2,snrsync)
parameter(N9SEC=9*12000,NMAX=15*12000,NDOWN=16,NZ=N9SEC/NDOWN)
complex crcvd0(NMAX) !Signal as received
complex c0(0:N9SEC-1) !For long FFT
complex c1(0:NZ-1)
complex c1sum(0:NZ-1)
complex z
real s(N9SEC/4)
real p(-1125:1125)
integer ipk(1)
integer isync(24)
baud=12000.0/1024.0
df2=fsample/N9SEC
fac=1.0/N9SEC
c0=fac*crcvd0(1:N9SEC)
call four2a(c0,N9SEC,1,-1,1) !Forward c2c FFT
iz=N9SEC/4
do i=1,iz
s(i)=real(c0(i))**2 + aimag(c0(i))**2
enddo
do i=1,4 !Smooth the spectrum a bit
call smo121(s,iz)
enddo
ia=nint((fsync-ftol)/df2)
ib=nint((fsync+ftol)/df2)
ipk=maxloc(s(ia:ib))
i0=ipk(1) + ia - 1
f2=df2*i0-750.0 ! f2 is the offset from nominal 750 Hz.
ia=nint(i0-baud/df2)
ib=nint(i0+baud/df2)
s1=0.0
s0=0.0
do i=ia,ib
s0=s0+s(i)
s1=s1+(i-i0)*s(i)
enddo
delta=s1/s0
i0=nint(i0+delta)
f2=i0*df2-750.0
c1=0.
ia=nint(i0-baud/df2)
ib=nint(i0+baud/df2)
do i=ia,ib
j=i-i0
if(j.ge.0) c1(j)=c0(i)
if(j.lt.0) c1(j+NZ)=c0(i)
enddo
call four2a(c1,NZ,1,1,1) !Reverse c2c FFT: back to time domain
c1sum(0)=c1(0)
do i=1,NZ-1
c1sum(i)=c1sum(i-1) + c1(i)
enddo
nspsd=1024/NDOWN
dt=NDOWN/12000.0
lagmax=1.5/dt
i0=nint(0.5*fsample/NDOWN) !Nominal start time is 0.5 s
pmax=0.
lagpk=0
do lag=-lagmax,lagmax
sp=0.
do j=1,24
i1=i0 + (isync(j)-1)*nspsd + lag
i2=i1 + nspsd
if(i1.lt.0 .or. i1.gt.NZ-1) cycle
if(i2.lt.0 .or. i2.gt.NZ-1) cycle
z=c1sum(i2)-c1sum(i1)
sp=sp + real(z)**2 + aimag(z)**2
enddo
if(sp.gt.pmax) then
pmax=sp
lagpk=lag
endif
p(lag)=sp
enddo
t2=lagpk*dt
snrsync=0.
sp=0.
sq=0.
nsum=0
tsym=1024/12000.0
do lag=-lagmax,lagmax
t=(lag-lagpk)*dt
if(abs(t).lt.tsym) cycle
nsum=nsum+1
sp=sp + p(lag)
sq=sq + p(lag)*p(lag)
enddo
ave=sp/nsum
rms=sqrt(sq/nsum-ave*ave)
snrsync=(pmax-ave)/rms
return
end subroutine qpc_sync