QRA65 now decodes using qra_loops() -- the same inner loops as QRA64. Very effective!

CMakeLists.txt

@@ -570,7 +570,6 @@ set (wsjt_FSRCS
   lib/softsym9w.f90
   lib/shell.f90
   lib/spec64.f90
-  lib/spec_qra65.f90
   lib/spec9f.f90
   lib/stdmsg.f90
   lib/subtract65.f90
@@ -1317,7 +1316,7 @@ if (${OPENMP_FOUND} OR APPLE)
   endif (APPLE)
   if (WIN32)
     set_target_properties (jt9 PROPERTIES
-      LINK_FLAGS -Wl,--stack,16777216
+      LINK_FLAGS -Wl,--stack,28000000
       )
   endif ()
   target_link_libraries (jt9 wsjt_fort_omp wsjt_cxx fort_qt)
@@ -1454,7 +1453,7 @@ else ()
     )
   if (WIN32)
     set_target_properties (wsjtx PROPERTIES
-      LINK_FLAGS -Wl,--stack,16777216
+      LINK_FLAGS -Wl,--stack,28000000
       )
   endif ()
 endif ()

lib/ana64.f90

@@ -3,13 +3,12 @@ subroutine ana64(dd,npts,c0)
   use timer_module, only: timer
 
   parameter (NMAX=60*12000)                  !Max size of raw data at 12000 Hz
-  parameter (NSPS=3456)                      !Samples per symbol at 6000 Hz
-  parameter (NSPC=7*NSPS)                    !Samples per Costas array
   real dd(NMAX)                              !Raw data
   complex c0(0:720000)                       !Complex spectrum of dd()
   save
 
-  nfft1=672000
+!  nfft1=672000
+  nfft1=720000
   nfft2=nfft1/2
   df1=12000.0/nfft1
   fac=2.0/nfft1

lib/qra64a.f90

