From d3ee8af01b318cf972d15f3c2ccaad44140e24b3 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 15 Dec 2016 18:42:33 +0000 Subject: [PATCH] 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 | 1 + lib/Makefile | 14 +++++++++----- lib/ana64.f90 | 24 ++++++++++++++++++++++++ lib/qra64a.f90 | 12 +++++++++--- lib/qratest.f90 | 4 +++- lib/spec64.f90 | 40 ++++++++++++++++++++++++++-------------- lib/sync64.f90 | 47 ++++++++++++++++++++++++++++++++--------------- 7 files changed, 104 insertions(+), 38 deletions(-) create mode 100644 lib/ana64.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 738110e6f..6e8cff36c 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/lib/Makefile b/lib/Makefile index a7fb74d15..00fae4dfe 100644 --- a/lib/Makefile +++ b/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 diff --git a/lib/ana64.f90 b/lib/ana64.f90 new file mode 100644 index 000000000..5681ee176 --- /dev/null +++ b/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 diff --git a/lib/qra64a.f90 b/lib/qra64a.f90 index b44262b4e..c9910cbae 100644 --- a/lib/qra64a.f90 +++ b/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 diff --git a/lib/qratest.f90 b/lib/qratest.f90 index 84e85600a..2b5de8c1c 100644 --- a/lib/qratest.f90 +++ b/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 diff --git a/lib/spec64.f90 b/lib/spec64.f90 index b44e5c12a..a0c0a04c7 100644 --- a/lib/spec64.f90 +++ b/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 - fac=1.0/nfft6 - do j=1,63 + nfft=nsps + fac=1.0/nfft + do j=1,NN jj=j+7 !Skip first Costas array - if(j.ge.32) jj=j+14 !Skip middle Costas array - ja=jpk + (jj-1)*nfft6 - jb=ja+nfft6-1 - cs(0:nfft6-1)=fac*c0(ja:jb) - call four2a(cs,nfft6,1,-1,1) - smax=0. + if(j.ge.33) jj=j+14 !Skip middle Costas array + ja=jpk + (jj-1)*nfft + jb=ja+nfft-1 + cs(0:nfft-1)=fac*c0(ja:jb) + call four2a(cs,nfft,1,-1,1) 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 diff --git a/lib/sync64.f90 b/lib/sync64.f90 index 4ac207dec..22ba4aa95 100644 --- a/lib/sync64.f90 +++ b/lib/sync64.f90 @@ -1,12 +1,10 @@ -subroutine sync64(dd,npts,nf1,nf2,nfqso,ntol,mode64,maxf1,dtx,f0,jpk,kpk, & - sync,c0) +subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,dtx,f0,jpk,sync,sync2,width) 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