More work on sync tests, etc., of QRA6[45].

lib/qra/qra64/qra64sim.f90

@@ -14,14 +14,18 @@ program qra64sim
   complex cdat(NMAX)                     !Generated complex waveform
   complex cspread(0:NFFT-1)              !Complex amplitude for Rayleigh fading
   complex z
+  complex c00(0:720000)                  !Analytic signal for dat()
   real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq
-  character msg*22,fname*11,csubmode*1,arg*12
+  character msg*22,fname*11,csubmode*1,arg*12,cd*1
   character msgsent*22
+  logical lsync
+  data lsync/.false./
 
   nargs=iargc()
-  if(nargs.ne. 7) then
+  if(nargs.ne.8) then
-     print *, 'Usage:   qra64sim         "msg"     A-E Nsigs fDop DT Nfiles SNR'
+     print*,'Usage:   qra64sim         "msg"     A-E Nsigs fDop DT Nfiles Sync SNR'
-     print *, 'Example  qra64sim "K1ABC W9XYZ EN37" A   10   0.2 0.0   1     0'
+     print*,'Example  qra64sim "K1ABC W9XYZ EN37" A   10   0.2 0.0   1      T  -26'
+     print*,'Sync = T to include sync test.'
      go to 999
   endif
   call getarg(1,msg)
@@ -36,6 +40,8 @@ program qra64sim
   call getarg(6,arg)
   read(arg,*) nfiles
   call getarg(7,arg)
+  if(arg(1:1).eq.'T' .or. arg(1:1).eq.'1') lsync=.true.
+  call getarg(8,arg)
   read(arg,*) snrdb
 
   if(mode64.ge.8) nsigs=1
@@ -54,6 +60,7 @@ program qra64sim
   write(*,1000) 
 1000 format('File  Sig    Freq  A-E   S/N   DT   Dop    Message'/60('-'))
 
+  nsync=0
   do ifile=1,nfiles                  !Loop over requested number of files
      write(fname,1002) ifile         !Output filename
 1002 format('000000_',i4.4)
@@ -165,6 +172,30 @@ program qra64sim
      if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts))
      write(10) h,iwave(1:npts)                !Save the .wav file
      close(10)
+
+     if(lsync) then
+        cd=' '
+        if(ifile.eq.nfiles) cd='d'
+        nf1=200
+        nf2=3000
+        nfqso=nint(f0)
+        ntol=100
+        minsync=0
+        emedelay=0.0
+        call ana64(dat,npts,c00)
+        call sync64(c00,nf1,nf2,nfqso,ntol,minsync,mode64,emedelay,xdt2,f02,  &
+             jpk0,sync,sync2,width)
+        terr=1.01/(8.0*baud)
+        ferr=1.01*mode64*baud
+        if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1
+        open(40,file='sync64.out',status='unknown',position='append')
+        write(40,1030) ifile,64,csubmode,snrdb,fspread,xdt2-xdt,f02-f0,   &
+             width,sync,sync2,nsync,cd
+1030    format(i4,i3,1x,a1,2f7.1,f7.2,4f8.1,i5,1x,a1)
+        close(40)
+     endif
   enddo
+  if(lsync) write(*,1040) snrdb,nfiles,nsync
+1040 format('SNR:',f6.1,'   nfiles:',i5,'   nsynced:',i5)
 
 999 end program qra64sim

lib/qra/qra65/qra65sim.f90

@@ -14,13 +14,16 @@ program qra65sim
   complex cspread(0:NMAX-1)              !Complex amplitude for Rayleigh fading
   complex z
   real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq
-  character msg*22,fname*17,csubmode*1,arg*12
+  character msg*22,fname*17,csubmode*1,arg*12,cd*1
   character msgsent*22
+  logical lsync
+  data lsync/.false./
   
   nargs=iargc()
