Temporary save, much work in progress.

This commit is contained in:
Joe Taylor 2020-11-30 09:52:47 -05:00
parent 64516e6abb
commit 9ff6f5b4d3
8 changed files with 114 additions and 108 deletions

View File

@ -496,6 +496,7 @@ set (wsjt_FSRCS
lib/polyfit.f90 lib/polyfit.f90
lib/prog_args.f90 lib/prog_args.f90
lib/ps4.f90 lib/ps4.f90
lib/q65_sync.f90
lib/qra64a.f90 lib/qra64a.f90
lib/qra_loops.f90 lib/qra_loops.f90
lib/qra/q65/q65_ap.f90 lib/qra/q65/q65_ap.f90
@ -529,7 +530,6 @@ set (wsjt_FSRCS
lib/sync4.f90 lib/sync4.f90
lib/sync64.f90 lib/sync64.f90
lib/sync65.f90 lib/sync65.f90
lib/sync_q65.f90
lib/ft4/getcandidates4.f90 lib/ft4/getcandidates4.f90
lib/ft4/get_ft4_bitmetrics.f90 lib/ft4/get_ft4_bitmetrics.f90
lib/ft8/sync8.f90 lib/ft8/sync8.f90

View File

@ -90,11 +90,22 @@ contains
this%callback => callback this%callback => callback
if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso !Silence warning if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso !Silence warning
nFadingModel=1 nFadingModel=1
dgen=0
call q65_enc(dgen,codewords) !Initialize Q65
! nQSOprogress=3 !###
dat4=0
call timer('sync_q65',0) call timer('sync_q65',0)
call sync_q65(iwave,ntrperiod*12000,mode65,nQSOprogress,nsps,nfqso, & call q65_sync(iwave,ntrperiod*12000,mode65,nQSOprogress,nsps,nfqso, &
ntol,xdt,f0,snr1,width) ntol,xdt,f0,snr1,dat4,snr2,irc)
call timer('sync_q65',1) call timer('sync_q65',1)
write(55,3055) nutc,xdt,f0,snr1,snr2,irc
3055 format(i4.4,4f9.2,i5)
if(irc.ge.0) then
xdt1=xdt
f1=f0
go to 100
endif
irc=-9 irc=-9
if(snr1.lt.2.8) go to 100 if(snr1.lt.2.8) go to 100
jpk0=(xdt+1.0)*6000 !### Is this OK? jpk0=(xdt+1.0)*6000 !### Is this OK?
@ -135,7 +146,7 @@ contains
endif endif
call timer('q65loops',0) call timer('q65loops',0)
call q65_loops(c00,nutc,npts/2,nsps/2,nmode,mode65,nsubmode, & call q65_loops(c00,nutc,npts/2,nsps/2,nmode,mode65,nsubmode, &
nFadingModel,ndepth,jpk0,xdt,f0,width,iaptype,apmask,apsymbols, & nFadingModel,ndepth,jpk0,xdt,f0,iaptype,apmask,apsymbols, &
codewords,snr1,xdt1,f1,snr2,irc,dat4) codewords,snr1,xdt1,f1,snr2,irc,dat4)
call timer('q65loops',1) call timer('q65loops',1)
snr2=snr2 + db(6912.0/nsps) snr2=snr2 + db(6912.0/nsps)

View File

@ -1,11 +1,11 @@
subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & subroutine q65_sync(iwave,nmax,mode_q65,nQSOprogress,nsps,nfqso,ntol, &
xdt,f0,snr1,width) xdt,f0,snr1,dat4,snr2,irc)
! Detect and align with the Q65 sync vector, returning time and frequency ! Detect and align with the Q65 sync vector, returning time and frequency
! offsets and SNR estimate. ! offsets and SNR estimate.
! Input: iwave(0:nmax-1) Raw data ! Input: iwave(0:nmax-1) Raw data
! mode65 Tone spacing 1 2 4 8 16 (A-E) ! mode_q65 Tone spacing 1 2 4 8 16 (A-E)
! nsps Samples per symbol at 12000 Sa/s ! nsps Samples per symbol at 12000 Sa/s
! nfqso Target frequency (Hz) ! nfqso Target frequency (Hz)
! ntol Search range around nfqso (Hz) ! ntol Search range around nfqso (Hz)
@ -19,17 +19,23 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
integer*2 iwave(0:nmax-1) !Raw data integer*2 iwave(0:nmax-1) !Raw data
integer isync(22) !Indices of sync symbols integer isync(22) !Indices of sync symbols
integer itone(85) integer itone(85)
integer codewords(63,64)
integer dat4(13)
integer ijpk(2)
real, allocatable :: s1(:,:) !Symbol spectra, 1/8-symbol steps real, allocatable :: s1(:,:) !Symbol spectra, 1/8-symbol steps
real, allocatable :: s3(:,:) !Data-symbol energies s3(LL,63)
real, allocatable :: ccf(:,:) !CCF(freq,lag) real, allocatable :: ccf(:,:) !CCF(freq,lag)
real, allocatable :: ccf1(:) !CCF(freq) at best lag real, allocatable :: ccf1(:) !CCF(freq) at best lag
real s3prob(0:63,63) !Symbol-value probabilities
real sync(85) !sync vector real sync(85) !sync vector
complex, allocatable :: c0(:) !Complex spectrum of symbol complex, allocatable :: c0(:) !Complex spectrum of symbol
data isync/1,9,12,13,15,22,23,26,27,33,35,38,46,50,55,60,62,66,69,74,76,85/ data isync/1,9,12,13,15,22,23,26,27,33,35,38,46,50,55,60,62,66,69,74,76,85/
data sync(1)/99.0/ data sync(1)/99.0/
save sync save sync
LL=64*(2+mode_q65)
nfft=nsps nfft=nsps
df=12000.0/nfft !Freq resolution = 0.5*baud df=12000.0/nfft !Freq resolution = baud
istep=nsps/NSTEP istep=nsps/NSTEP
iz=5000.0/df !Uppermost frequency bin, at 5000 Hz iz=5000.0/df !Uppermost frequency bin, at 5000 Hz
txt=85.0*nsps/12000.0 txt=85.0*nsps/12000.0
@ -38,6 +44,7 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
ia=ntol/df ia=ntol/df
allocate(s1(iz,jz)) allocate(s1(iz,jz))
allocate(s3(-64:LL-65,63))
allocate(c0(0:nfft-1)) allocate(c0(0:nfft-1))
allocate(ccf(-ia:ia,-53:214)) allocate(ccf(-ia:ia,-53:214))
allocate(ccf1(-ia:ia)) allocate(ccf1(-ia:ia))
@ -66,12 +73,13 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
s1(i,j)=real(c0(i))**2 + aimag(c0(i))**2 s1(i,j)=real(c0(i))**2 + aimag(c0(i))**2
enddo enddo
! For large Doppler spreads, should we smooth the spectra here? ! For large Doppler spreads, should we smooth the spectra here?
call smo121(s1(1:iz,j),iz) ! call smo121(s1(1:iz,j),iz)
enddo enddo
i0=nint(nfqso/df) !Target QSO frequency i0=nint(nfqso/df) !Target QSO frequency
call pctile(s1(i0-64:i0+192,1:jz),129*jz,40,base) call pctile(s1(i0-64:i0+192,1:jz),129*jz,40,base)
s1=s1/base - 1.0 ! s1=s1/base - 1.0
s1=s1/base
! Apply fast AGC ! Apply fast AGC
s1max=20.0 !Empirical choice s1max=20.0 !Empirical choice
@ -101,19 +109,9 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
enddo enddo
enddo enddo
ic=ntol/df ijpk=maxloc(ccf)
ccfmax=0. ipk=ijpk(1)-ia-1
ipk=0 jpk=ijpk(2)-53-1
jpk=0
do i=-ic,ic
do j=lag1,lag2
if(ccf(i,j).gt.ccfmax) then
ipk=i
jpk=j
ccfmax=ccf(i,j)
endif
enddo
enddo
f0=nfqso + ipk*df f0=nfqso + ipk*df
xdt=jpk*dtstep xdt=jpk*dtstep
@ -129,28 +127,10 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
smax=ccf(ipk,jpk) smax=ccf(ipk,jpk)
snr1=smax/rms snr1=smax/rms
! do j=lag1,lag2 !######################################################################
! write(55,3055) j,j*dtstep,ccf(ipk,j)/rms ! Experimental: Try early list decoding via "Deep Likelihood".
!3055 format(i5,f8.3,f10.3)
! enddo
! do i=-ia,ia
! write(56,3056) i*df,ccf(i,jpk)/rms
!3056 format(2f10.3)
! enddo
! flush(56)
ccf1=ccf(-ia:ia,jpk)
acf0=dot_product(ccf1,ccf1)
do i=1,ia
acf=dot_product(ccf1,cshift(ccf1,i))
if(acf.le.0.5*acf0) exit
enddo
width=i*1.414*df
!### Experimental:
if(nQSOprogress.lt.1) go to 900 if(nQSOprogress.lt.1) go to 900
! "Deep Likelihood" decode attempt
snr1a_best=0. snr1a_best=0.
do imsg=1,4 do imsg=1,4
ccf=0. ccf=0.
@ -159,7 +139,14 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
if(imsg.eq.3) msg='K1ABC W9XYZ 73' if(imsg.eq.3) msg='K1ABC W9XYZ 73'
if(imsg.eq.4) msg='CQ K9AN EN50' if(imsg.eq.4) msg='CQ K9AN EN50'
call genq65(msg,0,msgsent,itone,i3,n3) call genq65(msg,0,msgsent,itone,i3,n3)
j=0
do k=1,85
if(sync(k)>0.) cycle
j=j+1
codewords(j,imsg)=itone(k) - 1
enddo
! Compute 2D ccf using all 85 symbols in the list message
do lag=lag1,lag2 do lag=lag1,lag2
do k=1,85 do k=1,85
j=j0 + NSTEP*(k-1) + 1 + lag j=j0 + NSTEP*(k-1) + 1 + lag
@ -172,22 +159,12 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
enddo enddo
enddo enddo
ic=ntol/df ijpk=maxloc(ccf)
ccfmax=0. ipk=ijpk(1)-ia-1
ipk=0 jpk=ijpk(2)-53-1
jpk=0
do i=-ic,ic
do j=lag1,lag2
if(ccf(i,j).gt.ccfmax) then
ipk=i
jpk=j
ccfmax=ccf(i,j)
endif
enddo
enddo
f0a=nfqso + ipk*df f0a=nfqso + ipk*df
xdta=jpk*dtstep xdta=jpk*dtstep
sq=0. sq=0.
nsq=0 nsq=0
do j=lag1,lag2 do j=lag1,lag2
@ -205,27 +182,43 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, &
xdta_best=xdta xdta_best=xdta
f0a_best=f0a f0a_best=f0a
endif endif
! write(57,3001) imsg,xdt,xdta,f0,f0a,snr1,snr1a enddo ! imsg
!3001 format(i1,6f8.2)
! do j=lag1,lag2
! write(55,3055) j,j*dtstep,ccf(ipk,j)/rms
!3055 format(i5,f8.3,f10.3)
! enddo
! do i=-ia,ia
! write(56,3056) i*df,ccf(i,jpk)/rms
!3056 format(2f10.3)
! enddo
enddo
if(snr1a_best.gt.2.0) then if(snr1a_best.gt.2.0) then
xdt=xdta_best xdt=xdta_best
f0=f0a_best f0=f0a_best
snr1=1.4*snr1a_best snr1=1.4*snr1a_best
endif endif
! write(58,3006) xdta_best,f0a_best,snr1a_best,imsg_best ia=i0+ipk-63
!3006 format(3f8.2,i3) ib=ia+LL-1
j=j0+jpk-5
n=0
do k=1,85
j=j+8
if(sync(k).gt.0.0) then
cycle
endif
n=n+1
s3(-64:LL-65,n)=s1(ia:ib,j)
enddo
nsubmode=0
nFadingModel=1
baud=12000.0/nsps
dat4=0
irc=-2
do ibw=0,10
b90=1.72**ibw
call q65_intrinsics_ff(s3,nsubmode,b90/baud,nFadingModel,s3prob)
call q65_dec_fullaplist(s3,s3prob,codewords,4,esnodb,dat4,plog,irc)
if(irc.ge.0) then
xdt=xdta_best
f0=f0a_best
snr2=esnodb - db(2500.0/baud)
exit
endif
enddo
900 return 900 return
end subroutine sync_q65 end subroutine q65_sync

View File

@ -703,10 +703,11 @@ int q65_decode_fullaplist(q65_codec_ds *codec,
maxllh = llh; maxllh = llh;
maxcw = k; maxcw = k;
} }
// printf("BBB %d %f\n",k,llh);
// point to next codeword // point to next codeword
pCw+=nN; pCw+=nN;
} }
q65_llh=maxllh;
if (maxcw<0) // no llh larger than threshold found if (maxcw<0) // no llh larger than threshold found
return Q65_DECODE_FAILED; return Q65_DECODE_FAILED;

