OPtimize the

lib/q65_decode.f90

@@ -61,8 +61,6 @@ contains
     logical lapcqonly,unpk77_success
     complex, allocatable :: c00(:)        !Analytic signal, 6000 Sa/s
     complex, allocatable :: c0(:)         !Analytic signal, 6000 Sa/s
-    data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/,nsubmodez/-1/
-    save nc1z,nc2z,ng2z,maxaptypez,nsubmodez
 
     mode65=2**nsubmode
     nfft1=ntrperiod*12000
@@ -90,23 +88,16 @@ contains
     this%callback => callback
     if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso  !Silence warning
     nFadingModel=1
-
-! AP control could be done differently, but this works well:
-    maxaptype=0
-!    if(ndepth.eq.2) maxaptype=3
-!    if(ndepth.eq.3) maxaptype=5
-    if(ndepth.ge.2) maxaptype=5       !###
-    minsync=-2
     call qra_params(ndepth,maxaptype,idfmax,idtmax,ibwmin,ibwmax,maxdist)
-    naptype=maxaptype
 
     call timer('sync_q65',0)
-    call sync_q65(iwave,ntrperiod*12000,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
+    call sync_q65(iwave,ntrperiod*12000,mode65,nsps,nfqso,ntol,xdt,f0,   &
+         snr1,width)
     call timer('sync_q65',1)
 
     irc=-1
     if(snr1.lt.2.5) go to 100
-    jpk0=(xdt+1.0)*6000                      !###
+    jpk0=(xdt+1.0)*6000                      !### Is this OK?
     if(ntrperiod.le.30) jpk0=(xdt+0.5)*6000  !###
     if(jpk0.lt.0) jpk0=0
     fac=1.0/32767.0
@@ -122,7 +113,6 @@ contains
     if(nQSOprogress.eq.5) npasses=3
     if(lapcqonly) npasses=1
     do ipass=0,npasses
-!       print*,'A',nQSOprogress,ipass,npasses
        apmask=0
        apsymbols=0
        if(ipass.ge.1) then
@@ -135,9 +125,6 @@ contains
           write(c78,1050) apsymbols1
           read(c78,1060) apsymbols
           apsymbols(13)=apsymbols(13)/2              !Fixup for c77-->c78
-!          write(72,3060) 'A',ipass,apmask,apmask
-!3060      format(a1,i1,1x,13b6.6/3x,13i6)
-!          write(72,3060) 'B',ipass,apsymbols,apsymbols
        endif
        call timer('q65loops',0)
        call q65_loops(c00,npts/2,nsps/2,nmode,mode65,nsubmode,nFadingModel,  &
@@ -149,24 +136,19 @@ contains
 
 100 decoded='                                     '
     if(irc.ge.0) then
-!###       
+!###
-!       irc=(irc/100) * 100                  !### TEMPORARY ??? ###
        navg=irc/100
        irc=ipass
 !###
        write(c77,1000) dat4
 1000   format(12b6.6,b5.5)
 
-!       write(72,3080) 'C',ipass,c77,'0'
-!3080   format(a1,i1,1x,a77,a1)
-!       write(72,3060) 'C',ipass,dat4,dat4
-       
        call unpack77(c77,0,decoded,unpk77_success) !Unpack to get msgsent
        nsnr=nint(snr2)
        call this%callback(nutc,sync,nsnr,xdt,f0,decoded,              &
             irc,qual,ntrperiod,fmid,w50)
     else
-       ! Report sync, even if no decode.
+! Report sync, even if no decode.
        nsnr=db(snr1) - 35.0
        call this%callback(nutc,sync,nsnr,xdt,f0,decoded,              &
             irc,qual,ntrperiod,fmid,w50)

lib/qra/q65/q65_loops.f90

@@ -21,6 +21,7 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
   data nap/0,2,3,2,3,4,2,3,6,4,6,6/,nsave/0/
   save nsave,s3avg
 
+  ircbest=9999
   allocate(c0(0:npts2-1))
   irc=-99
   s3lim=20.
@@ -32,7 +33,7 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
   LL=64*(mode64+2)
   NN=63
   napmin=99
-  ncall=0
+  baud=6000.0/nsps
 
   do iavg=0,1
      if(iavg.eq.1) then
@@ -43,13 +44,13 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
         ndf=idf/2
         if(mod(idf,2).eq.0) ndf=-ndf
         a=0.
-        a(1)=-(f0+0.4*ndf)
+        a(1)=-(f0+0.5*baud*ndf)
         call twkfreq(c00,c0,npts2,6000.0,a)
         do idt=1,idtmax
            ndt=idt/2
            if(iavg.eq.0) then
               if(mod(idt,2).eq.0) ndt=-ndt
-              jpk=jpk0 + 240*ndt                  !240/6000 = 0.04 s = tsym/32
+              jpk=jpk0 + nsps*ndt/16              !tsym/16
               if(jpk.lt.0) jpk=0
               call timer('spec64  ',0)
               call spec64(c0,nsps,mode,mode64,jpk,s3,LL,NN)
@@ -60,13 +61,14 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
            else
               s3(1:LL*NN)=s3avg(1:LL*NN)
            endif
-           do ibw=ibwmax,ibwmin,-2
+           do ibw=ibwmin,ibwmax
-              ndist=ndf**2 + ndt**2 + ((ibwmax-ibw)/2)**2
+              nbw=ibw
+              ndist=ndf**2 + ndt**2 + ((nbw-2))**2
               if(ndist.gt.maxdist) cycle
-              b90=1.728**ibw
+!              b90=1.728**ibw
+              b90=3.0**nbw
               if(b90.gt.230.0) cycle
 !              if(b90.lt.0.15*width) exit
-              ncall=ncall+1
               call timer('q65_intr',0)
               call q65_intrinsics_ff(s3,nsubmode,b90,nFadingModel,s3prob)
               call timer('q65_intr',1)
@@ -74,11 +76,11 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
               call timer('q65_dec ',0)
               call q65_dec(s3,s3prob,APmask,APsymbols,esnodb,dat4,irc)
               call timer('q65_dec ',1)
+              if(irc.ge.0) go to 100
               ! irc > 0 ==> number of iterations required to decode
               !  -1 = invalid params
               !  -2 = decode failed
               !  -3 = CRC mismatch
-              if(irc.ge.0) go to 100
            enddo  ! ibw (b90 loop)
         enddo  ! idt (DT loop)
      enddo  ! idf (f0 loop)
@@ -86,7 +88,7 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
         a=0.
         a(1)=-f0
         call twkfreq(c00,c0,npts2,6000.0,a)
-        jpk=3000                       !###  These definitions need work ###
+        jpk=3000                       !###  Are these definitions OK?
         if(nsps.ge.3600) jpk=6000      !###  TR >= 60 s
         call spec64(c0,nsps,mode,mode64,jpk,s3,LL,NN)
         call pctile(s3,LL*NN,40,base)
@@ -98,26 +100,24 @@ subroutine q65_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel,   &
      if(iavg.eq.0 .and. nsave.lt.2) exit
   enddo  ! iavg
 
-100 if(mode.eq.65 .and. nsps.eq.7200/2) xdt=xdt+0.4 !### Empirical -- WHY ??? ###
+100  if(irc.ge.0) then
-
-  if(irc.ge.0) then
      navg=nsave
-     baud=6000.0/nsps
      snr2=esnodb - db(2500.0/baud)
      if(iavg.eq.0) navg=0
-     !### For tests only:
+!### For tests only:
      open(53,file='fort.53',status='unknown',position='append')
      write(c77,1100) dat4
 1100 format(12b6.6,b5.5)
      call unpack77(c77,0,decoded,unpk77_success) !Unpack to get msgsent
-     write(53,3053) idf,idt,ibw,b90,xdt,f0,snr2,ndist,irc,ipass,navg,  &
+     write(53,3053) ndf,ndt,nbw,b90,xdt,f0,snr2,ndist,irc,ipass,navg,  &
           trim(decoded)
 3053 format(3i5,f7.1,f7.2,2f7.1,4i4,2x,a)
      close(53)
-     !###  
+!###  
      nsave=0
      s3avg=0.
      irc=irc + 100*navg
   endif
+
   return
 end subroutine q65_loops

lib/qra64a.f90

@@ -22,7 +22,7 @@ subroutine qra64a(dd,npts,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
   nft=99
   if(nfqso.lt.nf1 .or. nfqso.gt.nf2) go to 900
 
-  mycall=mycall_12(1:6)                     !### May need fixing ###
+  mycall=mycall_12(1:6)                     !### May need fixing? ###
   hiscall=hiscall_12(1:6)
   hisgrid=hisgrid_6(1:4)
   call packcall(mycall,nc1,ltext)

lib/sync_q65.f90

@@ -1,4 +1,4 @@
-subroutine sync_q65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
+subroutine sync_q65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1,width)
 
 ! Detect and align with the Q65 sync vector, returning time and frequency
 ! offsets and SNR estimate.
@@ -118,9 +118,17 @@ subroutine sync_q65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1)
 !  enddo
 
 !  do i=-ia,ia
-!     write(56,3056) i*df,ccf(i,0)/rms
+!     write(56,3056) i*df,ccf(i,jpk)/rms
 !3056 format(2f10.3)
 !  enddo
+!  flush(56)
+
+  acf0=dot_product(ccf(-ia:ia,jpk),ccf(-ia:ia,jpk))
+  do i=1,ia
+     acf=dot_product(ccf(-ia:ia,jpk),ccf(-ia+i:ia+i,jpk))
+     if(acf.le.0.5*acf0) exit
+  enddo
+  width=i*1.414*df
 
   return
 end subroutine sync_q65