Do a final sync peakup in frequency using longer FFTs.

lib/superfox/sfox_sync.f90

@@ -5,13 +5,12 @@ subroutine sfox_sync(iwave,fsample,isync,f,t)
   integer*2 iwave(0:NMAX-1)
   integer isync(44)
   integer ipeak(2)
+  integer ipeak2(1)
   complex, allocatable :: c(:)             !Work array
-  real x(171)
   real, allocatable :: s(:,:)              !Symbol spectra, stepped by NSTEP 
   real, allocatable :: savg(:)             !Average spectrum
   real, allocatable :: ccf(:,:)
-!  character*1 line(-15:15),mark(0:6),c1
+  real, allocatable :: s2(:)               !Fine spectrum of sync tone
-!  data mark/' ','.','-','+','X','$','#'/
 
   nfft=nsps
   nh=nfft/2
@@ -27,11 +26,6 @@ subroutine sfox_sync(iwave,fsample,isync,f,t)
   lag1=-lagmax
   lag2=lagmax
 
-  x=0.
-  do i=1,NS
-     x(isync(i))=1.0
-  enddo
-
   allocate(s(0:nh/2,jz))
   allocate(savg(0:nh/2))
   allocate(c(0:nfft-1))
@@ -71,12 +65,12 @@ subroutine sfox_sync(iwave,fsample,isync,f,t)
      ccfmax=0.
      do lag=lag1,lag2
         ccft=0.
-        do kk=1,NS
+        do m=1,NS
-           k=isync(kk)
+           k=isync(m)
            n=NSTEP*(k-1) + 1
            j=n+lag+j0
            if(j.ge.1 .and. j.le.jz) ccft=ccft + s(i,j)
-        enddo  ! kk
+        enddo  ! m
         ccft=ccft - NS*savg(i)
         ccf(i,lag)=ccft
         if(ccft.gt.ccfmax) then
@@ -96,34 +90,51 @@ subroutine sfox_sync(iwave,fsample,isync,f,t)
   ipk=ipeak(1)-1+ia
   jpk=ipeak(2)-1+lag1
 
-  dxi=0.
   dxj=0.
-  if(ipk.gt.ia .and. ipk.lt.ib) then
-     call peakup(ccf(ipk-1,jpk),ccf(ipk,jpk),ccf(ipk+1,jpk),dxi)
-  endif
   if(jpk.gt.lag1 .and. jpk.lt.lag2) then
      call peakup(ccf(ipk,jpk-1),ccf(ipk,jpk),ccf(ipk,jpk+1),dxj)
   endif
 
   f=ibest*df + bw/2 + dxi*df
-  t=lagbest*dtstep + dxj*dtstep
+  t=(lagbest+dxj)*dtstep
-!  write(*,4100) ibest,lagbest,f,dxi*df,t,dxj*dtstep
+  t=t-0.01                               !### Why is this needed? ###
-!4100 format(2i6,2f10.1,2f10.3)
 
-  nsum=0
+  nfft2=4*NSPS
-  sq=0.
+  deallocate(c)
-  do lag=lag1,lag2
+  allocate(c(0:nfft2-1))
-     if(abs(lag-lagbest).gt.3) then
+  allocate(s2(0:nfft2-1))
-        sq=sq + ccf(ibest,lag)**2
+
-        nsum=nsum+1
+  i0=(t+0.5)*fsample
-     endif
+  s2=0.
-     write(51,3051) lag*dtstep,ccf(ibest,lag)
+  df2=fsample/nfft2
-3051 format(2f12.4)
+  do m=1,NS
+     i1=i0+(isync(m)-1)*NSPS
+     i2=i1+NSPS-1
+     k=-1
+     do i=i1,i2,2          !Load iwave data into complex array c0, for r2c FFT
+        if(i.gt.0) then
+           xx=iwave(i)
+           yy=iwave(i+1)
+        else
+           xx=0.
+           yy=0.
+        endif
+        k=k+1
+        c(k)=fac*cmplx(xx,yy)
+     enddo
+     c(k+1:)=0.
+     call four2a(c,nfft2,1,-1,0)              !r2c FFT
+     do i=1,nfft2/4
+        s2(i)=s2(i) + real(c(i))**2 + aimag(c(i))**2
+     enddo
   enddo
