QRA66 now works with qra66sim and jt9 and TR periods 15 30 60 120 300 s.

lib/jt9.f90

@@ -249,7 +249,7 @@ program jt9
               call timer('symspec ',1)
            endif
            nhsym0=nhsym
-           if(nhsym.ge.181 .and. mode.ne.240 .and. mode.ne.241) exit
+           if(nhsym.ge.181 .and. mode.ne.240 .and. mode.ne.241 .and. mode.ne.66) exit
         endif
      enddo
      close(unit=wav%lun)

lib/qra/qra66/qra66sim.f90

@@ -74,7 +74,7 @@ program qra66sim
   h=default_header(12000,npts)
 
   write(*,1000) 
-1000 format('File     Freq  A-E   S/N   DT    Dop  Message'/60('-'))
+1000 format('File    TR   Freq Mode  S/N   DT    Dop  Message'/60('-'))
 
   do ifile=1,nfiles                  !Loop over requested number of files
      write(fname,1002) ifile         !Output filename
@@ -91,11 +91,12 @@ program qra66sim
      bandwidth_ratio=2500.0/6000.0
      sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
      if(snrdb.gt.90.0) sig=1.0
-     write(*,1020) ifile,f0,csubmode,xsnr,xdt,fspread,msgsent
+     write(*,1020) ifile,ntrperiod,f0,csubmode,xsnr,xdt,fspread,msgsent
-1020    format(i4,f10.3,2x,a1,2x,f5.1,f6.2,f6.1,1x,a22)
+1020    format(i4,i6,f7.1,2x,a1,2x,f5.1,f6.2,f6.1,1x,a22)
      phi=0.d0
      dphi=0.d0
-     k=(xdt+0.5)*12000                   !Start audio at t = xdt + 0.5 s
+     k=(xdt+0.5)*12000                   !Start audio at t=xdt+0.5 s (TR=15 and 30 s)
+     if(ntrperiod.ge.60) k=(xdt+1.0)*12000   !TR 60+ at t = xdt + 1.0 s
      isym0=-99
      do i=1,npts                         !Add this signal into cdat()
         isym=i/nsps + 1
@@ -137,16 +138,6 @@ program qra66sim
            cspread(nfft-i)=z
         enddo
 
-!        do i=0,nfft-1
-!           f=i*df
-!           if(i.gt.nh) f=(i-nfft)*df
-!           s=real(cspread(i))**2 + aimag(cspread(i))**2
-!          write(13,3000) i,f,s,cspread(i)
-!3000      format(i5,f10.3,3f12.6)
-!        enddo
-!        s=real(cspread(0))**2 + aimag(cspread(0))**2
-!        write(13,3000) 1024,0.0,s,cspread(0)
-
         call four2a(cspread,nfft,1,1,1)             !Transform to time domain
 
         sum=0.

lib/qra66_decode.f90

@@ -34,7 +34,7 @@ contains
     use timer_module, only: timer
     use packjt
     use, intrinsic :: iso_c_binding
-    parameter (NMAX=60*12000)             !### Needs to be 300*12000 ###
+    parameter (NMAX=300*12000)             !### Needs to be 300*12000 ###
     class(qra66_decoder), intent(inout) :: this
     procedure(qra66_decode_callback) :: callback
     character(len=12) :: mycall, hiscall
@@ -43,7 +43,7 @@ contains
     integer*2 iwave(NMAX)                 !Raw data
     integer dat4(12)
     logical lapdx,ltext
-    complex c0(0:NMAX-1)                  !Analytic signal, 6000 S/s
+    complex, allocatable :: c0(:)         !Analytic signal, 6000 S/s
     real s3(-64:127,63)
     real s3a(-64:127,63)
     real a(5)
@@ -52,6 +52,8 @@ contains
 
     nfft1=ntrperiod*12000
     nfft2=ntrperiod*6000
+    allocate (c0(0:nfft1-1))
+    
     if(ntrperiod.eq.15) then
        nsps=1800
     else if(ntrperiod.eq.30) then
@@ -67,11 +69,12 @@ contains
     endif
     baud=12000.0/nsps
     df1=12000.0/nfft1
-    print*,'aaa',ntrperiod,hisgrid,nsps,baud
+
-    do i=1,NMAX
+!    do i=1,NMAX
-       write(61,3061) i/12000.0,iwave(i)/32767.0
+!       write(61,3061) i/12000.0,iwave(i)/32767.0
-3061   format(2f12.6)
+!3061   format(2f12.6)
-    enddo
+!    enddo
+
     this%callback => callback
     
     if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso  !Silence warning
@@ -107,19 +110,20 @@ contains
 
 ! Downsample to give complex data at 6000 S/s
     fac=2.0/nfft1
-    c0=fac*iwave
+    c0=fac*iwave(1:nfft1)
     call four2a(c0,nfft1,1,-1,1)           !Forward c2c FFT
     c0(nfft2/2+1:nfft2)=0.                 !Zero the top half
     c0(0)=0.5*c0(0)
     call four2a(c0,nfft2,1,1,1)            !Inverse c2c FFT
     call timer('sync66  ',0)