-  if(nargs.ne.8) then
+  if(nargs.ne.9) then
-     print *, 'Usage:   qra65sim         "msg"     A-E freq fDop DT TRp Nfiles SNR'
+     print *, 'Usage:   qra65sim         "msg"     A-E freq fDop DT TRp Nfiles Sync SNR'
-     print *, 'Example: qra65sim "K1ABC W9XYZ EN37" A  1500 0.2 0.0  15   1    -10'
+     print *, 'Example: qra65sim "K1ABC W9XYZ EN37" A  1500 0.0 0.0  60   1      T  -26'
+     print*,'Sync = T to include sync test.'
      go to 999
   endif
   call getarg(1,msg)
@@ -37,8 +40,15 @@ program qra65sim
   call getarg(7,arg)
   read(arg,*) nfiles
   call getarg(8,arg)
+  if(arg(1:1).eq.'T' .or. arg(1:1).eq.'1') lsync=.true.
+  call getarg(9,arg)
   read(arg,*) snrdb
 
+  if(nfiles.lt.0) then
+     nfiles=-nfiles
+     lsync=.true.
+  endif
+
   if(ntrperiod.eq.15) then
      nsps=1800
   else if(ntrperiod.eq.30) then
@@ -75,6 +85,7 @@ program qra65sim
   write(*,1000) 
 1000 format('File    TR   Freq Mode  S/N   DT    Dop  Message'/60('-'))
 
+  nsync=0
   do ifile=1,nfiles                  !Loop over requested number of files
      if(ntrperiod.lt.60) then
         write(fname,1002) ifile         !Output filename
@@ -169,10 +180,23 @@ program qra65sim
      write(10) h,iwave(1:npts)                !Save the .wav file
      close(10)
 
-!     do i=1,NMAX
+     if(lsync) then
-!        write(15,3020) i/12000.0,iwave(i)
+        cd=' '
-!3020    format(f10.6,i8)
+        if(ifile.eq.nfiles) cd='d'
-!     enddo
+        nfqso=nint(f0)
+        ntol=100
+        call sync_qra65(iwave,npts,mode65,nsps,nfqso,ntol,xdt2,f02,snr2)
+        terr=1.01/(8.0*baud)
+        ferr=1.01*mode65*baud
+        if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1
+        open(40,file='sync65.out',status='unknown',position='append')
+        write(40,1030) ifile,65,csubmode,snrdb,fspread,xdt2-xdt,f02-f0,   &
+             snr2,nsync,cd
+1030    format(i4,i3,1x,a1,2f7.1,f7.2,2f8.1,i5,1x,a1)
+        close(40)
+     endif
   enddo
+  if(lsync) write(*,1040) snrdb,nfiles,nsync
+1040 format('SNR:',f6.1,'   nfiles:',i5,'   nsynced:',i5)
 
 999 end program qra65sim

lib/qra65_decode.f90

@@ -93,8 +93,9 @@ contains
 
 ! AP control could be done differently, but this works well:
     maxaptype=0
-    if(ndepth.eq.2) maxaptype=3
+!    if(ndepth.eq.2) maxaptype=3
-    if(ndepth.eq.3) maxaptype=5
+!    if(ndepth.eq.3) maxaptype=5
+    if(ndepth.ge.2) maxaptype=5       !###
 
     if(nc1.ne.nc1z .or. nc2.ne.nc2z .or. ng2.ne.ng2z .or.            &
          maxaptype.ne.maxaptypez) then

lib/qra_loops.f90

@@ -69,7 +69,6 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
            endif
         enddo  ! ibw (b90 loop)
         if(iand(ndepth,3).lt.3 .and. irc.ge.0) go to 100
-        if(irc.eq.0) go to 100
      enddo  ! idt (DT loop)
   enddo  ! idf (f0 loop)
 
@@ -85,5 +84,9 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
      ibw=ibwkeep
   endif
 
-200 return
+200 continue
+  write(53,3053) idt,idf,ibw,irc,b90,xdt,f0,snr2
+3053 format(4i5,f7.1,f7.2,2f7.1)
+  
+  return
 end subroutine qra_loops

lib/sync_qra65.f90

@@ -1,6 +1,8 @@
 subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
 