-
+  ipeak2=maxloc(s2)
-  rms=sqrt(sq/nsum)
+  ipk=ipeak2(1)-1
-  snrsync=ccf(ibest,lagbest)/rms
+  dxi=0.
-!  print*,'snr:',snrsync
+  if(ipk.gt.1 .and. ipk.lt.nfft/4) then
+     call peakup(s2(ipk-1),s2(ipk),s2(ipk+1),dxi)
+  endif
+  f=(ipk+dxi)*df2 + bw/2.0
 
   return
 end subroutine sfox_sync

lib/superfox/sfoxtest.f90

@@ -31,7 +31,6 @@ program sfoxtest
   integer, allocatable :: rxdat2(:)
   integer, allocatable :: rxprob2(:)
   integer, allocatable :: correct(:)
-  logical hard_sync
   character fname*17,arg*12,itu*2
 
   data isync/ 1,   2,   5,  11,  19,  24,  26,  28,  29,  35,  &
@@ -39,13 +38,12 @@ program sfoxtest
              80,  82,  84,  85,  92,  98, 103, 107, 109, 111,  &
             116, 122, 130, 131, 134, 136, 137, 140, 146, 154,  &
             159, 161, 163, 165/
-
   data nsb/1,2,4,7,11,16,22,29,37,39/
 
   nargs=iargc()
   if(nargs.ne.11) then
-     print*,'Usage:   sfoxtest  f0   DT  ITU M  N   K NS v hs nfiles snr'
+     print*,'Usage:   sfoxtest  f0   DT  ITU M  N   K NS v st nfiles snr'
-     print*,'Example: sfoxtest 1500 0.15  MM 7 127 48 33 0  F   10   -10'
+     print*,'Example: sfoxtest 1500 0.15  MM 7 127 48 33 0  3   10   -10'
      print*,'         f0=0 means f0, DT will assume suitable random values'
      print*,'         LQ: Low Latitude Quiet'
      print*,'         MM: Mid Latitude Moderate'
@@ -53,7 +51,7 @@ program sfoxtest
      print*,'         ... and similarly for LM LD MQ MD HQ HM'
      print*,'         NS: number of sync symbols'
      print*,'         v=1 for .wav files, 2 for verbose output, 3 for both'
-     print*,'         hs = T for hard-wired sync'
+     print*,'         st: Sync type, 0 for hard-wired, otherwise 1-3'
      print*,'         snr=0 means loop over SNRs'
      go to 999
   endif
@@ -73,7 +71,7 @@ program sfoxtest
   call getarg(8,arg)
   read(arg,*) nv
   call getarg(9,arg)
-  hard_sync=arg(1:1).eq.'T'
+  read(arg,*) nstype
   call getarg(10,arg)
   read(arg,*) nfiles
   call getarg(11,arg)
@@ -86,8 +84,6 @@ program sfoxtest
   call sfox_init(mm0,nn0,kk0,itu,fspread,delay,fsample,ns0)
   tsync=NSYNC/fsample
   txt=(NN+NS)*NSPS/fsample
-  nstype=nv/10
-  nv=mod(nv,10)
 
   write(*,1000) MM,NN,KK,NSPS,baud,bw,itu,tsync,txt,nstype
 1000 format('M:',i2,'   Base code: (',i3,',',i3,')   NSPS:',i5,   &
@@ -202,13 +198,12 @@ program sfoxtest
         if(snr.lt.90.0) iwave(1:NMAX)=nint(rms*dat(1:NMAX))
 
         crcvd=sig*crcvd+cnoise
-        if(hard_sync) then
+        if(nstype.eq.0) then
-           f=f1  ! + 5.0*(ran1(idum)-0.5)
+           f=f1                                   !Hard-wired sync
-           t=xdt ! + 0.01*(ran1(idum)-0.5)
+           t=xdt
         else
-! Find signal freq and DT
            call timer('sync    ',0)
-           call sfox_sync(iwave,fsample,isync,f,t)
+           call sfox_sync(iwave,fsample,isync,f,t) ! Find signal freq and DT
            call timer('sync    ',1)
         endif
         ferr=f-f1
@@ -223,7 +218,7 @@ program sfoxtest
         endif
 
         a=0.
-        a(1)=1500.0-f - baud
+        a(1)=1500.0-f - baud                !Shift frequencies down by one bin
         call timer('twkfreq ',0)
         call twkfreq(crcvd,crcvd,NMAX,fsample,a)
         call timer('twkfreq ',1)