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
This commit is contained in:
Joe Taylor 2016-12-15 18:42:33 +00:00
parent f26bc2b9cc
commit 27fe255e70
7 changed files with 104 additions and 38 deletions

View File

@ -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

View File

@ -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

24
lib/ana64.f90 Normal file
View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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