From 42a3e246bf2417c685622da2f680cff4d7fc06fa Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Fri, 18 Jun 2021 14:16:16 -0400 Subject: [PATCH] Compute pol angle from sync symbols. Display recommended TxPol. --- map65/libm65/CMakeLists.txt | 3 +- map65/libm65/map65a.f90 | 69 ++++------------- map65/libm65/polfit.f90 | 86 ++++++++++++++++++++++ map65/libm65/q65b.f90 | 44 +++++++---- map65/libm65/synctest.f90 | 2 +- map65/libm65/txpol.f90 | 33 +++++++++ map65/libm65/wideband_sync.f90 | 131 +++++++++++++++++++++------------ 7 files changed, 248 insertions(+), 120 deletions(-) create mode 100644 map65/libm65/polfit.f90 create mode 100644 map65/libm65/txpol.f90 diff --git a/map65/libm65/CMakeLists.txt b/map65/libm65/CMakeLists.txt index eff6dd9a8..238180322 100644 --- a/map65/libm65/CMakeLists.txt +++ b/map65/libm65/CMakeLists.txt @@ -60,8 +60,8 @@ set (libm65_FSRCS nchar.f90 noisegen.f90 packjt.f90 -# pctile.f90 pfxdump.f90 + polfit.f90 recvpkt.f90 rfile3a.f90 s3avg.f90 @@ -80,6 +80,7 @@ set (libm65_FSRCS trimlist.f90 twkfreq.f90 twkfreq_xy.f90 + txpol.f90 wavhdr.f90 f77_wisdom.f diff --git a/map65/libm65/map65a.f90 b/map65/libm65/map65a.f90 index 091d270d4..e8359a5a0 100644 --- a/map65/libm65/map65a.f90 +++ b/map65/libm65/map65a.f90 @@ -19,7 +19,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & real*8 fcenter character*22 msg(MAXMSG) character*3 shmsg0(4) - character mycall*12,hiscall*12,mygrid*6,hisgrid*6,grid*6,cp*1,cm*1 + character mycall*12,hiscall*12,mygrid*6,hisgrid*6,cp*1,cm*1 integer indx(MAXMSG),nsiz(MAXMSG) logical done(MAXMSG) logical xpol,bq65,q65b_called @@ -43,12 +43,13 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & mode65=mod(nmode,10) if(mode65.eq.3) mode65=4 mode_q65=nmode/10 + xpol=(nxpol.ne.0) nts_jt65=2**(mode65-1) !JT65 tone separation factor nts_q65=2**(mode_q65) !Q65 tone separation factor if(nagain.eq.0) then call timer('get_cand',0) - call get_candidates(ss,savg,mfa,mfb,nts_jt65,nts_q65,cand,ncand) + call get_candidates(ss,savg,xpol,mfa,mfb,nts_jt65,nts_q65,cand,ncand) call timer('get_cand',1) candec=.false. endif @@ -66,7 +67,6 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & mcall3a=mcall3b mousefqso0=mousefqso - xpol=(nxpol.ne.0) if(.not.xpol) ndphi=0 nsum=0 @@ -318,28 +318,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & if(npol.lt.0) npol=npol+180 endif -! If Tx station's grid is in decoded message, compute optimum TxPol - i1=index(decoded,' ') - i2=index(decoded(i1+1:),' ') + i1 - grid=' ' - if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm' - ntxpol=0 - cp=' ' - if(xpol) then - if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and. & - grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and. & - grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and. & - grid(4:4).ge.'0' .and. grid(4:4).le.'9') then - ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180) - if(nxant.eq.0) then - cp='H' - if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V' - else - cp='/' - if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\' - endif - endif - endif + call txpol(xpol,decoded,mygrid,npol,nxant,ntxpol,cp) if(ndphi.eq.0) then write(*,1010) nkHz,ndf,npol,nutc,dt,nsync2, & @@ -368,9 +347,9 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & q65b_called=.true. f0=cand(icand)%f call timer('q65b ',0) - call q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol, & - xpol,mycall,hiscall,hisgrid,mode_q65,f0,fqso,newdat, & - nagain,max_drift,idec) + call q65b(nutc,nqd,nxant,fcenter,nfcal,nfsample,ikhz,mousedf, & + ntol,xpol,mycall,mygrid, hiscall,hisgrid,mode_q65,f0,fqso, & + newdat,nagain,max_drift,idec) call timer('q65b ',1) if(idec.ge.0) candec(icand)=.true. enddo @@ -379,9 +358,9 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & ikhz=mousefqso f0=freq - (nkhz_center-48.0-1.27046) !### ??? ### call timer('q65b ',0) - call q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol, & - xpol,mycall,hiscall,hisgrid,mode_q65,f0,fqso,newdat, & - nagain,max_drift,idec) + call q65b(nutc,nqd,nxant,fcenter,nfcal,nfsample,ikhz,mousedf, & + ntol,xpol,mycall,mygrid,hiscall,hisgrid,mode_q65,f0,fqso, & + newdat,nagain,max_drift,idec) call timer('q65b ',1) endif endif @@ -416,8 +395,8 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & ikhz=nint(freq) f0=cand(icand)%f call timer('q65b ',0) - call q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol, & - xpol,mycall,hiscall,hisgrid,mode_q65,f0,fqso,newdat, & + call q65b(nutc,nqd,nxant,fcenter,nfcal,nfsample,ikhz,mousedf,ntol, & + xpol,mycall,mygrid,hiscall,hisgrid,mode_q65,f0,fqso,newdat, & nagain,max_drift,idec) call timer('q65b ',1) if(idec.ge.0) candec(icand)=.true. @@ -494,28 +473,8 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & if(npol.lt.0) npol=npol+180 endif -! If Tx station's grid is in decoded message, compute optimum TxPol - i1=index(decoded,' ') - i2=index(decoded(i1+1:),' ') + i1 - grid=' ' - if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm' - ntxpol=0 - cp=' ' - if(xpol) then - if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and. & - grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and. & - grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and. & - grid(4:4).ge.'0' .and. grid(4:4).le.'9') then - ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180) - if(nxant.eq.0) then - cp='H' - if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V' - else - cp='/' - if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\' - endif - endif - endif + call txpol(xpol,decoded,mygrid,npol,nxant,ntxpol,cp) + cmode='#A' if(mode65.eq.2) cmode='#B' if(mode65.eq.4) cmode='#C' diff --git a/map65/libm65/polfit.f90 b/map65/libm65/polfit.f90 new file mode 100644 index 000000000..7013432e2 --- /dev/null +++ b/map65/libm65/polfit.f90 @@ -0,0 +1,86 @@ +subroutine polfit(y,npts,a) + +! Input: y(npts) !Expect npts=4 +! Output: a(1) = baseline +! a(2) = amplitude +! a(3) = theta (deg) + + real y(npts) + real a(3) + real deltaa(3) + integer ipk(1) + save + +! Set starting values: + a(1)=minval(y) + a(2)=maxval(y)-a(1) + ipk=maxloc(y) + a(3)=(ipk(1)-1)*45.0 + + deltaa(1:2)=0.1*(a(2)-a(1)) + deltaa(3)=10.0 + nterms=3 + +! Start the iteration + chisqr=0. + chisqr0=1.e6 + iters=10 + + do iter=1,iters + do j=1,nterms + chisq1=fchisq_pol(y,npts,a) + fn=0. + delta=deltaa(j) +10 a(j)=a(j)+delta + chisq2=fchisq_pol(y,npts,a) + if(chisq2.eq.chisq1) go to 10 + if(chisq2.gt.chisq1) then + delta=-delta !Reverse direction + a(j)=a(j)+delta + tmp=chisq1 + chisq1=chisq2 + chisq2=tmp + endif +20 fn=fn+1.0 + a(j)=a(j)+delta + chisq3=fchisq_pol(y,npts,a) + if(chisq3.lt.chisq2) then + chisq1=chisq2 + chisq2=chisq3 + go to 20 + endif + +! Find minimum of parabola defined by last three points + delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5) + a(j)=a(j)-delta + deltaa(j)=deltaa(j)*fn/3. +! write(*,4000) iter,j,a,deltaa,chisq2 +!4000 format(2i2,2(2x,3f8.2),f12.5) + enddo + chisqr=fchisq_pol(y,npts,a) +! write(*,4000) 0,0,a,chisqr + if(deltaa(1).lt.0.01*(a(2)-a(1)) .and. deltaa(2).lt.0.01*(a(2)-a(1)) & + .and. deltaa(3).lt.1.0) exit + if(chisqr/chisqr0.gt.0.99) exit + a(3)=mod(a(3)+360.0,180.0) + chisqr0=chisqr + enddo + + return +end subroutine polfit + +real function fchisq_pol(y,npts,a) + + real y(npts),a(3) + data rad/57.2957795/ + + chisq = 0. + do i=1,npts + theta=(i-1)*45.0 + yfit=a(1) + a(2)*cos((theta-a(3))/rad)**2 + chisq=chisq + (y(i) - yfit)**2 + enddo + fchisq_pol=chisq + + return +end function fchisq_pol diff --git a/map65/libm65/q65b.f90 b/map65/libm65/q65b.f90 index c85e3be3d..1f1840b30 100644 --- a/map65/libm65/q65b.f90 +++ b/map65/libm65/q65b.f90 @@ -1,5 +1,6 @@ -subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, & - mycall0,hiscall0,hisgrid,mode_q65,f0,fqso,newdat,nagain,max_drift,idec) +subroutine q65b(nutc,nqd,nxant,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, & + mycall0,mygrid,hiscall0,hisgrid,mode_q65,f0,fqso,newdat,nagain, & + max_drift,idec) ! This routine provides an interface between MAP65 and the Q65 decoder ! in WSJT-X. All arguments are input data obtained from the MAP65 GUI. @@ -16,6 +17,7 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, & parameter (MAXFFT1=5376000) !56*96000 parameter (MAXFFT2=336000) !56*6000 (downsampled by 1/16) parameter (NMAX=60*12000) + parameter (RAD=57.2957795) ! type(hdr) h !Header for the .wav file integer*2 iwave(60*12000) complex ca(MAXFFT1),cb(MAXFFT1) !FFTs of raw x,y data @@ -25,10 +27,11 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, & real*8 fcenter,freq0 character*12 mycall0,hiscall0 character*12 mycall,hiscall - character*6 hisgrid + character*6 mygrid,hisgrid character*4 grid4 character*80 line character*80 wsjtx_dir + character*1 cp,cmode*2 common/cacb/ca,cb save @@ -88,10 +91,14 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, & ! 96000 5376000 0.017857143 336000 6000.000 ! 95238 5120000 0.018601172 322560 5999.994 - if(ipol.eq.1) cz(0:MAXFFT2-1)=cx - if(ipol.eq.2) cz(0:MAXFFT2-1)=0.707*(cx+cy) - if(ipol.eq.3) cz(0:MAXFFT2-1)=cy - if(ipol.eq.4) cz(0:MAXFFT2-1)=0.707*(cx-cy) + poldeg=0. + if(xpol) then + poldeg=sync(ipk)%pol + cz(0:MAXFFT2-1)=cos(poldeg/RAD)*cx + sin(poldeg/RAD)*cy + else + cz(0:MAXFFT2-1)=cx + endif + cz(MAXFFT2)=0. ! Roll off below 500 Hz and above 2500 Hz. ja=nint(500.0/df) @@ -136,21 +143,30 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, & if(nsnr0.gt.-99) then nq65df=nint(1000*(0.001*k0*df+nkhz_center-48.0+1.000-1.27046-ikhz))-nfcal nq65df=nq65df + nfreq0 - 1000 - + npol=nint(poldeg) + if(nxant.ne.0) then + npol=npol-45 + if(npol.lt.0) npol=npol+180 + endif + call txpol(xpol,msg0(1:22),mygrid,npol,nxant,ntxpol,cp) if(nqd.eq.1 .and. abs(nq65df-mousedf).lt.ntol) then - write(line,1020) ikhz,nq65df,45*(ipol-1),nutc,xdt0,nsnr0,msg0(1:27),cq0 -1020 format('!',i3.3,i5,i4,i6.4,f5.1,i5,' : ',a27,a3) + write(line,1020) ikhz,nq65df,npol,nutc,xdt0,nsnr0,msg0(1:27),cq0, & + ntxpol,cp +1020 format('!',i3.3,i5,i4,i6.4,f5.1,i5,' : ',a27,a3,i4,1x,a1) write(*,1100) trim(line) 1100 format(a) endif ! Write to lu 26, for Messages and Band Map windows - write(26,1014) freq0,nq65df,0,0,0,xdt0,45*(ipol-1),0, & - nsnr0,nutc,msg0(1:22),':',char(ichar('A') + mode_q65-1) -1014 format(f8.3,i5,3i3,f5.1,i4,i3,i4,i5.4,4x,a22,2x,a1,3x,':',a1) + + cmode=': ' + cmode(2:2)=char(ichar('A') + mode_q65-1) + write(26,1014) freq0,nq65df,0,0,0,xdt0,npol,0, & + nsnr0,nutc,msg0(1:22),':',cp,cmode +1014 format(f8.3,i5,3i3,f5.1,i4,i3,i4,i5.4,4x,a22,1x,2a1,2x,a2) ! Write to file map65_rx.log: - write(21,1110) freq0,nq65df,xdt0,45*(ipol-1),nsnr0,nutc,msg0(1:28),cq0 + write(21,1110) freq0,nq65df,xdt0,npol,nsnr0,nutc,msg0(1:28),cq0 1110 format(f8.3,i5,f5.1,2i4,i5.4,2x,a28,': A',2x,a3) endif diff --git a/map65/libm65/synctest.f90 b/map65/libm65/synctest.f90 index f84042eae..a4a851c5b 100644 --- a/map65/libm65/synctest.f90 +++ b/map65/libm65/synctest.f90 @@ -41,7 +41,7 @@ program synctest call timer('synctest',0) call timer('get_cand',0) - call get_candidates(ss,savg,nfa,nfb,nts_jt65,nts_q65,cand,ncand) + call get_candidates(ss,savg,.true.,nfa,nfb,nts_jt65,nts_q65,cand,ncand) call timer('get_cand',1) do k=1,ncand diff --git a/map65/libm65/txpol.f90 b/map65/libm65/txpol.f90 new file mode 100644 index 000000000..042d3e549 --- /dev/null +++ b/map65/libm65/txpol.f90 @@ -0,0 +1,33 @@ +subroutine txpol(xpol,decoded,mygrid,npol,nxant,ntxpol,cp) + +! If Tx station's grid is in decoded message, compute optimum TxPol + character*22 decoded + character*6 mygrid,grid + character*1 cp + logical xpol + + ntxpol=0 + i1=index(decoded,' ') + i2=index(decoded(i1+1:),' ') + i1 + grid=' ' + if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm' + ntxpol=0 + cp=' ' + if(xpol .and.grid(1:4).ne.'RR73') then + if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and. & + grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and. & + grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and. & + grid(4:4).ge.'0' .and. grid(4:4).le.'9') then + ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180) + if(nxant.eq.0) then + cp='H' + if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V' + else + cp='/' + if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\' + endif + endif + endif + + return +end subroutine txpol diff --git a/map65/libm65/wideband_sync.f90 b/map65/libm65/wideband_sync.f90 index 78fa9812f..352d985fc 100644 --- a/map65/libm65/wideband_sync.f90 +++ b/map65/libm65/wideband_sync.f90 @@ -4,12 +4,14 @@ module wideband_sync real :: snr !Relative S/N of sync detection real :: f !Freq of sync tone, 0 to 96000 Hz real :: xdt !DT of matching sync pattern, -1.0 to +4.0 s + real :: pol !Polarization angle, degrees integer :: ipol !Polarization angle, 1 to 4 ==> 0, 45, 90, 135 deg integer :: iflip !Sync type: JT65 = +/- 1, Q65 = 0 end type candidate type sync_dat real :: ccfmax real :: xdt + real :: pol integer :: ipol integer :: iflip logical :: birdie @@ -17,12 +19,13 @@ module wideband_sync parameter (NFFT=32768) parameter (MAX_CANDIDATES=50) + parameter (SNR1_THRESHOLD=4.5) type(sync_dat) :: sync(NFFT) integer nkhz_center contains -subroutine get_candidates(ss,savg,nfa,nfb,nts_jt65,nts_q65,cand,ncand) +subroutine get_candidates(ss,savg,xpol,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 @@ -33,7 +36,7 @@ subroutine get_candidates(ss,savg,nfa,nfb,nts_jt65,nts_q65,cand,ncand) real ss(4,322,NFFT),savg(4,NFFT) real pavg(-20:20) integer indx(NFFT) - logical skip + logical xpol,skip type(candidate) :: cand(MAX_CANDIDATES) do j=322,1,-1 !Find end of data in ss() @@ -41,7 +44,7 @@ subroutine get_candidates(ss,savg,nfa,nfb,nts_jt65,nts_q65,cand,ncand) enddo jz=j -call wb_sync(ss,savg,jz,nfa,nfb) +call wb_sync(ss,savg,xpol,jz,nfa,nfb) tstep=2048.0/11025.0 !0.185760 s: 0.5*tsym_jt65, 0.3096*tsym_q65 df3=96000.0/NFFT @@ -57,7 +60,7 @@ call wb_sync(ss,savg,jz,nfa,nfb) f0=0.001*(n-1)*df3 snr1=sync(n)%ccfmax ! print*,'=A',f0,snr1 - if(snr1.lt.4.5) exit + 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 @@ -94,6 +97,7 @@ call wb_sync(ss,savg,jz,nfa,nfb) 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) if(k.ge.MAX_CANDIDATES) exit @@ -103,18 +107,21 @@ call wb_sync(ss,savg,jz,nfa,nfb) return end subroutine get_candidates -subroutine wb_sync(ss,savg,jz,nfa,nfb) +subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) ! Compute "orange sync curve" using the Q65 sync pattern + use timer_module, only: timer parameter (NFFT=32768) parameter (LAGMAX=30) real ss(4,322,NFFT) real savg(4,NFFT) real savg_med(4) - logical first + real ccf4(4),ccf4best(4),a(3) + logical first,xpol 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/ @@ -147,12 +154,14 @@ subroutine wb_sync(ss,savg,jz,nfa,nfb) df3=96000.0/NFFT ia=nint(1000*nfa/df3) + 1 !Flat frequency range for WSE converters ib=nint(1000*nfb/df3) + 1 + npol=1 + if(xpol) npol=4 - do i=1,4 + do i=1,npol call pctile(savg(i,ia:ib),ib-ia+1,50,savg_med(i)) enddo ! do i=ia,ib -! write(14,3014) 0.001*(i-1)*df3,savg(1:4,i) +! write(14,3014) 0.001*(i-1)*df3,savg(1:npol,i) !3014 format(5f10.3) ! enddo @@ -162,57 +171,82 @@ subroutine wb_sync(ss,savg,jz,nfa,nfb) do i=ia,ib ccfmax=0. - do ipol=1,4 - do lag=0,LAGMAX + do lag=0,LAGMAX - ccf=0. - do j=1,22 !Test for Q65 sync - k=isync(j) + lag - ccf=ccf + ss(ipol,k,i+1) + ss(ipol,k+1,i+1) + ss(ipol,k+2,i+1) - enddo - ccf=ccf - savg(ipol,i+1)*3*22/float(jz) - if(ccf.gt.ccfmax) then - ipolbest=ipol - lagbest=lag - ccfmax=ccf - flip=0. - endif + ccf=0. + ccf4=0. + do j=1,22 !Test for Q65 sync + k=isync(j) + lag + ccf4=ccf4 + ss(1:npol,k,i+1) + ss(1:npol,k+1,i+1) + ss(1:npol,k+2,i+1) + enddo + ccf4=ccf4 - savg(1:npol,i+1)*3*22/float(jz) + ccf=maxval(ccf4) + ip=maxloc(ccf4) + ipol=ip(1) + if(ccf.gt.ccfmax) then + ipolbest=ipol + lagbest=lag + ccfmax=ccf + ccf4best=ccf4 + flip=0. + endif - ccf=0. - do j=1,63 !Test for JT65 sync, std msg - k=jsync0(j) + lag - ccf=ccf + ss(ipol,k,i+1) + ss(ipol,k+1,i+1) - enddo - ccf=ccf - savg(ipol,i+1)*2*63/float(jz) - if(ccf.gt.ccfmax) then - ipolbest=ipol - lagbest=lag - ccfmax=ccf - flip=1.0 - endif + ccf=0. + ccf4=0. + do j=1,63 !Test for JT65 sync, std msg + k=jsync0(j) + lag + ccf4=ccf4 + ss(1:npol,k,i+1) + ss(1:npol,k+1,i+1) + enddo + ccf4=ccf4 - savg(1:npol,i+1)*2*63/float(jz) + ccf=maxval(ccf4) + ip=maxloc(ccf4) + ipol=ip(1) + if(ccf.gt.ccfmax) then + ipolbest=ipol + lagbest=lag + ccfmax=ccf + ccf4best=ccf4 + flip=1.0 + endif - ccf=0. - do j=1,63 !Test for JT65 sync, OOO msg - k=jsync1(j) + lag - ccf=ccf + ss(ipol,k,i+1) + ss(ipol,k+1,i+1) - enddo - ccf=ccf - savg(ipol,i+1)*2*63/float(jz) - if(ccf.gt.ccfmax) then - ipolbest=ipol - lagbest=lag - ccfmax=ccf - 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(1:npol,k,i+1) + ss(1:npol,k+1,i+1) + enddo + ccf4=ccf4 - savg(1:npol,i+1)*2*63/float(jz) + ccf=maxval(ccf4) + ip=maxloc(ccf4) + ipol=ip(1) + if(ccf.gt.ccfmax) then + ipolbest=ipol + lagbest=lag + ccfmax=ccf + ccf4best=ccf4 + flip=-1.0 + endif - enddo ! lag - enddo !ipol + enddo ! lag + poldeg=0. + if(xpol .and. ccfmax.ge.SNR1_THRESHOLD) then + call polfit(ccf4best,4,a) + poldeg=a(3) + endif 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(ipolbest,i)/savg_med(ipolbest)).lt.3.0) sync(i)%birdie=.true. + if(sync(i)%iflip.eq.0 .and. sync(i)%ccfmax .gt. 20.0) then + write(50,3050) i,lagbest,sync(i)%ccfmax,sync(i)%xdt,sync(i)%ipol, & + sync(i)%birdie,ccf4best +3050 format(2i5,f10.3,f8.2,i5,1x,L3,4f7.1) + endif + enddo ! i (frequency bin) ! do i=ia,ib @@ -223,7 +257,6 @@ subroutine wb_sync(ss,savg,jz,nfa,nfb) call pctile(sync(ia:ib)%ccfmax,ib-ia+1,50,base) sync(ia:ib)%ccfmax=sync(ia:ib)%ccfmax/base -! print*,base return end subroutine wb_sync