From 15ddaf9e8a267153fb15a9eaaf9229df95ac0040 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 17 Mar 2016 18:52:06 +0000 Subject: [PATCH] Various changes to JT65 decoding, all potentially temporary. 1. Measure Doppler width by fitting a (modified) Lorentzian. 2. Don't call "slope" in sync65(). 3. New definition of "sync1". 4. Get snr from sync1. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6536 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- CMakeLists.txt | 2 + lib/fchisq0.f90 | 23 ++++++++++ lib/fer65.f90 | 5 ++- lib/jt65_decode.f90 | 54 ++++++++--------------- lib/jt65_test.f90 | 4 +- lib/lorentzian.f90 | 105 ++++++++++++++++++++++++++++++++++++++++++++ lib/slope.f90 | 2 +- lib/sync65.f90 | 22 ++++++++-- 8 files changed, 175 insertions(+), 42 deletions(-) create mode 100644 lib/fchisq0.f90 create mode 100644 lib/lorentzian.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index f68efa717..03dc1a977 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -309,6 +309,7 @@ set (wsjt_FSRCS lib/fast9.f90 lib/fast_decode.f90 lib/fchisq.f90 + lib/fchisq0.f90 lib/fchisq65.f90 lib/fftw3mod.f90 lib/fil3.f90 @@ -362,6 +363,7 @@ set (wsjt_FSRCS lib/jtmsk_decode.f90 lib/jtmsk_short.f90 lib/libration.f90 + lib/lorentzian.f90 lib/lpf1.f90 lib/mixlpf.f90 lib/moondopjpl.f90 diff --git a/lib/fchisq0.f90 b/lib/fchisq0.f90 new file mode 100644 index 000000000..e01872542 --- /dev/null +++ b/lib/fchisq0.f90 @@ -0,0 +1,23 @@ +real function fchisq0(y,npts,a) + + real y(npts),a(4) + + rewind 51 + chisq = 0. + do i=1,npts + x=i + z=(x-a(3))/(0.5*a(4)) + yfit=a(1) + if(abs(z).lt.3.0) then + d=1.0 + z*z + yfit=a(1) + a(2) * (1.0/d - 0.1) + endif + chisq=chisq + (y(i) - yfit)**2 +! write(51,3001) i,y(i),yfit,y(i)-yfit +!3001 format(i5,3f10.4) + enddo + fchisq0=chisq + + return +end function fchisq0 + diff --git a/lib/fer65.f90 b/lib/fer65.f90 index 6d6604a1a..e56c75829 100644 --- a/lib/fer65.f90 +++ b/lib/fer65.f90 @@ -151,7 +151,10 @@ program fer65 dsnr=xsnr-snr dfreq=xfreq-1500.0 - if(ngood.eq.0) dsnr=0. + if(ngood.eq.0) then + dsnr=0. + dfreq=0. + endif write(20,1100) snr,nsync,ngood,nbad,xsync,esync,dsnr,esnr, & xdt,edt,dfreq,efreq,xdrift,edrift,xwidth,ewidth 1100 format(f5.1,2i6i4,2f6.1,f6.1,f5.1,f6.2,f5.2,6f5.1) diff --git a/lib/jt65_decode.f90 b/lib/jt65_decode.f90 index ce5b82ab9..6d996ba15 100644 --- a/lib/jt65_decode.f90 +++ b/lib/jt65_decode.f90 @@ -124,6 +124,19 @@ contains nfb=min(4000,nfqso+ntol) thresh0=1.0 endif + df=12000.0/8192.0 !df = 1.465 Hz + if(single_decode) then + ia=max(1,nint(nfa/df)-100) + ib=min(NSZ,nint(nfb/df)+100) + nz=ib-ia+1 + call lorentzian(savg(ia),nz,a) + baseline=a(1) + amp=a(2) + f0=(a(3)+ia-1)*df + width=a(4)*df +! write(*,3001) baseline,amp,f0,width +!3001 format(4f10.3) + endif ! robust = .false.: use float ccf. Only if ncand>50 fall back to robust (1-bit) ccf ! robust = .true. : use only robust (1-bit) ccf @@ -143,39 +156,7 @@ contains ! If a candidate was found within +/- ntol of nfqso, move it into ca(1). call fqso_first(nfqso,ntol,ca,ncand) - df=12000.0/8192.0 !df = 1.465 Hz - width=0. - if(single_decode) then - ncand=1 - smax=-1.e30 - do i=151,NSZ-150 - if(savg(i).gt.smax) then - smax=savg(i) - ipk=i - endif -! write(50,3001) i*df,savg(i) -!3001 format(2f12.3) - enddo - base=(sum(savg(ipk-149:ipk-50)) + sum(savg(ipk+51:ipk+150)))/200.0 - - stest=smax - 0.5*(smax-base) - ssum=savg(ipk) - do i=1,50 - if(savg(ipk+i).lt.stest) exit - ssum=ssum + savg(ipk+i) - enddo - do i=1,50 - if(savg(ipk-i).lt.stest) exit - ssum=ssum + savg(ipk-i) - enddo - ww=ssum/savg(ipk) - width=2 - t=ww*ww - 5.67 - if(t.gt.0.0) width=sqrt(t) - width=df*width -! print*,'Width:',width - endif - + if(single_decode) ncand=1 nvec=ntrials if(ncand.gt.75) then ! write(*,*) 'Pass ',ipass,' ncandidates too large ',ncand @@ -220,9 +201,10 @@ contains nfreq=nint(freq+a(1)) ndrift=nint(2.0*a(2)) - s2db=10.0*log10(sync2) - 35 !### empirical ### - if(width.gt.3) s2db=s2db + 2.1*sqrt(width-3.0) + 1.5 + & - 0.11*(width-7.0) !### empirical^2 ### +! s2db=10.0*log10(sync2) - 35 !### empirical ### +! if(width.gt.3) s2db=s2db + 2.1*sqrt(width-3.0) + 1.5 + & +! 0.11*(width-7.0) !### empirical^2 ### + s2db=sync1 - 30.0 nsnr=nint(s2db) if(nsnr.lt.-30) nsnr=-30 if(nsnr.gt.-1) nsnr=-1 diff --git a/lib/jt65_test.f90 b/lib/jt65_test.f90 index 40fbc43ba..cb5b93e31 100644 --- a/lib/jt65_test.f90 +++ b/lib/jt65_test.f90 @@ -58,8 +58,10 @@ contains integer, intent(in) :: submode integer, intent(in) :: aggression integer nwidth + real t - nwidth=max(nint(width),2) + t=max(0.0,width*width-7.2) + nwidth=max(nint(sqrt(t)),2) write(*,1010) utc,snr,dt,freq,decoded 1010 format(i4.4,i4,f5.1,i5,1x,'#',1x,a22) write(13,1012) utc,nint(sync),snr,dt,freq,drift,nwidth, & diff --git a/lib/lorentzian.f90 b/lib/lorentzian.f90 new file mode 100644 index 000000000..4ea349d7b --- /dev/null +++ b/lib/lorentzian.f90 @@ -0,0 +1,105 @@ +subroutine lorentzian(y,npts,a) + +! Input: y(npts); assume x(i)=i, i=1,npts +! Output: a(4) +! a(1) = baseline +! a(2) = amplitude +! a(3) = x0 +! a(4) = width + + real y(npts) + real a(4) + real deltaa(4) + real x(1000) + + if(npts.gt.1000) stop 'Error in lorentzian' + + a=0. + df=12000.0/8192.0 !df = 1.465 Hz + width=0. + ymax=-1.e30 + do i=1,npts + if(y(i).gt.ymax) then + ymax=y(i) + ipk=i + endif +! write(50,3001) i,i*df,y(i) +!3001 format(i6,2f12.3) + enddo +! base=(sum(y(ipk-149:ipk-50)) + sum(y(ipk+51:ipk+150)))/200.0 + base=(sum(y(1:20)) + sum(y(npts-19:npts)))/40.0 + stest=ymax - 0.5*(ymax-base) + ssum=y(ipk) + do i=1,50 + if(ipk+i.gt.npts) exit + if(y(ipk+i).lt.stest) exit + ssum=ssum + y(ipk+i) + enddo + do i=1,50 + if(ipk-i.lt.1) exit + if(y(ipk-i).lt.stest) exit + ssum=ssum + y(ipk-i) + enddo + ww=ssum/y(ipk) + width=2 + t=ww*ww - 5.67 + if(t.gt.0.0) width=sqrt(t) + a(1)=base + a(2)=ymax-base + a(3)=ipk + a(4)=width + +! Now find Lorentzian parameters + + deltaa(1)=0.1 + deltaa(2)=0.1 + deltaa(3)=1.0 + deltaa(4)=1.0 + nterms=4 + +! Start the iteration + chisqr=0. + chisqr0=1.e6 + do iter=1,5 + do j=1,nterms + chisq1=fchisq0(y,npts,a) + fn=0. + delta=deltaa(j) +10 a(j)=a(j)+delta + chisq2=fchisq0(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=fchisq0(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,chisq2 +!4000 format(i1,i2,4f10.4,f11.3) + enddo + chisqr=fchisq0(y,npts,a) +! write(*,4000) 0,0,a,chisqr + if(chisqr/chisqr0.gt.0.99) go to 30 + chisqr0=chisqr + enddo + +30 ccfbest=ccfmax * (1378.125/fsample)**2 + dtbest=dtmax + + return +end subroutine lorentzian + diff --git a/lib/slope.f90 b/lib/slope.f90 index f39463217..df85242c5 100644 --- a/lib/slope.f90 +++ b/lib/slope.f90 @@ -1,7 +1,7 @@ subroutine slope(y,npts,xpk) ! Remove best-fit slope from data in y(i). When fitting the straight line, -! ignore the peak around xpk +/- 2. +! ignore the peak around xpk +/- nhw real y(npts) diff --git a/lib/sync65.f90 b/lib/sync65.f90 index 67645f40c..9b08066da 100644 --- a/lib/sync65.f90 +++ b/lib/sync65.f90 @@ -31,7 +31,7 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust) do i=ia,ib call xcor(ss,i,nhsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk0,flip,fdot,nrobust) ! Remove best-fit slope from ccfblue and normalize so baseline rms=1.0 - call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0) +! call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0) ccfred(i)=ccfblue(lagpk0) if(ccfred(i).gt.ccfmax) then ccfmax=ccfred(i) @@ -43,6 +43,12 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust) ccfred(ia-1)=ccfred(ia) ccfred(ib+1)=ccfred(ib) +! do i=1,NSZ +! ssum=sum(ss(1:322,i))/322.0 +! write(61,3001) i*df,ccfred(i),ssum +! enddo + + do i=ia,ib freq=i*df itry=0 @@ -59,7 +65,7 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust) endif if(itry.ne.0) then call xcor(ss,i,nhsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot,nrobust) - call slope(ccfblue(lag1),lag2-lag1+1,lagpk-lag1+1.0) +! call slope(ccfblue(lag1),lag2-lag1+1,lagpk-lag1+1.0) xlag=lagpk if(lagpk.gt.lag1 .and. lagpk.lt.lag2) then call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2) @@ -70,10 +76,20 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust) ccfblue(lag2)=0. ca(ncand)%freq=freq ca(ncand)%dt=dtx - ca(ncand)%sync=ccfred(i) +! ca(ncand)%sync=ccfred(i) + ca(ncand)%sync=db(ccfred(i)) - 16.0 endif if(ncand.eq.MAXCAND) return enddo +! do i=1,NSZ +! write(52,3001) i*df,ccfred(i) +!3001 format(3f12.3) +! enddo +! do i=-10,58 +! write(53,3002) i,i*2048.0/11025.0,ccfblue(i) +!3002 format(i6,2f12.3) +! enddo + return end subroutine sync65