@@ -81,8 +81,6 @@ subroutine qra64a(dd,npts,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
      snr2=0.
   endif
   nfreq=nint(f0)
-  write(71,3071) idf,idt,ncall,irc,nsnr,dtx,nfreq,decoded
-3071 format(5i5,f7.2,i6,2x,a22)
 
 900 if(irc.lt.0) then
      sy=max(1.0,sync)

lib/qra65_decode.f90

@@ -41,18 +41,20 @@ contains
     character(len=6) :: hisgrid
     character*37 decoded
     integer*2 iwave(NMAX)                 !Raw data
+    real dd(NMAX)                         !Raw data
     integer dat4(12)
     logical lapdx,ltext
+    complex, allocatable :: c00(:)        !Analytic signal, 6000 S/s
     complex, allocatable :: c0(:)         !Analytic signal, 6000 S/s
     real, allocatable, save :: s3(:,:)    !Symbol spectra
     real, allocatable, save :: s3a(:,:)   !Symbol spectra for avg messages
-    real a(5)
     data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/,nsubmodez/-1/
     save nc1z,nc2z,ng2z,maxaptypez,nsave,nsubmodez
 
     mode65=2**nsubmode
     nfft1=ntrperiod*12000
     nfft2=ntrperiod*6000
+    allocate (c00(0:nfft1-1))
     allocate (c0(0:nfft1-1))
 
     if(nsubmode.ne.nsubmodez) then
@@ -75,16 +77,10 @@ contains
     else
       stop 'Invalid TR period'
    endif
+   npts=ntrperiod*12000
    baud=12000.0/nsps
    df1=12000.0/nfft1
-
-!    do i=1,12000*ntrperiod
-!       write(61,3061) i/12000.0,iwave(i)/32767.0
-!3061   format(2f12.6)
-!    enddo
-
    this%callback => callback    
-    
    if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso  !Silence warning
 
 ! Prime the QRA decoder for possible use of AP
@@ -119,61 +115,22 @@ contains
    call sync_qra65(iwave,ntrperiod*12000,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
    call timer('sync_q65',1)
 
-! Downsample to give complex data at 6000 S/s
+   jpk0=(xdt+1.0)*6000   !###
-    call timer('down_q65',0)
+   if(jpk0.lt.0) jpk0=0
-    fac=2.0/nfft1
-    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
-    a=0.
-    a(1)=-(f0 + mode65*baud)             !Data tones start mode65 bins higher
-    call twkfreq(c0,c0,ntrperiod*6000,6000.0,a)
-    call timer('down_q65',1)
 
-    jpk=(xdt+0.5)*6000 - 384                       !### Empirical ###
+   fac=1.0/32767.0
-    if(ntrperiod.ge.60) jpk=(xdt+1.0)*6000 - 384   !### TBD ??? ###
+   dd=fac*iwave
-    if(jpk.lt.0) jpk=0
+!   npts=648000
-    xdt=jpk/6000.0 - 0.5
+   minsync=-2
-    LL=64*(mode65+2)
+   nmode=65
-    NN=63
-    call timer('spec_q65',0)
-    call spec_qra65(c0(jpk:),nsps/2,s3,LL,NN)  !Compute synced symbol spectra
-    call timer('spec_q65',1)
 
-    do j=1,63                              !Normalize to symbol baseline
+   call ana64(dd,npts,c00)
-       call pctile(s3(:,j),LL,40,base)
-       s3(:,j)=s3(:,j)/base
-    enddo
 
-    LL2=64*(mode65+1)-1
+   call timer('qraloops',0)
-    s3max=20.0
+   call qra_loops(c00,npts/2,nmode,mode65,nsubmode,nFadingModel,minsync,     &
-    do j=1,63                              !Apply AGC to suppress pings
+        ndepth,nc1,nc2,ng2,naptype,jpk0,xdt,f0,width,snr2,s3,irc,dat4)
-     xx=maxval(s3(-64:LL2,j))
+   if(nmode.eq.65) xdt=xdt+0.4        !### Empirical  WHY ??? ###
-     if(xx.gt.s3max) s3(-64:LL2,j)=s3(-64:LL2,j)*s3max/xx
+   call timer('qraloops',1)
-    enddo
-
-! Call Nico's QRA64 decoder
-    call timer('qra64_de',0)
-    call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90,      &
-         nFadingModel,dat4,snr2,irc)
-    call timer('qra64_de',1)
-
-    if(irc.lt.0) then
-! No luck so far. Try for an average decode.
-       call timer('qra64_av',0)
-       s3a=s3a+s3
-       nsave=nsave+1
-       if(nsave.ge.2) then
-          call qra64_dec(s3a,nc1,nc2,ng2,naptype,0,nSubmode,b90,      &
-               nFadingModel,dat4,snr2,irc)
-          if(irc.ge.0) irc=100*nsave + irc
-       endif
-       call timer('qra64_av',1)
-    endif
-    snr2=snr2 + db(6912.0/nsps)
-    if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2)
 
     decoded='                                     '
     if(irc.ge.0) then

lib/qra_loops.f90

@@ -1,5 +1,5 @@
 subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
-     ndepth,nc1,nc2,ng2,naptype,jpk0,dtx,f0,width,snr2,s3,irc,dat4)
+     ndepth,nc1,nc2,ng2,naptype,jpk0,xdt,f0,width,snr2,s3,irc,dat4)
 
   use timer_module, only: timer
   parameter (LN=1152*63)
@@ -10,7 +10,6 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
   integer dat4(12),dat4x(12)       !Decoded message (as 12 integers)
   integer nap(0:11)                !AP return codes
   data nap/0,2,3,2,3,4,2,3,6,4,6,6/
-!  save
 
   irc=-99
   s3lim=20.
@@ -26,6 +25,8 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
   NN=63
   napmin=99
   ncall=0
+  nsps=3456                                   !QRA64
+  if(mode.eq.65) nsps=3840                    !QRA65
 
   do idf0=1,11
      idf=idf0/2
@@ -37,24 +38,8 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
         idt=idt0/2
         if(mod(idt0,2).eq.0) idt=-idt
         jpk=jpk0 + 750*idt
-        if(mode.eq.64) then
-           call spec64(c0,jpk,s3a,LL,NN)
-        else
         if(jpk.lt.0) jpk=0
