diff --git a/CMakeLists.txt b/CMakeLists.txt index d4f630f47..b910350da 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -496,6 +496,7 @@ set (wsjt_FSRCS lib/polyfit.f90 lib/prog_args.f90 lib/ps4.f90 + lib/q65_sync.f90 lib/qra64a.f90 lib/qra_loops.f90 lib/qra/q65/q65_ap.f90 @@ -529,7 +530,6 @@ set (wsjt_FSRCS lib/sync4.f90 lib/sync64.f90 lib/sync65.f90 - lib/sync_q65.f90 lib/ft4/getcandidates4.f90 lib/ft4/get_ft4_bitmetrics.f90 lib/ft8/sync8.f90 diff --git a/lib/q65_decode.f90 b/lib/q65_decode.f90 index 304ca45b4..819aba0d1 100644 --- a/lib/q65_decode.f90 +++ b/lib/q65_decode.f90 @@ -90,11 +90,22 @@ contains this%callback => callback if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso !Silence warning nFadingModel=1 + dgen=0 + call q65_enc(dgen,codewords) !Initialize Q65 +! nQSOprogress=3 !### + dat4=0 call timer('sync_q65',0) - call sync_q65(iwave,ntrperiod*12000,mode65,nQSOprogress,nsps,nfqso, & - ntol,xdt,f0,snr1,width) + call q65_sync(iwave,ntrperiod*12000,mode65,nQSOprogress,nsps,nfqso, & + ntol,xdt,f0,snr1,dat4,snr2,irc) 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 if(snr1.lt.2.8) go to 100 jpk0=(xdt+1.0)*6000 !### Is this OK? @@ -135,7 +146,7 @@ contains endif call timer('q65loops',0) 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) call timer('q65loops',1) snr2=snr2 + db(6912.0/nsps) diff --git a/lib/sync_q65.f90 b/lib/q65_sync.f90 similarity index 74% rename from lib/sync_q65.f90 rename to lib/q65_sync.f90 index 201d2a080..d81cbd5f8 100644 --- a/lib/sync_q65.f90 +++ b/lib/q65_sync.f90 @@ -1,11 +1,11 @@ -subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & - xdt,f0,snr1,width) +subroutine q65_sync(iwave,nmax,mode_q65,nQSOprogress,nsps,nfqso,ntol, & + xdt,f0,snr1,dat4,snr2,irc) ! Detect and align with the Q65 sync vector, returning time and frequency ! offsets and SNR estimate. ! 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 ! nfqso Target frequency (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 isync(22) !Indices of sync symbols integer itone(85) + integer codewords(63,64) + integer dat4(13) + integer ijpk(2) 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 :: ccf1(:) !CCF(freq) at best lag + real s3prob(0:63,63) !Symbol-value probabilities real sync(85) !sync vector 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 sync(1)/99.0/ save sync + LL=64*(2+mode_q65) nfft=nsps - df=12000.0/nfft !Freq resolution = 0.5*baud + df=12000.0/nfft !Freq resolution = baud istep=nsps/NSTEP iz=5000.0/df !Uppermost frequency bin, at 5000 Hz txt=85.0*nsps/12000.0 @@ -38,6 +44,7 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & ia=ntol/df allocate(s1(iz,jz)) + allocate(s3(-64:LL-65,63)) allocate(c0(0:nfft-1)) allocate(ccf(-ia:ia,-53:214)) 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 enddo ! 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 i0=nint(nfqso/df) !Target QSO frequency 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 s1max=20.0 !Empirical choice @@ -101,19 +109,9 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & enddo enddo - ic=ntol/df - ccfmax=0. - ipk=0 - 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 + ijpk=maxloc(ccf) + ipk=ijpk(1)-ia-1 + jpk=ijpk(2)-53-1 f0=nfqso + ipk*df xdt=jpk*dtstep @@ -129,28 +127,10 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & smax=ccf(ipk,jpk) snr1=smax/rms -! do j=lag1,lag2 -! write(55,3055) j,j*dtstep,ccf(ipk,j)/rms -!3055 format(i5,f8.3,f10.3) -! enddo +!###################################################################### +! Experimental: Try early list decoding via "Deep Likelihood". -! 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 -! "Deep Likelihood" decode attempt snr1a_best=0. do imsg=1,4 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.4) msg='CQ K9AN EN50' 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 k=1,85 j=j0 + NSTEP*(k-1) + 1 + lag @@ -172,22 +159,12 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & enddo enddo - ic=ntol/df - ccfmax=0. - ipk=0 - 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 + ijpk=maxloc(ccf) + ipk=ijpk(1)-ia-1 + jpk=ijpk(2)-53-1 f0a=nfqso + ipk*df xdta=jpk*dtstep - + sq=0. nsq=0 do j=lag1,lag2 @@ -205,27 +182,43 @@ subroutine sync_q65(iwave,nmax,mode65,nQSOprogress,nsps,nfqso,ntol, & xdta_best=xdta f0a_best=f0a endif -! write(57,3001) imsg,xdt,xdta,f0,f0a,snr1,snr1a -!3001 format(i1,6f8.2) + enddo ! imsg -! 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 xdt=xdta_best f0=f0a_best snr1=1.4*snr1a_best endif -! write(58,3006) xdta_best,f0a_best,snr1a_best,imsg_best -!3006 format(3f8.2,i3) + ia=i0+ipk-63 + 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 -end subroutine sync_q65 +end subroutine q65_sync diff --git a/lib/qra/q65/q65.c b/lib/qra/q65/q65.c index c91571ba6..c92335007 100644 --- a/lib/qra/q65/q65.c +++ b/lib/qra/q65/q65.c @@ -703,10 +703,11 @@ int q65_decode_fullaplist(q65_codec_ds *codec, maxllh = llh; maxcw = k; } + // printf("BBB %d %f\n",k,llh); // point to next codeword pCw+=nN; } - + q65_llh=maxllh; if (maxcw<0) // no llh larger than threshold found return Q65_DECODE_FAILED; diff --git a/lib/qra/q65/q65.h b/lib/qra/q65/q65.h index 2e764a32b..f48c40da9 100644 --- a/lib/qra/q65/q65.h +++ b/lib/qra/q65/q65.h @@ -39,7 +39,7 @@ // maximum number of weights for the fast-fading metric evaluation #define Q65_FASTFADING_MAXWEIGTHS 65 - +float q65_llh; typedef struct { const qracode *pQraCode; // qra code to be used by the codec float decoderEsNoMetric; // value for which we optimize the decoder metric diff --git a/lib/qra/q65/q65_loops.f90 b/lib/qra/q65/q65_loops.f90 index bb1033771..5c48f85a1 100644 --- a/lib/qra/q65/q65_loops.f90 +++ b/lib/qra/q65/q65_loops.f90 @@ -1,5 +1,5 @@ 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) use packjt77 @@ -91,7 +91,6 @@ subroutine q65_loops(c00,nutc,npts2,nsps,mode,mode_q65,nsubmode,nFadingModel, & ! b90=1.728**ibw b90=3.0**nbw if(b90.gt.230.0) cycle -! if(b90.lt.0.15*width) exit call timer('q65_intr',0) b90ts = b90/baud 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 codewords(1:63,4)=cw4 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) else 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) f1=f0 + 0.5*baud*ndf !### For tests only: - open(53,file='fort.53',status='unknown',position='append') - write(c77,1100) dat4(1:12),dat4(13)/2 -1100 format(12b6.6,b5.5) - call unpack77(c77,0,decoded,unpk77_success) !Unpack to get msgsent - m=nutc - if(nsps.ge.3600) m=100*m - ihr=m/10000 - imin=mod(m/100,100) - isec=mod(m,100) - hours=ihr + imin/60.0 + isec/3600.0 - write(53,3053) m,hours,ndf,ndt,nbw,ndist,irc,iaptype,kavg,snr1, & - xdt1,f1,snr2,trim(decoded) -3053 format(i6.6,f8.4,4i3,i4,2i3,f6.1,f6.2,f7.1,f6.1,1x,a) - close(53) +! open(53,file='fort.53',status='unknown',position='append') +! write(c77,1100) dat4(1:12),dat4(13)/2 +!1100 format(12b6.6,b5.5) +! call unpack77(c77,0,decoded,unpk77_success) !Unpack to get msgsent +! m=nutc +! if(nsps.ge.3600) m=100*m +! ihr=m/10000 +! imin=mod(m/100,100) +! isec=mod(m,100) +! hours=ihr + imin/60.0 + isec/3600.0 +! write(53,3053) m,hours,ndf,ndt,nbw,ndist,irc,iaptype,kavg,snr1, & +! xdt1,f1,snr2,trim(decoded) +!3053 format(i6.6,f8.4,4i3,i4,2i3,f6.1,f6.2,f7.1,f6.1,1x,a) +! close(53) !### nsave=0 s3avg=0. diff --git a/lib/qra/q65/q65_subs.c b/lib/qra/q65/q65_subs.c index 9f29da01c..e55fe927d 100644 --- a/lib/qra/q65/q65_subs.c +++ b/lib/qra/q65/q65_subs.c @@ -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[], - int* ncw, float* esnodb0, int xdec[], int* rc0) + int* ncw, float* esnodb0, int xdec[], float* plog, int* rc0) { /* Input: s3[LL,NN] Symbol spectra * s3prob[LL,NN] Symbol-value intrinsic probabilities @@ -128,6 +128,7 @@ void q65_dec_fullaplist_(float s3[], float s3prob[], int codewords[], float esnodb; rc = q65_decode_fullaplist(&codec,ydec,xdec,s3prob,codewords,*ncw); + *plog=q65_llh; *rc0=rc; // rc = -1: Invalid params diff --git a/lib/qra/q65/q65sim.f90 b/lib/qra/q65/q65sim.f90 index f99ab4f1d..241d0adb8 100644 --- a/lib/qra/q65/q65sim.f90 +++ b/lib/qra/q65/q65sim.f90 @@ -193,21 +193,21 @@ program q65sim write(10) h,iwave(1:npts) !Save the .wav file close(10) - if(lsync) then - cd=' ' - if(ifile.eq.nfiles) cd='d' - nfqso=nint(f0) - ntol=100 - call sync_q65(iwave,npts,mode65,nsps,nfqso,ntol,xdt2,f02,snr2) - terr=1.01/(8.0*baud) - ferr=1.01*mode65*baud - if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1 - open(40,file='sync65.out',status='unknown',position='append') - write(40,1030) ifile,65,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, & - snr2,nsync,cd -1030 format(i4,i3,1x,a1,2f7.1,f7.2,2f8.1,i5,1x,a1) - close(40) - endif +! if(lsync) then +! cd=' ' +! if(ifile.eq.nfiles) cd='d' +! nfqso=nint(f0) +! ntol=100 +! call q65_sync(iwave,npts,mode65,nsps,nfqso,ntol,xdt2,f02,snr2) +! terr=1.01/(8.0*baud) +! ferr=1.01*mode65*baud +! if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1 +! open(40,file='sync65.out',status='unknown',position='append') +! write(40,1030) ifile,65,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, & +! snr2,nsync,cd +!1030 format(i4,i3,1x,a1,2f7.1,f7.2,2f8.1,i5,1x,a1) +! close(40) +! endif enddo if(lsync) write(*,1040) snrdb,nfiles,nsync 1040 format('SNR:',f6.1,' nfiles:',i5,' nsynced:',i5)