From fc040d435a8217fea77f39b513bce685846cd1f9 Mon Sep 17 00:00:00 2001
From: Joe Taylor <joe@princeton.edu>
Date: Wed, 4 Jan 2023 12:02:10 -0500
Subject: [PATCH] Sort out some stuff having to do with multiple decodes in Q65
mode.
---
lib/q65_decode.f90 | 16 +++++++++++--
lib/qra/q65/q65.f90 | 57 +++++++++++++++++++++++++++++++++++++--------
2 files changed, 61 insertions(+), 12 deletions(-)
diff --git a/lib/q65_decode.f90 b/lib/q65_decode.f90
index 515cdab6e..0834395bb 100644
--- a/lib/q65_decode.f90
+++ b/lib/q65_decode.f90
@@ -70,6 +70,7 @@ contains
character*80 fmt
integer*2 iwave(NMAX) !Raw data
real, allocatable :: dd(:) !Raw data
+ real f0decodes(100)
integer dat4(13) !Decoded message as 12 6-bit integers
integer dgen(13)
logical lclearave,lnewdat0,lapcqonly,unpk77_success
@@ -83,8 +84,13 @@ contains
call sec0(0,tdecode)
ndecodes=0
decodes=' '
+ f0decodes=0.
nfa=nfa0
nfb=nfb0
+ if(single_decode) then
+ nfa=nfqso-ntol
+ nfb=nfqso+ntol
+ endif
nqd=nqd0
lnewdat=lnewdat0
max_drift=max_drift0
@@ -169,8 +175,6 @@ contains
call q65_dec0(iavg,nutc,iwave,ntrperiod,nfqso,ntol,ndepth,lclearave, &
emedelay,xdt,f0,snr1,width,dat4,snr2,idec,stageno)
call timer('q65_dec0',1)
-! write(*,3001) '=a',nfqso,ntol,ndepth,xdt,f0,idec
-!3001 format(a2,3i5,f7.2,f7.1,i5)
if(idec.ge.0) then
dtdec=xdt !We have a q3 or q0 decode at nfqso
@@ -280,6 +284,7 @@ contains
if(idupe.eq.0) then
ndecodes=min(ndecodes+1,100)
decodes(ndecodes)=decoded
+ f0decodes(ndecodes)=f0dec
call q65_snr(dat4,dtdec,f0dec,mode_q65,nused,snr2)
nsnr=nint(snr2)
call this%callback(nutc,snr1,nsnr,dtdec,f0dec,decoded, &
@@ -317,6 +322,10 @@ contains
snr1=candidates(icand,1)
xdt= candidates(icand,2)
f0 = candidates(icand,3)
+ do i=1,ndecodes
+ fdiff=f0-f0decodes(i)
+ if(fdiff.gt.-baud*mode_q65 .and. fdiff.lt.65*baud*mode_q65) go to 800
+ enddo
jpk0=(xdt+1.0)*6000 !Index of nominal start of signal
if(ntrperiod.le.30) jpk0=(xdt+0.5)*6000 !For shortest sequences
if(jpk0.lt.0) jpk0=0
@@ -365,6 +374,7 @@ contains
if(idupe.eq.0) then
ndecodes=min(ndecodes+1,100)
decodes(ndecodes)=decoded
+ f0decodes(ndecodes)=f0dec
call q65_snr(dat4,dtdec,f0dec,mode_q65,nused,snr2)
nsnr=nint(snr2)
call this%callback(nutc,snr1,nsnr,dtdec,f0dec,decoded, &
@@ -373,6 +383,7 @@ contains
if(iand(ndepth,128).ne.0 .and. .not.lagain .and. &
int(abs(f0dec-nfqso)).le.ntol ) call q65_clravg !AutoClrAvg
call sec0(1,tdecode)
+ ios=1
! open(22,file=trim(data_dir)//'/q65_decodes.dat',status='unknown',&
! position='append',iostat=ios)
if(ios.eq.0) then
@@ -394,6 +405,7 @@ contains
endif
endif
endif
+800 continue
enddo ! icand
if(iavg.eq.0 .and.navg(iseq).ge.2 .and. iand(ndepth,16).ne.0) go to 50
900 return
diff --git a/lib/qra/q65/q65.f90 b/lib/qra/q65/q65.f90
index aa104bbfc..fa548e18f 100644
--- a/lib/qra/q65/q65.f90
+++ b/lib/qra/q65/q65.f90
@@ -148,15 +148,15 @@ subroutine q65_dec0(iavg,nutc,iwave,ntrperiod,nfqso,ntol,ndepth,lclearave, &
ii1=max(1,i0-64)
ii2=i0-65+LL
call pctile(s1(ii1:ii2,1:jz),ii2-ii1+1*jz,45,base)
- s1=s1/base
+! s1=s1/base
s1raw=s1
! Apply fast AGC to the symbol spectra
- s1max=20.0 !Empirical choice
- do j=1,jz !### Maybe wrong way? ###
- smax=maxval(s1(ii1:ii2,j))
- if(smax.gt.s1max) s1(ii1:ii2,j)=s1(ii1:ii2,j)*s1max/smax
- enddo
+! s1max=20.0 !Empirical choice
+! do j=1,jz !### Maybe wrong way? ###
+! smax=maxval(s1(ii1:ii2,j))
+! if(smax.gt.s1max) s1(ii1:ii2,j)=s1(ii1:ii2,j)*s1max/smax
+! enddo
dat4=0
if(ncw.gt.0 .and. iavg.le.1) then
@@ -219,7 +219,7 @@ subroutine q65_dec0(iavg,nutc,iwave,ntrperiod,nfqso,ntol,ndepth,lclearave, &
width=df*(i2-i1)
if(ncw.eq.0) ccf1=0.
call q65_write_red(iz,xdt,ccf2_avg,ccf2) !### Need this call for WSJT-X
-
+
if(idec.lt.0 .and. (iavg.eq.0 .or. iavg.eq.2)) then
call q65_dec_q012(s3,LL,snr2,dat4,idec,decoded)
endif
@@ -489,11 +489,14 @@ subroutine q65_ccf_22(s1,iz,jz,nfqso,ntol,ndepth,ntrperiod,iavg,ipk,jpk, &
real ccf2(iz) !Orange sync curve
real, allocatable :: xdt2(:)
real, allocatable :: s1avg(:)
+ real, allocatable :: ccf_lag(:)
integer, allocatable :: indx(:)
+ integer ipk1(1)
allocate(xdt2(iz))
allocate(s1avg(iz))
allocate(indx(iz))
+ allocate(ccf_lag(lag1:lag2))
ia=max(nfa,100)/df
ib=min(nfb,4900)/df
@@ -507,12 +510,15 @@ subroutine q65_ccf_22(s1,iz,jz,nfqso,ntol,ndepth,ntrperiod,iavg,ipk,jpk, &
s1avg(i)=sum(s1(i,1:jz))
enddo
+ call pctile(s1avg(ia:ib),ib-ia+1,40,base0)
ccfbest=0.
ibest=0
lagpk=0
lagbest=0
idrift_best=0
do i=ia,ib
+ stest=s1avg(i)/(1.015*base0)
+ if(stest.lt.1.4) cycle
ccfmax=0.
do lag=lag1,lag2
do idrift=-max_drift,max_drift
@@ -533,6 +539,40 @@ subroutine q65_ccf_22(s1,iz,jz,nfqso,ntol,ndepth,ntrperiod,iavg,ipk,jpk, &
endif
enddo ! idrift
enddo ! lag
+
+! Look at SNR on the "blue curve"
+ idrift=idrift_max
+ do lag=lag1,lag2
+ ccft=0.
+ do kk=1,22
+ k=isync(kk)
+ ii=i + nint(idrift*(k-43)/85.0)
+ if(ii.lt.1 .or. ii.gt.iz) cycle
+ n=NSTEP*(k-1) + 1
+ j=n+lag+j0
+ if(j.ge.1 .and. j.le.jz) ccft=ccft + s1(ii,j)
+ enddo ! kk
+ ccf_lag(lag)=ccft - (22.0/jz)*s1avg(i)
+ enddo
+
+ pk=maxval(ccf_lag)
+ ipk1=maxloc(ccf_lag)
+ lag0=ipk1(1)-1+lag1
+ xsum=0.
+ sq=0.
+ nsum=0
+ do j=lag1,lag2
+ if(abs(j-lag0).gt.7) then
+ xsum=xsum+ccf_lag(j)
+ sq=sq+ccf_lag(j)**2
+ nsum=nsum+1
+ endif
+ enddo
+ ave=xsum/nsum
+ rms=sqrt(sq/nsum - ave*ave)
+ snr=(pk-ave)/rms
+ if(snr.lt.5.0) cycle !Blue SNR must be at least 5
+
ccf2(i)=ccfmax
xdt2(i)=lagpk*dtstep
if(ccfmax.gt.ccfbest .and. abs(i*df-nfqso).le.ftol) then
@@ -554,8 +594,6 @@ subroutine q65_ccf_22(s1,iz,jz,nfqso,ntol,ndepth,ntrperiod,iavg,ipk,jpk, &
! Save parameters for best candidates
jzz=ib-ia+1
- call pctile(ccf2(ia:ib),jzz,40,base)
- ccf2=ccf2/base
call indexx(ccf2(ia:ib),jzz,indx)
ncand=0
maxcand=20
@@ -563,7 +601,6 @@ subroutine q65_ccf_22(s1,iz,jz,nfqso,ntol,ndepth,ntrperiod,iavg,ipk,jpk, &
k=jzz-j+1
if(k.lt.1 .or. k.gt.iz) cycle
i=indx(k)+ia-1
- if(ccf2(i).lt.3.3) exit !Candidate limit
f=i*df
i3=max(1, i-mode_q65)
i4=min(iz,i+mode_q65)