-! Look for the sync vector in a QRA65 signal.
+! Detect and align with the QRA65 sync vector, returning time and frequency
+! offsets and SNR estimate.
+
 ! Input:  iwave(0:nmax-1)        Raw data
 !         mode65                 Tone spacing 1 2 4 8 16 (A-E)
 !         nsps                   Samples per symbol at 12000 Sa/s
@@ -10,7 +12,7 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
 !         f0                     Frequency of sync tone
 !         snr1                   Relative SNR of sync signal
   
-  parameter (NSTEP=4)                    !Quarter-symbol steps
+  parameter (NSTEP=8)                    !Step size nsps/NSTEP
   integer*2 iwave(0:nmax-1)              !Raw data
   integer isync(22)                      !Indices of sync symbols
   integer ijpk(2)                        !Indices i and j at peak of ccf
@@ -28,7 +30,7 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
   iz=5000.0/df                           !Uppermost frequency bin, at 5000 Hz
   txt=85.0*nsps/12000.0
   jz=(txt+1.0)*12000.0/istep             !Number of quarter-symbol steps
-  if(nsps.ge.7680) jz=(txt+2.0)*12000.0/istep   !For TR 60 s and higher
+  if(nsps.ge.6912) jz=(txt+2.0)*12000.0/istep   !For TR 60 s and higher
 
   allocate(s1(iz,jz))
   allocate(c0(0:nfft-1))
@@ -41,11 +43,11 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
   endif
 
   fac=1/32767.0
-  do j=1,jz                     !Compute symbol spectra at quarter-symbol steps
+  do j=1,jz                              !Compute symbol spectra at step size
      ia=(j-1)*istep
      ib=ia+nsps-1
      k=-1
-     do i=ia,ib,2
+     do i=ia,ib,2          !Load iwave data into complex array c0, for r2c FFT
         xx=iwave(i)
         yy=iwave(i+1)
         k=k+1
@@ -57,11 +59,12 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
         s1(i,j)=real(c0(i))**2 + aimag(c0(i))**2
      enddo
 ! For large Doppler spreads, should we smooth the spectra here?
+     call smo121(s1(1:iz,j),iz)
   enddo
 
   i0=nint(nfqso/df)                           !Target QSO frequency
   call pctile(s1(i0-64:i0+192,1:jz),129*jz,40,base)
-  s1=s1/base                                  !Maybe should subtract 1.0 here?
+  s1=s1/base - 1.0
 
 ! Apply fast AGC
   s1max=20.0                                  !Empirical choice
@@ -70,32 +73,30 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
      if(smax.gt.s1max) s1(i0-64:i0+192,j)=s1(i0-64:i0+192,j)*s1max/smax
   enddo
 
-  dt4=nsps/(NSTEP*12000.0)                    !1/4 of symbol duration
+  dtstep=nsps/(NSTEP*12000.0)                 !Step size in seconds
-  j0=0.5/dt4
+  j0=0.5/dtstep
-  if(nsps.ge.7680) j0=1.0/dt4                 !Nominal index for start of signal
+  if(nsps.ge.6192) j0=1.0/dtstep              !Nominal index for start of signal
   
   ccf=0.
   ia=min(64,nint(ntol/df))
-  lag1=-1.0/dt4
+  lag1=-1.0/dtstep
-  lag2=4.0/dt4 + 0.9999
+  lag2=4.0/dtstep + 0.9999
 
-  do i=-ia,ia
   do lag=lag1,lag2
      do k=1,85
         n=NSTEP*(k-1) + 1
         j=n+lag+j0
         if(j.ge.1 .and. j.le.jz) then
-              ccf(i,lag)=ccf(i,lag) + sync(k)*s1(i0+i,j)
+           ccf(-ia:ia,lag)=ccf(-ia:ia,lag) + sync(k)*s1(i0-ia:i0+ia,j)
         endif
      enddo
   enddo
-  enddo
 
   ijpk=maxloc(ccf)
   ipk=ijpk(1)-65
   jpk=ijpk(2)-27
   f0=nfqso + ipk*df
-  xdt=jpk*dt4
+  xdt=jpk*dtstep
   snr1=maxval(ccf)/22.0
 
   return