diff --git a/q65w/libq65/CMakeLists.txt b/q65w/libq65/CMakeLists.txt index e0ac18ede..b5c21ee62 100644 --- a/q65w/libq65/CMakeLists.txt +++ b/q65w/libq65/CMakeLists.txt @@ -1,6 +1,6 @@ set (libq65_FSRCS # Modules come first: - wideband_sync.f90 +wideband_sync.f90 # Non-module Fortran routines: astro.f90 @@ -14,14 +14,16 @@ set (libq65_FSRCS four2a.f90 ftninit.f90 ftnquit.f90 - q65b.f90 geocentric.f90 + getcand2.f90 grid2deg.f90 indexx.f90 lorentzian.f90 moon2.f90 moondop.f90 + q65b.f90 q65c.f90 + q65_sync.f90 q65wa.f90 recvpkt.f90 sun.f90 diff --git a/q65w/libq65/getcand2.f90 b/q65w/libq65/getcand2.f90 new file mode 100644 index 000000000..7ebcc35aa --- /dev/null +++ b/q65w/libq65/getcand2.f90 @@ -0,0 +1,62 @@ +subroutine getcand2(ss,savg0,nts_q65,cand,ncand) + + use wideband_sync +! parameter(NFFT=32768) + real ss(322,NFFT) + real savg0(NFFT),savg(NFFT) + integer ipk1(1) + logical sync_ok + type(candidate) :: cand(MAX_CANDIDATES) + data nseg/16/,npct/40/ + + savg=savg0 + nlen=NFFT/nseg + do iseg=1,nseg + ja=(iseg-1)*nlen + 1 + jb=ja + nlen - 1 + call pctile(savg(ja),nlen,npct,base) + savg(ja:jb)=savg(ja:jb)/(1.015*base) + savg0(ja:jb)=savg0(ja:jb)/(1.015*base) + enddo + + df=96000.0/NFFT + bw=65*nts_q65*1.666666667 + nbw=bw/df + 1 + nb0=2*nts_q65 + smin=1.4 + nguard=5 + + j=0 + sync(1:NFFT)%ccfmax=0. + + do i=1,NFFT-nbw-nguard + if(savg(i).lt.smin) cycle + spk=maxval(savg(i:i+nb0)) + ipk1=maxloc(savg(i:i+nb0)) + i0=ipk1(1) + i - 1 + fpk=0.001*i0*df +! Check to see if sync tone is present. + call q65_sync(ss,i0,nts_q65,sync_ok,snr_sync,xdt) + if(.not.sync_ok) cycle + j=j+1 +! write(73,3073) j,fpk+32.0-2.270,snr_sync,xdt +!3073 format(i3,3f10.3) + cand(j)%f=fpk + cand(j)%xdt=xdt + cand(j)%snr=snr_sync + cand(j)%iflip=0 + sync(i0)%ccfmax=snr_sync + ia=min(i,i0-nguard) + ib=i0+nbw+nguard + savg(ia:ib)=0. + if(j.ge.30) exit + enddo + ncand=j + +! do i=1,NFFT +! write(72,3072) i,0.001*i*df+32.0,savg0(i),savg(i),sync(i)%ccfmax +!3072 format(i6,f15.6,3f15.3) +! enddo + + return +end subroutine getcand2 diff --git a/q65w/libq65/q65_sync.f90 b/q65w/libq65/q65_sync.f90 new file mode 100644 index 000000000..bd5f9aa77 --- /dev/null +++ b/q65w/libq65/q65_sync.f90 @@ -0,0 +1,56 @@ +subroutine q65_sync(ss,i0,nts_q65,sync_ok,snr,xdt) + + parameter (NFFT=32768) + parameter (LAGMAX=33) + real ss(322,NFFT) + real ccf(0:LAGMAX) + logical sync_ok + logical first + integer isync(22),ipk(1) + +! Q65 sync symbols + 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 first/.true./ + save first,isync + + tstep=2048.0/11025.0 !0.185760 s: 0.5*tsym_jt65, 0.3096*tsym_q65 + if(first) then + fac=0.6/tstep !3.230 + do i=1,22 !Expand the Q65 sync stride + isync(i)=nint((isync(i)-1)*fac) + 1 + enddo + first=.false. + endif + + m=nts_q65/2 + ccf=0. + do lag=0,LAGMAX + do j=1,22 !Test for Q65 sync + k=isync(j) + lag +! ccf=ccf + ss(k,i0) + ss(k+1,i0) + ss(k+2,i0) + ccf(lag)=ccf(lag) + sum(ss(k,i0-m:i0+m)) + sum(ss(k+1,i0-m:i0+m)) & + + sum(ss(k+2,i0-m:i0+m)) + enddo + enddo + ccfmax=maxval(ccf) + ipk=maxloc(ccf) + lagbest=ipk(1)-1 + xdt=lagbest*tstep - 1.0 + + xsum=0. + sq=0. + nsum=0 + do i=0,lagmax + if(abs(i-lagbest).gt.2) then + xsum=xsum+ccf(i) + sq=sq+ccf(i)**2 + nsum=nsum+1 + endif + enddo + ave=xsum/nsum + rms=sqrt(sq/nsum - ave*ave) + snr=(ccfmax-ave)/rms + sync_ok=snr.ge.5.0 + + return +end subroutine q65_sync diff --git a/q65w/libq65/q65wa.f90 b/q65w/libq65/q65wa.f90 index 0a949fc1b..af0ad6b2d 100644 --- a/q65w/libq65/q65wa.f90 +++ b/q65w/libq65/q65wa.f90 @@ -86,123 +86,3 @@ subroutine q65wa(dd,ss,savg,newdat,nutc,fcenter,ntol,nfa,nfb, & return end subroutine q65wa - -subroutine getcand2(ss,savg0,nts_q65,cand,ncand) - - use wideband_sync -! parameter(NFFT=32768) - real ss(322,NFFT) - real savg0(NFFT),savg(NFFT) - integer ipk1(1) - logical sync_ok - type(candidate) :: cand(MAX_CANDIDATES) - data nseg/16/,npct/40/ - - savg=savg0 - nlen=NFFT/nseg - do iseg=1,nseg - ja=(iseg-1)*nlen + 1 - jb=ja + nlen - 1 - call pctile(savg(ja),nlen,npct,base) - savg(ja:jb)=savg(ja:jb)/(1.015*base) - savg0(ja:jb)=savg0(ja:jb)/(1.015*base) - enddo - - df=96000.0/NFFT - bw=65*nts_q65*1.666666667 - nbw=bw/df + 1 - nb0=2*nts_q65 - smin=1.4 - nguard=5 - - j=0 - sync(1:NFFT)%ccfmax=0. - - do i=1,NFFT-nbw-nguard - if(savg(i).lt.smin) cycle - spk=maxval(savg(i:i+nb0)) - ipk1=maxloc(savg(i:i+nb0)) - i0=ipk1(1) + i - 1 - fpk=0.001*i0*df -! Check to see if sync tone is present. - call q65_sync(ss,i0,nts_q65,sync_ok,snr_sync,xdt) - if(.not.sync_ok) cycle - j=j+1 -! write(73,3073) j,fpk+32.0-2.270,snr_sync,xdt -!3073 format(i3,3f10.3) - cand(j)%f=fpk - cand(j)%xdt=xdt - cand(j)%snr=snr_sync - cand(j)%iflip=0 - sync(i0)%ccfmax=snr_sync - ia=min(i,i0-nguard) - ib=i0+nbw+nguard - savg(ia:ib)=0. - if(j.ge.30) exit - enddo - ncand=j - -! do i=1,NFFT -! write(72,3072) i,0.001*i*df+32.0,savg0(i),savg(i),sync(i)%ccfmax -!3072 format(i6,f15.6,3f15.3) -! enddo - - return -end subroutine getcand2 - -subroutine q65_sync(ss,i0,nts_q65,sync_ok,snr,xdt) - - parameter (NFFT=32768) - parameter (LAGMAX=33) - real ss(322,NFFT) - real ccf(0:LAGMAX) - logical sync_ok - logical first - integer isync(22),ipk(1) - -! Q65 sync symbols - 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 first/.true./ - save first,isync - - tstep=2048.0/11025.0 !0.185760 s: 0.5*tsym_jt65, 0.3096*tsym_q65 - if(first) then - fac=0.6/tstep !3.230 - do i=1,22 !Expand the Q65 sync stride - isync(i)=nint((isync(i)-1)*fac) + 1 - enddo - first=.false. - endif - - m=nts_q65/2 - ccf=0. - do lag=0,LAGMAX - do j=1,22 !Test for Q65 sync - k=isync(j) + lag -! ccf=ccf + ss(k,i0) + ss(k+1,i0) + ss(k+2,i0) - ccf(lag)=ccf(lag) + sum(ss(k,i0-m:i0+m)) + sum(ss(k+1,i0-m:i0+m)) & - + sum(ss(k+2,i0-m:i0+m)) - enddo - enddo - ccfmax=maxval(ccf) - ipk=maxloc(ccf) - lagbest=ipk(1)-1 - xdt=lagbest*tstep - 1.0 - - xsum=0. - sq=0. - nsum=0 - do i=0,lagmax - if(abs(i-lagbest).gt.2) then - xsum=xsum+ccf(i) - sq=sq+ccf(i)**2 - nsum=nsum+1 - endif - enddo - ave=xsum/nsum - rms=sqrt(sq/nsum - ave*ave) - snr=(ccfmax-ave)/rms - sync_ok=snr.ge.5.0 - - return -end subroutine q65_sync diff --git a/q65w/libq65/wideband_sync.f90 b/q65w/libq65/wideband_sync.f90 index 57ee103fc..6ef154f46 100644 --- a/q65w/libq65/wideband_sync.f90 +++ b/q65w/libq65/wideband_sync.f90 @@ -24,237 +24,4 @@ module wideband_sync type(sync_dat) :: sync(NFFT) integer nkhz_center - contains - -subroutine get_candidates(ss,savg,jz,nfa,nfb,nts_jt65,nts_q65,cand,ncand) - -! Search symbol spectra ss() over frequency range nfa to nfb (in kHz) for -! JT65 and Q65 sync patterns. The nts_* variables are the submode tone -! spacings: 1 2 4 8 16 for A B C D E. Birdies are detected and -! excised. Candidates are returned in the structure array cand(). - - parameter (MAX_PEAKS=100) - real ss(322,NFFT),savg(NFFT) - real pavg(-20:20) - integer indx(NFFT) - logical skip - type(candidate) :: cand(MAX_CANDIDATES) - type(candidate) :: cand0(MAX_CANDIDATES) - - call wb_sync(ss,savg,jz,nfa,nfb,nts_q65) !Output to sync() array - - tstep=2048.0/11025.0 !0.185760 s: 0.5*tsym_jt65, 0.3096*tsym_q65 - df3=96000.0/NFFT - ia=nint(1000*nfa/df3) + 1 - ib=nint(1000*nfb/df3) + 1 - if(ia.lt.1) ia=1 - if(ib.gt.NFFT-1) ib=NFFT-1 - iz=ib-ia+1 - - call indexx(sync(ia:ib)%ccfmax,iz,indx) !Sort by relative snr - - k=0 - do i=1,MAX_PEAKS - n=indx(iz+1-i) + ia - 1 - f0=0.001*(n-1)*df3 - snr1=sync(n)%ccfmax - if(snr1.lt.SNR1_THRESHOLD) exit - flip=sync(n)%iflip - if(flip.ne.0.0 .and. nts_jt65.eq.0) cycle - if(flip.eq.0.0 .and. nts_q65.eq.0) cycle - if(sync(n)%birdie) cycle - -! Test for signal outside of TxT range and set bw for this signal type - j1=(sync(n)%xdt + 1.0)/tstep - 1.0 - j2=(sync(n)%xdt + 52.0)/tstep + 1.0 - if(flip.ne.0) j2=(sync(n)%xdt + 47.811)/tstep + 1.0 - ipol=sync(n)%ipol - pavg=0. - do j=1,j1 - pavg=pavg + ss(j,n-20:n+20) - enddo - do j=j2,jz - pavg=pavg + ss(j,n-20:n+20) - enddo - jsum=j1 + (jz-j2+1) - pmax=maxval(pavg(-2:2)) !### Why not just pavg(0) ? - base=(sum(pavg)-pmax)/jsum - pmax=pmax/base - if(pmax.gt.5.0) cycle - skip=.false. - do m=1,k !Skip false syncs within signal bw - diffhz=1000.0*(f0-cand(m)%f) - bw=nts_q65*110.0 - if(cand(m)%iflip.ne.0) bw=nts_jt65*178.0 - if(diffhz.gt.-0.03*bw .and. diffhz.lt.1.03*bw) skip=.true. - enddo - if(skip) cycle - k=k+1 - cand(k)%snr=snr1 - cand(k)%f=f0 - cand(k)%xdt=sync(n)%xdt - cand(k)%pol=sync(n)%pol - cand(k)%ipol=sync(n)%ipol - cand(k)%iflip=nint(flip) - cand(k)%indx=n -! write(50,3050) i,k,m,f0+32.0,diffhz,bw,snr1,db(snr1) -!3050 format(3i5,f8.3,2f8.0,2f8.2) - if(k.ge.MAX_CANDIDATES) exit - enddo - ncand=k - - cand0(1:ncand)=cand(1:ncand) - call indexx(cand0(1:ncand)%f,ncand,indx) !Sort by relative snr - do i=1,ncand - k=indx(i) - cand(i)=cand0(k) - enddo - - return -end subroutine get_candidates - -subroutine wb_sync(ss,savg,jz,nfa,nfb,nts_q65) - -! Compute "orange sync curve" using the Q65 sync pattern - - use timer_module, only: timer - parameter (NFFT=32768) - parameter (LAGMAX=30) - real ss(322,NFFT) - real savg(NFFT) - real savg_med - logical first - integer isync(22) - integer jsync0(63),jsync1(63) - integer ip(1) - -! Q65 sync symbols - 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 jsync0/ & - 1, 4, 5, 9, 10, 11, 12, 13, 14, 16, 18, 22, 24, 25, 28, 32, & - 33, 34, 37, 38, 39, 40, 42, 43, 45, 46, 47, 48, 52, 53, 55, 57, & - 59, 60, 63, 64, 66, 68, 70, 73, 80, 81, 89, 90, 92, 95, 97, 98, & - 100,102,104,107,108,111,114,119,120,121,122,123,124,125,126/ - data jsync1/ & - 2, 3, 6, 7, 8, 15, 17, 19, 20, 21, 23, 26, 27, 29, 30, 31, & - 35, 36, 41, 44, 49, 50, 51, 54, 56, 58, 61, 62, 65, 67, 69, 71, & - 72, 74, 75, 76, 77, 78, 79, 82, 83, 84, 85, 86, 87, 88, 91, 93, & - 94, 96, 99,101,103,105,106,109,110,112,113,115,116,117,118/ - data first/.true./ - save first,isync,jsync0,jsync1 - - tstep=2048.0/11025.0 !0.185760 s: 0.5*tsym_jt65, 0.3096*tsym_q65 - if(first) then - fac=0.6/tstep - do i=1,22 !Expand the Q65 sync stride - isync(i)=nint((isync(i)-1)*fac) + 1 - enddo - do i=1,63 - jsync0(i)=2*(jsync0(i)-1) + 1 - jsync1(i)=2*(jsync1(i)-1) + 1 - enddo - first=.false. - endif - - df3=96000.0/NFFT - ia=nint(1000*nfa/df3) + 1 !Flat frequency range for WSE converters - ib=nint(1000*nfb/df3) + 1 - if(ia.lt.1) ia=1 - if(ib.gt.NFFT-1) ib=NFFT-1 - - call pctile(savg(ia:ib),ib-ia+1,50,savg_med) - - lagbest=0 - ipolbest=1 - flip=0. - - do i=ia,ib - ccfmax=0. - do lag=0,LAGMAX - - ccf=0. - ccf4=0. - do j=1,22 !Test for Q65 sync - k=isync(j) + lag - ccf4=ccf4 + ss(k,i+1) + ss(k+1,i+1) & - + ss(k+2,i+1) - enddo - ccf4=ccf4 - savg(i+1)*3*22/float(jz) - ccf=ccf4 - ipol=1 - if(ccf.gt.ccfmax) then - ipolbest=ipol - lagbest=lag - ccfmax=ccf - ccf4best=ccf4 - flip=0. - endif - - ccf=0. - ccf4=0. - do j=1,63 !Test for JT65 sync, std msg - k=jsync0(j) + lag - ccf4=ccf4 + ss(k,i+1) + ss(k+1,i+1) - enddo - ccf4=ccf4 - savg(i+1)*2*63/float(jz) - ccf=ccf4 - ipol=1 - if(ccf.gt.ccfmax) then - ipolbest=ipol - lagbest=lag - ccfmax=ccf - ccf4best=ccf4 - flip=1.0 - endif - - ccf=0. - ccf4=0. - do j=1,63 !Test for JT65 sync, OOO msg - k=jsync1(j) + lag - ccf4=ccf4 + ss(k,i+1) + ss(k+1,i+1) - enddo - ccf4=ccf4 - savg(i+1)*2*63/float(jz) - ccf=ccf4 - ipol=1 - if(ccf.gt.ccfmax) then - ipolbest=ipol - lagbest=lag - ccfmax=ccf - ccf4best=ccf4 - flip=-1.0 - endif - - enddo ! lag - - poldeg=0. - sync(i)%ccfmax=ccfmax - sync(i)%xdt=lagbest*tstep-1.0 - sync(i)%pol=poldeg - sync(i)%ipol=ipolbest - sync(i)%iflip=flip - sync(i)%birdie=.false. - if(ccfmax/(savg(i)/savg_med).lt.3.0) sync(i)%birdie=.true. - enddo ! i (frequency bin) - - call pctile(sync(ia:ib)%ccfmax,ib-ia+1,50,base) - sync(ia:ib)%ccfmax=sync(ia:ib)%ccfmax/base - - bw=65*nts_q65*1.66666667 !Q65-60x bandwidth - nbw=bw/df3 + 1 !Number of bins to blank - syncmin=2.0 - nguard=10 - do i=ia,ib - if(sync(i)%ccfmax.lt.syncmin) cycle - spk=maxval(sync(i:i+nbw)%ccfmax) - ip =maxloc(sync(i:i+nbw)%ccfmax) - i0=ip(1)+i-1 - ja=min(i,i0-nguard) - jb=i0+nbw+nguard - sync(ja:jb)%ccfmax=0. - sync(i0)%ccfmax=spk - enddo - - return -end subroutine wb_sync - end module wideband_sync