Many adjustments to QRA64 decoder. Best performance yet on the 114 benchmark files.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7381 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

CMakeLists.txt

@@ -300,6 +300,7 @@ endif (WIN32)
 set (wsjt_FSRCS
   lib/afc65b.f90
   lib/afc9.f90
+  lib/ana64.f90
   lib/ana932.f90
   lib/analytic.f90
   lib/astro.f90

lib/Makefile

@@ -17,13 +17,14 @@ CFLAGS = -O2 -I.
 %.o: %.F90
 	${FC} ${FFLAGS} -c $<
 
-all:	 qratest.exe
+all:	 qratest.exe tplt.exe
 
-OBJS1 = qratest.o qra64a.o sync64.o four2a.o smo.o smo121.o averms.o \
+OBJS1 = qratest.o qra64a.o ana64.o sync64a.o four2a.o smo.o smo121.o averms.o \
 	timer_module.o packjt.o twkfreq.o spec64.o fmtmsg.o pctile.o \
 	grid2deg.o deg2grid.o shell.o badmsg.o qra64_subs.o \
 	qracodes.o npfwht.o pdmath.o qra12_63_64_irr_b.o \
-	qra13_64_64_irr_e.o qra64.o
+	qra13_64_64_irr_e.o qra64.o image.o \
+	zplt64a.o zplt64b.o lorentzian.o fchisq0.o peakup.o sync64.o
 
 qra64_subs.o: qra/qra64/qra64_subs.c
 	gcc -c -O2 -o qra64_subs.o qra/qra64/qra64_subs.c
@@ -46,11 +47,14 @@ qra12_63_64_irr_b.o: qra/qracodes/qra12_63_64_irr_b.c
 qra13_64_64_irr_e.o: qra/qracodes/qra13_64_64_irr_e.c
 	gcc -c -O2 -o qra13_64_64_irr_e.o  qra/qracodes/qra13_64_64_irr_e.c
 
-
 qratest.exe: $(OBJS1)
 	$(FC) -o qratest.exe $(OBJS1) C:\JTSDK\fftw3f\libfftw3f-3.dll
 
+OBJS2 = tplt.o zplt64.o image.o
+tplt.exe: $(OBJS2)
+	$(FC) -o tplt.exe $(OBJS2)
+
 .PHONY : clean
 
 clean:
-	$(RM) *.o qratest.exe 
+	$(RM) *.o qratest.exe tplt.exe

lib/ana64.f90

@@ -0,0 +1,24 @@
+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
+  nfft2=nfft1/2
+  df1=12000.0/nfft1
+  fac=2.0/nfft1
+  c0(0:npts-1)=fac*dd(1:npts)
+  c0(npts:nfft1)=0.
+  call four2a(c0,nfft1,1,-1,1)             !Forward c2c FFT
+  c0(nfft2/2+1:nfft2)=0.
+  c0(0)=0.5*c0(0)
+  call four2a(c0,nfft2,1,1,1)              !Inverse c2c FFT; c0 is analytic sig
+
+  return
+end subroutine ana64

lib/qra64a.f90

@@ -58,12 +58,13 @@ subroutine qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
   endif
   naptype=maxaptype
 
-  maxf1=0
+  call ana64(dd,npts,c00)
+  
   call timer('sync64  ',0)
-  call sync64(dd,npts,nf1,nf2,nfqso,ntol,mode64,maxf1,dtx,f0,jpk0,kpk,sync,c00)
+  call sync64(c00,nf1,nf2,nfqso,ntol,mode64,dtx,f0,jpk0,sync,sync2,width)
   call timer('sync64  ',1)
   nfreq=nint(f0)
-  if((sync-3.4).lt.float(minsync)) go to 900
+  if((sync-3.4).lt.float(minsync) .or. sync2.lt.-2.0 .or.width.gt.240.0) go to 900
   a=0.
   a(1)=-f0
   npts2=npts/2
@@ -91,6 +92,7 @@ subroutine qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
      do iter=itz,0,-2
         b90=1.728**iter
         if(b90.gt.230.0) cycle
+        if(b90.lt.0.15*width) exit
         s3(1:LL*NN)=s3a(1:LL*NN)
         call timer('qra64_de',0)
         call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90,      &
@@ -148,5 +150,9 @@ subroutine qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
   endif
   call timer('qra64a  ',1)
 
+  write(83,4001) nutc,sync,sync2,b90,width,b90/width,nsnr,dtx,nfreq,decoded,nft-100
+4001 format(i4.4,5f6.1,i4,f6.2,i5,1x,a22,i3)
+
+
   return
 end subroutine qra64a

lib/qratest.f90

@@ -4,6 +4,7 @@ program qratest
   real dd(NMAX)
   character arg*8,mycall*12,hiscall*12,hisgrid*6,decoded*22
   character c*1
+  logical loop
 
   nargs=iargc()
   if(nargs.lt.1 .or. nargs.gt.4) then
@@ -12,6 +13,7 @@ program qratest
   endif
   call getarg(1,arg)
   read(arg,*) nfile
+  loop=arg(1:1).eq.'+'
   minsync0=-1
   nfqso0=-1
   ntol0=-1
@@ -46,7 +48,7 @@ program qratest
      if(mode64.eq.16) c='e'
      write(*,1000) ifile,c,nutc,nsnr,dtx,nfreq,decoded,nft-100,sync-3.4
 1000 format(i4,1x,a1,1x,i4.4,i4,f6.2,i5,1x,a22,i3,f6.2)
-     if(ifile.eq.nfile) exit
+     if(ifile.eq.nfile .and. (.not.loop)) exit
   enddo
 
 999 end program qratest

lib/spec64.f90

@@ -4,27 +4,39 @@ subroutine spec64(c0,npts2,mode64,jpk,s3,LL,NN)
   complex c0(0:360000)                       !Complex spectrum of dd()
   complex cs(0:NSPS-1)                       !Complex symbol spectrum
   real s3(LL,NN)                             !Synchronized symbol spectra
+  real xbase0(LL),xbase(LL)
 
-  nfft6=nsps
+  nfft=nsps
-  fac=1.0/nfft6
+  fac=1.0/nfft
-  do j=1,63
+  do j=1,NN
      jj=j+7                                  !Skip first Costas array
-     if(j.ge.32) jj=j+14                     !Skip middle Costas array
+     if(j.ge.33) jj=j+14                     !Skip middle Costas array
-     ja=jpk + (jj-1)*nfft6
+     ja=jpk + (jj-1)*nfft
-     jb=ja+nfft6-1
+     jb=ja+nfft-1
-     cs(0:nfft6-1)=fac*c0(ja:jb)
+     cs(0:nfft-1)=fac*c0(ja:jb)
-     call four2a(cs,nfft6,1,-1,1)
+     call four2a(cs,nfft,1,-1,1)
-     smax=0.
      do ii=1,LL
         i=ii-65
-        if(i.lt.0) i=i+nfft6
+        if(i.lt.0) i=i+nfft
         s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2
-        if(s3(ii,j).gt.smax) then
-           smax=s3(ii,j)
-           ipk=ii-65
-        endif
      enddo
   enddo
 
+  df=6000.0/nfft
+  do i=1,LL
+     call pctile(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape
+  enddo
+  
+  nh=25
+  xbase(1:nh-1)=sum(xbase0(1:nh-1))/(nh-1.0)
+  xbase(LL-nh+1:LL)=sum(xbase0(LL-nh+1:LL))/(nh-1.0)
+  do i=nh,LL-nh
+     xbase(i)=sum(xbase0(i-nh+1:i+nh))/(2*nh+1)  !Smoothed passband shape
+  enddo
+  
+  do i=1,LL
+     s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization
+  enddo
+
   return
 end subroutine spec64

lib/sync64.f90

@@ -1,12 +1,10 @@
-subroutine sync64(dd,npts,nf1,nf2,nfqso,ntol,mode64,maxf1,dtx,f0,jpk,kpk,   &
+subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,dtx,f0,jpk,sync,sync2,width)
-     sync,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
   real s1(0:NSPC-1)                          !Power spectrum of Costas 1
   real s2(0:NSPC-1)                          !Power spectrum of Costas 2
   real s3(0:NSPC-1)                          !Power spectrum of Costas 3
@@ -14,7 +12,7 @@ subroutine sync64(dd,npts,nf1,nf2,nfqso,ntol,mode64,maxf1,dtx,f0,jpk,kpk,   &
   real s0a(0:NSPC-1)                         !Best synchromized spectrum (saved)
   real s0b(0:NSPC-1)                         !tmp
   real s0c(0:NSPC-1)                         !tmp
-  logical old_qra_sync
+  real a(5)
   integer icos7(0:6)                         !Costas 7x7 tones
   integer ipk0(1)
   complex cc(0:NSPC-1)                       !Costas waveform
@@ -43,17 +41,6 @@ subroutine sync64(dd,npts,nf1,nf2,nfqso,ntol,mode64,maxf1,dtx,f0,jpk,kpk,   &
      mode64z=mode64
   endif
 
-  nfft1=672000
-  nfft2=nfft1/2
-  df1=12000.0/nfft1
-  fac=2.0/nfft1
-  c0(0:npts-1)=fac*dd(1:npts)
-  c0(npts:nfft1)=0.
-  call four2a(c0,nfft1,1,-1,1)             !Forward c2c FFT
-  c0(nfft2/2+1:nfft2)=0.
-  c0(0)=0.5*c0(0)
-  call four2a(c0,nfft2,1,1,1)              !Inverse c2c FFT; c0 is analytic sig
-  npts2=npts/2                            !Downsampled complex data length
   nfft3=NSPC
   nh3=nfft3/2
   df3=6000.0/nfft3
@@ -145,5 +132,35 @@ subroutine sync64(dd,npts,nf1,nf2,nfqso,ntol,mode64,maxf1,dtx,f0,jpk,kpk,   &
   write(17) ia,ib,s0a(ia:ib)                !Save data for red curve
   close(17)
 
+  nskip=50
+  call lorentzian(s0a(ia+nskip:ib-nskip),iz-2*nskip,a)
+  f0a=(a(3)+ia+49)*df3
+  w1=df3*a(4)
+  w2=2*nadd*df3
+  width=w1
+  if(w1.gt.1.2*w2) width=sqrt(w1**2 - w2**2)
+
+  sq=0.
+  do i=1,20
+     j=ia+nskip+1
+     k=ib-nskip-21+i
+     sq=sq + (s0a(j)-a(1))**2 + (s0a(k)-a(1))**2
+  enddo
+  rms=sqrt(sq/40.0)
+  sync2=10.0*log10(a(2)/rms)
+
+!  do i=1,iz-2*nskip
+!     x=i
+!     z=(x-a(3))/(0.5*a(4))
+!     yfit=a(1)
+!     if(abs(z).lt.3.0) then
+!        d=1.0 + z*z
+!        yfit=a(1) + a(2)*(1.0/d - 0.1)
+!     endif
+!     j=i+ia+49
+!     write(76,1110) j*df3,s0a(j),yfit
+!1110 format(3f10.3)
+!  enddo
+
   return
 end subroutine sync64