-    call sync66(iwave,60*12000,nsps,nfqso,ntol,xdt,f0,snr1)    !### 300*12000 ###
+    call sync66(iwave,ntrperiod*12000,nsps,nfqso,ntol,xdt,f0,snr1)
     call timer('sync66  ',1)
-    jpk=(xdt+0.5)*6000 - 384               !### Empirical ###
+    jpk=(xdt+0.5)*6000 - 384                       !### Empirical ###
+    if(ntrperiod.ge.60) jpk=(xdt+1.0)*6000 - 384   !### TBD ###
     if(jpk.lt.0) jpk=0
     a=0.
     a(1)=-(f0 + 2.0*baud)                  !Data tones start 2 bins higher
-    call twkfreq(c0,c0,15*6000,6000.0,a)
+    call twkfreq(c0,c0,ntrperiod*6000,6000.0,a)
     xdt=jpk/6000.0 - 0.5
     call spec66(c0(jpk:),nsps/2,s3)
 

lib/sync66.f90

@@ -5,7 +5,7 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1)
   integer b11(11)                          !Barker 11 code
   integer ijpk(2)                          !Indices i and j at peak of sync_sig
   real, allocatable :: s1(:,:)             !Symbol spectra
-  real, allocatable :: x(:)                !Work array; 2FSK sync modulation
+  real, allocatable :: x(:)                !Work array; demoduated 2FSK sync signal
   real sync(4*85)                          !sync vector
   real sync_sig(-64:64,-15:15)
   complex, allocatable :: c0(:)            !Complex spectrum of symbol
@@ -13,11 +13,13 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1)
   data sync(1)/99.0/
   save sync
 
-  nfft=2*NSPS
+  nfft=2*nsps
   df=12000.0/nfft
   istep=nsps/NSTEP
-  iz=5000.0/df
+  iz=5000.0/df                             !Uppermost frequency bin, at 5000 Hz
-  jz=352
+  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
   allocate(s1(iz,jz))
   allocate(x(jz))
   allocate(c0(0:nsps))
@@ -27,12 +29,13 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1)
      do k=1,22
         kk=k
         if(kk.gt.11) kk=k-11
-        sync(16*k-15)=2.0*b11(kk) - 1.0
+        i=4*NSTEP*(k-1) + 1
+        sync(i)=2.0*b11(kk) - 1.0
      enddo
   endif
 
   fac=1/32767.0
-  do j=1,JZ                                  !Compute symbol spectra
+  do j=1,jz                          !Compute symbol spectra at quarter-symbol steps
      ia=(j-1)*istep
      ib=ia+nsps-1
      k=-1
@@ -44,34 +47,43 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1)
      enddo
      c0(k+1:nfft/2)=0.
      call four2a(c0,nfft,1,-1,0)              !r2c FFT
-     do i=1,IZ
+     do i=1,iz
         s1(i,j)=real(c0(i))**2 + aimag(c0(i))**2
      enddo
   enddo
 
   i0=nint(nfqso/df)
-  call pctile(s1(i0-64:i0+192,1:JZ),129*JZ,40,base)
+  call pctile(s1(i0-64:i0+192,1:jz),129*jz,40,base)
   s1=s1/base
   s1max=20.0
 
 ! Apply AGC
-  do j=1,JZ
+  do j=1,jz
      x(j)=maxval(s1(i0-64:i0+192,j))
      if(x(j).gt.s1max) s1(i0-64:i0+192,j)=s1(i0-64:i0+192,j)*s1max/x(j)
   enddo
 
-  dt4=nsps/(NSTEP*12000.0)
-  j0=0.5/dt4
-
   sync_sig=0.
   ia=min(64,nint(ntol/df))
+  dt4=nsps/(NSTEP*12000.0)                      !duration of 1/4 symbol
+  lag2=nint(0.5/dt4)
+  if(nsps.ge.7680) lag2=nint(1.0/dt4)
+  lag1=-lag2
+
+  jadd=11
+  if(nsps.ge.3600) jadd=7
+  if(nsps.ge.7680) jadd=6
+  if(nsps.ge.16000) jadd=3
+  if(nsps.ge.41472) jadd=1
+  
   do i=-ia,ia
-     x=s1(i0+2+i,:)-s1(i0+i,:)
+     x=s1(i0+2+i,:)-s1(i0+i,:)                  !Do the 2FSK demodulation
+!     do lag=lag1,lag2
      do lag=-15,15
-! Make this simpler: just add the 22 nonzero values?
+        do k=1,22
-        do n=1,4*85
+           n=4*NSTEP*(k-1) + 1
-           j=n+lag+11
+           j=n+lag+jadd
-           if(j.ge.1 .and. j.le.JZ) sync_sig(i,lag)=sync_sig(i,lag) + sync(n)*x(j)
+           if(j.ge.1 .and. j.le.jz) sync_sig(i,lag)=sync_sig(i,lag) + sync(n)*x(j)
         enddo
      enddo
   enddo
@@ -83,6 +95,10 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1)
 ! Use peakup() here?
   f0=nfqso + ii*df
   jdt=jj
+  
+!  j0=0.5/dt4
+!  if(nsps.ge.7680) j0=1.0/dt4
+
   tsym=nsps/12000.0
   xdt=jdt*tsym/4.0
 
@@ -100,4 +116,3 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1)
 
   return
 end subroutine sync66
-