-           call timer('spec_q65',0)
+        call spec64(c0,nsps,mode,jpk,s3a,LL,NN)
-           call spec_qra65(c0(jpk:),nsps2,s3,LL,NN) !Get synced symbol spectra
-           call timer('spec_q65',1)
-!           do j=1,63                              !Normalize to symbol baseline
-!              call pctile(s3(:,j),LL,40,base)
-!              s3(:,j)=s3(:,j)/base
-!           enddo
-!           LL2=64*(mode65+1)-1
-!           s3max=20.0
-!           do j=1,63                              !Apply AGC to suppress pings
-!              xx=maxval(s3(-64:LL2,j))
-!              if(xx.gt.s3max) s3(-64:LL2,j)=s3(-64:LL2,j)*s3max/xx
-!           enddo
-        endif
         call pctile(s3a,LL*NN,40,base)
         s3a=s3a/base
         where(s3a(1:LL*NN)>s3lim) s3a(1:LL*NN)=s3lim
@@ -77,7 +62,7 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
               snr2x=snr2
               napmin=nap(iirc)
               irckeep=irc
-              dtxkeep=jpk/6000.0 - 1.0
+              xdtkeep=jpk/6000.0 - 1.0
               f0keep=-a(1)
               idfkeep=idf
               idtkeep=idt
@@ -94,7 +79,7 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync,  &
      b90=b90x
      snr2=snr2x
      irc=irckeep
-     dtx=dtxkeep
+     xdt=xdtkeep
      f0=f0keep
      idt=idtkeep
      idf=idfkeep

lib/spec64.f90

@@ -1,13 +1,17 @@
-subroutine spec64(c0,jpk,s3,LL,NN)
+subroutine spec64(c0,nsps,mode,jpk,s3,LL,NN)
 
-  parameter (NSPS=3456)                      !Samples per symbol at 6000 Hz
+  parameter (MAXFFT=3840)
   complex c0(0:360000)                       !Complex spectrum of dd()
-  complex cs(0:NSPS-1)                       !Complex symbol spectrum
+  complex cs(0:MAXFFT-1)                     !Complex symbol spectrum
   real s3(LL,NN)                             !Synchronized symbol spectra
   real xbase0(LL),xbase(LL)
+  integer isync(22)                          !Indices of sync symbols
+  data isync/1,9,12,13,15,22,23,26,27,33,35,38,46,50,55,60,62,66,69,74,76,85/
 
   nfft=nsps
   fac=1.0/nfft
+
+  if(mode.eq.64) then
      do j=1,NN
         jj=j+7                                  !Skip first Costas array
         if(j.ge.33) jj=j+14                     !Skip middle Costas array
@@ -21,6 +25,26 @@ subroutine spec64(c0,jpk,s3,LL,NN)
            s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2
         enddo
      enddo
+  else
+     j=0
+     n=1
+     do k=1,84
+        if(k.eq.isync(n)) then
+           n=n+1
+           cycle
+        endif
+        j=j+1
+        ja=(k-1)*nsps + jpk
+        jb=ja+nsps-1
+        cs(0:nfft-1)=fac*c0(ja:jb)
+        call four2a(cs,nsps,1,-1,1)             !c2c FFT to frequency
+        do ii=1,LL
+           i=ii-65
+           if(i.lt.0) i=i+nsps
+           s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2
+        enddo
+     enddo
+  endif
 
   df=6000.0/nfft
   do i=1,LL
@@ -38,5 +62,18 @@ subroutine spec64(c0,jpk,s3,LL,NN)
      s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization
   enddo
 
+!  print*,'a',LL,NN,jpk
+!  df=6000.0/nfft
+!  do i=1,LL
+!     write(71,3071) i,i-65,i*df,(s3(i,j),j=1,4)
+!3071 format(2i8,f10.3,4e12.3)
+!  enddo
+!
+!  do j=1,NN
+!     write(72,3072) j,maxloc(s3(1:LL,j)),maxloc(s3(1:LL,j))-65
+!3072 format(3i8)
+!  enddo
+!  if(nfft.ne.-999) stop
+
   return
 end subroutine spec64