View File

@ -39,7 +39,7 @@
// maximum number of weights for the fast-fading metric evaluation // maximum number of weights for the fast-fading metric evaluation
#define Q65_FASTFADING_MAXWEIGTHS 65 #define Q65_FASTFADING_MAXWEIGTHS 65
float q65_llh;
typedef struct { typedef struct {
const qracode *pQraCode; // qra code to be used by the codec const qracode *pQraCode; // qra code to be used by the codec
float decoderEsNoMetric; // value for which we optimize the decoder metric float decoderEsNoMetric; // value for which we optimize the decoder metric

View File

@ -1,5 +1,5 @@
subroutine q65_loops(c00,nutc,npts2,nsps,mode,mode_q65,nsubmode,nFadingModel, & subroutine q65_loops(c00,nutc,npts2,nsps,mode,mode_q65,nsubmode,nFadingModel, &
ndepth,jpk0,xdt0,f0,width,iaptype,APmask,APsymbols,codewords,snr1, & ndepth,jpk0,xdt0,f0,iaptype,APmask,APsymbols,codewords,snr1, &
xdt1,f1,snr2,irc,dat4) xdt1,f1,snr2,irc,dat4)
use packjt77 use packjt77
@ -91,7 +91,6 @@ subroutine q65_loops(c00,nutc,npts2,nsps,mode,mode_q65,nsubmode,nFadingModel, &
! b90=1.728**ibw ! b90=1.728**ibw
b90=3.0**nbw b90=3.0**nbw
if(b90.gt.230.0) cycle if(b90.gt.230.0) cycle
! if(b90.lt.0.15*width) exit
call timer('q65_intr',0) call timer('q65_intr',0)
b90ts = b90/baud b90ts = b90/baud
call q65_intrinsics_ff(s3,nsubmode,b90ts,nFadingModel,s3prob) call q65_intrinsics_ff(s3,nsubmode,b90ts,nFadingModel,s3prob)
@ -99,7 +98,8 @@ subroutine q65_loops(c00,nutc,npts2,nsps,mode,mode_q65,nsubmode,nFadingModel, &
if(iaptype.eq.4) then if(iaptype.eq.4) then
codewords(1:63,4)=cw4 codewords(1:63,4)=cw4
call timer('q65_apli',0) call timer('q65_apli',0)
call q65_dec_fullaplist(s3,s3prob,codewords,4,esnodb,dat4,irc) call q65_dec_fullaplist(s3,s3prob,codewords,4,esnodb, &
dat4,plog,irc)
call timer('q65_apli',1) call timer('q65_apli',1)
else else
call timer('q65_dec ',0) call timer('q65_dec ',0)
@ -140,20 +140,20 @@ subroutine q65_loops(c00,nutc,npts2,nsps,mode,mode_q65,nsubmode,nFadingModel, &
xdt1=xdt0 + nsps*ndt/(16.0*6000.0) xdt1=xdt0 + nsps*ndt/(16.0*6000.0)
f1=f0 + 0.5*baud*ndf f1=f0 + 0.5*baud*ndf
!### For tests only: !### For tests only:
open(53,file='fort.53',status='unknown',position='append') ! open(53,file='fort.53',status='unknown',position='append')
write(c77,1100) dat4(1:12),dat4(13)/2 ! write(c77,1100) dat4(1:12),dat4(13)/2
1100 format(12b6.6,b5.5) !1100 format(12b6.6,b5.5)
call unpack77(c77,0,decoded,unpk77_success) !Unpack to get msgsent ! call unpack77(c77,0,decoded,unpk77_success) !Unpack to get msgsent
m=nutc ! m=nutc
if(nsps.ge.3600) m=100*m ! if(nsps.ge.3600) m=100*m
ihr=m/10000 ! ihr=m/10000
imin=mod(m/100,100) ! imin=mod(m/100,100)
isec=mod(m,100) ! isec=mod(m,100)
hours=ihr + imin/60.0 + isec/3600.0 ! hours=ihr + imin/60.0 + isec/3600.0
write(53,3053) m,hours,ndf,ndt,nbw,ndist,irc,iaptype,kavg,snr1, & ! write(53,3053) m,hours,ndf,ndt,nbw,ndist,irc,iaptype,kavg,snr1, &
xdt1,f1,snr2,trim(decoded) ! xdt1,f1,snr2,trim(decoded)
3053 format(i6.6,f8.4,4i3,i4,2i3,f6.1,f6.2,f7.1,f6.1,1x,a) !3053 format(i6.6,f8.4,4i3,i4,2i3,f6.1,f6.2,f7.1,f6.1,1x,a)
close(53) ! close(53)
!### !###
nsave=0 nsave=0
s3avg=0. s3avg=0.

View File

@ -111,7 +111,7 @@ void q65_dec_(float s3[], float s3prob[], int APmask[], int APsymbols[],
} }
void q65_dec_fullaplist_(float s3[], float s3prob[], int codewords[], void q65_dec_fullaplist_(float s3[], float s3prob[], int codewords[],
int* ncw, float* esnodb0, int xdec[], int* rc0) int* ncw, float* esnodb0, int xdec[], float* plog, int* rc0)
{ {
/* Input: s3[LL,NN] Symbol spectra /* Input: s3[LL,NN] Symbol spectra
* s3prob[LL,NN] Symbol-value intrinsic probabilities * s3prob[LL,NN] Symbol-value intrinsic probabilities
@ -128,6 +128,7 @@ void q65_dec_fullaplist_(float s3[], float s3prob[], int codewords[],
float esnodb; float esnodb;
rc = q65_decode_fullaplist(&codec,ydec,xdec,s3prob,codewords,*ncw); rc = q65_decode_fullaplist(&codec,ydec,xdec,s3prob,codewords,*ncw);
*plog=q65_llh;
*rc0=rc; *rc0=rc;
// rc = -1: Invalid params // rc = -1: Invalid params

View File

@ -193,21 +193,21 @@ program q65sim
write(10) h,iwave(1:npts) !Save the .wav file write(10) h,iwave(1:npts) !Save the .wav file
close(10) close(10)
if(lsync) then ! if(lsync) then
cd=' ' ! cd=' '
if(ifile.eq.nfiles) cd='d' ! if(ifile.eq.nfiles) cd='d'
nfqso=nint(f0) ! nfqso=nint(f0)
ntol=100 ! ntol=100
call sync_q65(iwave,npts,mode65,nsps,nfqso,ntol,xdt2,f02,snr2) ! call q65_sync(iwave,npts,mode65,nsps,nfqso,ntol,xdt2,f02,snr2)
terr=1.01/(8.0*baud) ! terr=1.01/(8.0*baud)
ferr=1.01*mode65*baud ! ferr=1.01*mode65*baud
if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1 ! if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1
open(40,file='sync65.out',status='unknown',position='append') ! open(40,file='sync65.out',status='unknown',position='append')
write(40,1030) ifile,65,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, & ! write(40,1030) ifile,65,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, &
snr2,nsync,cd ! snr2,nsync,cd
1030 format(i4,i3,1x,a1,2f7.1,f7.2,2f8.1,i5,1x,a1) !1030 format(i4,i3,1x,a1,2f7.1,f7.2,2f8.1,i5,1x,a1)
close(40) ! close(40)
endif ! endif
enddo enddo
if(lsync) write(*,1040) snrdb,nfiles,nsync if(lsync) write(*,1040) snrdb,nfiles,nsync
1040 format('SNR:',f6.1,' nfiles:',i5,' nsynced:',i5) 1040 format('SNR:',f6.1,' nfiles:',i5,' nsynced:',i5)