From 9f388b63dcefb56c17318720283a7d86d6ae7516 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 17 Mar 2016 13:28:57 +0000 Subject: [PATCH] Many changes in aid of decoding signals with significant Doppler spread in submodes JT65B, C. More changes still to come! git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6535 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/decode65a.f90 | 11 ++++++--- lib/decoder.f90 | 7 +++--- lib/fer65.f90 | 54 +++++++++++++++++++++++++++------------------ lib/jt65.f90 | 2 +- lib/jt65_decode.f90 | 52 +++++++++++++++++++++++++++++++++++-------- lib/jt65_mod.f90 | 1 + lib/jt65_test.f90 | 12 ++++++---- lib/jt65sim.f90 | 2 +- 8 files changed, 99 insertions(+), 42 deletions(-) diff --git a/lib/decode65a.f90 b/lib/decode65a.f90 index 60aef1242..714966e4e 100644 --- a/lib/decode65a.f90 +++ b/lib/decode65a.f90 @@ -83,11 +83,14 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials, & call timer('dec65b ',0) qualbest=0. + minsmo=0 maxsmo=0 - if(mode65.eq.2) maxsmo=5 - if(mode65.eq.4) maxsmo=10 + if(mode65.ge.2) then + minsmo=nint(width/df) + maxsmo=2*minsmo + endif nn=0 - do ismo=0,maxsmo + do ismo=minsmo,maxsmo if(ismo.gt.0) then do j=1,126 call smo121(s1(-255,j),512) @@ -127,6 +130,8 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials, & endif enddo +! print*,width,minsmo,maxsmo,nsmo,nn + if(nft.eq.2) then decoded=decoded_best qual=qualbest diff --git a/lib/decoder.f90 b/lib/decoder.f90 index ef192511c..21b11e2f9 100644 --- a/lib/decoder.f90 +++ b/lib/decoder.f90 @@ -198,7 +198,7 @@ contains 1000 format(a1,i5.4,f6.1,f6.2,i6,1x,a1) end subroutine jt4_average - subroutine jt65_decoded(this,utc,sync,snr,dt,freq,drift,decoded,ft, & + subroutine jt65_decoded(this,utc,sync,snr,dt,freq,drift,width,decoded,ft, & qual,nsmo,nsum,minsync,nsubmode,naggressive) use jt65_decode @@ -211,6 +211,7 @@ contains real, intent(in) :: dt integer, intent(in) :: freq integer, intent(in) :: drift + real, intent(in) :: width character(len=22), intent(in) :: decoded integer, intent(in) :: ft integer, intent(in) :: qual @@ -243,8 +244,8 @@ contains 1010 format(i4.4,i4,f5.1,i5,1x,a1,1x,a22,a3) endif - write(13,1012) utc,nint(sync),snr,dt,float(freq),drift,decoded,ft,nsmo -1012 format(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22,' JT65',i4,i2) + write(13,1012) utc,nint(sync),snr,dt,float(freq),drift,decoded,ft,nsum,nsmo +1012 format(i4.4,i4,i5,f6.2,f8.0,i4,3x,a22,' JT65',3i3) call flush(6) !$omp end critical(decode_results) diff --git a/lib/fer65.f90 b/lib/fer65.f90 index 62d36838a..6d6604a1a 100644 --- a/lib/fer65.f90 +++ b/lib/fer65.f90 @@ -18,7 +18,7 @@ program fer65 ! -s S/N in 2500 Hz -s single-decode mode implicit real*8 (a-h,o-z) - real*8 s(5),sq(5) + real*8 s(6),sq(6) character arg*12,cmnd*100,decoded*22,submode*1,csync*1 logical syncok @@ -51,11 +51,11 @@ program fer65 open(21,file='fer65.21',status='unknown') write(20,1000) submode,iters,ntrials,naggressive,d -1000 format('JT65',a1,' Iters:',i6,' T:',i7,' Aggressive:',i3, & +1000 format(/'JT65',a1,' Iters:',i6,' T:',i7,' Aggressive:',i3, & ' Doppler:',f6.1) write(20,1002) 1002 format(/' dB nsync ngood nbad sync dsnr ', & - 'DT Freq Drift'/77('-')) + 'DT Freq Drift Width'/85('-')) do isnr=0,20 snr=snr1+isnr @@ -77,37 +77,41 @@ program fer65 isync=0 nsnr=0 dt=0. - freq=0. + nfreq=0 ndrift=0 + nwidth=0 cmnd='./jt65 -m A -a 10 -f 1500 -n 1000 -d 3 -s -X 32 000000_0001.wav > decoded.txt' cmnd(11:11)=submode ! print*,cmnd call system(cmnd) open(13,file='fort.13',status='old',err=20) - read(13,1012) nutc,isync,nsnr,dt,freq,ndrift,decoded,nft,nsum,nsmo -1012 format(i4,i4,i5,f6.2,f8.0,i4,3x,a22,5x,3i3) + read(13,1012) nutc,isync,nsnr,dt,nfreq,ndrift,nwidth,decoded, & + nft,nsum,nsmo +1012 format(i4,i4,i5,f6.2,i5,i4,i3,1x,a22,5x,3i3) close(13) - syncok=abs(dt).lt.0.2 .and. abs(freq-1500.0).lt.dfmax + syncok=abs(dt).lt.0.2 .and. float(abs(nfreq-1500)).lt.dfmax csync=' ' if(syncok) csync='*' - write(21,1014) nutc,isync,nsnr,dt,freq,ndrift,nft,nsum,nsmo,csync, & - decoded(1:16) -1014 format(i4,i4,i5,f6.2,f8.0,i4,3x,3i3,1x,a1,1x,a16) + write(21,1014) nutc,isync,nsnr,dt,nfreq,ndrift,nwidth, & + nft,nsum,nsmo,csync,decoded(1:16) +1014 format(i4,i4,i5,f6.2,i5,i4,3x,4i3,1x,a1,1x,a16) if(syncok) then nsync=nsync+1 + s(1)=s(1) + isync + sq(1)=sq(1) + isync*isync + s(6)=s(6) + nwidth + sq(6)=sq(6) + nwidth*nwidth if(decoded.eq.'K1ABC W9XYZ EN37 ') then ngood=ngood+1 - s(1)=s(1) + isync s(2)=s(2) + nsnr s(3)=s(3) + dt - s(4)=s(4) + freq + s(4)=s(4) + nfreq s(5)=s(5) + ndrift - sq(1)=sq(1) + isync*isync sq(2)=sq(2) + nsnr*nsnr sq(3)=sq(3) + dt*dt - sq(4)=sq(4) + freq*freq + sq(4)=sq(4) + nfreq*nfreq sq(5)=sq(5) + ndrift*ndrift else if(decoded.ne.' ') then nbad=nbad+1 @@ -118,20 +122,27 @@ program fer65 fsync=float(nsync)/iter fgood=float(ngood)/iter fbad=float(nbad)/iter - write(*,1020) iter,isync,nsnr,dt,int(freq),ndrift,fsync,fgood,fbad, & - decoded(1:16) -1020 format(i8,2i4,f7.2,i6,i4,2f7.3,f10.6,1x,a16) + write(*,1020) iter,isync,nsnr,dt,nfreq,ndrift,nwidth,fsync,fgood, & + fbad,decoded(1:18) +1020 format(i8,2i4,f7.2,i6,i4,i3,2f7.3,f8.4,1x,a18) enddo + if(nsync.ge.1) then + xsync=s(1)/nsync + xwidth=s(6)/nsync + endif + if(nsync.ge.2) then + esync=sqrt(sq(1)/nsync - xsync**2) + ewidth=sqrt(sq(6)/nsync - xwidth**2) + endif + if(ngood.ge.1) then - xsync=s(1)/ngood xsnr=s(2)/ngood xdt=s(3)/ngood xfreq=s(4)/ngood xdrift=s(5)/ngood endif if(ngood.ge.2) then - esync=sqrt(sq(1)/ngood - xsync**2) esnr=sqrt(sq(2)/ngood - xsnr**2) edt=sqrt(sq(3)/ngood - xdt**2) efreq=sqrt(sq(4)/ngood - xfreq**2) @@ -139,10 +150,11 @@ program fer65 endif dsnr=xsnr-snr + dfreq=xfreq-1500.0 if(ngood.eq.0) dsnr=0. write(20,1100) snr,nsync,ngood,nbad,xsync,esync,dsnr,esnr, & - xdt,edt,xfreq,efreq,xdrift,edrift -1100 format(f5.1,2i6i4,2f6.1,f6.1,f5.1,f6.2,f5.2,f7.1,3f5.1) + 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) flush(20) enddo diff --git a/lib/jt65.f90 b/lib/jt65.f90 index 71c720877..d3c51cba5 100644 --- a/lib/jt65.f90 +++ b/lib/jt65.f90 @@ -79,7 +79,7 @@ program jt65 case ('X') read (optarg(:narglen), *) nexp_decoded case ('s') - ntol=10 + ntol=100 nlow=nfqso-ntol nhigh=nfqso+ntol n2pass=1 diff --git a/lib/jt65_decode.f90 b/lib/jt65_decode.f90 index 484b6d108..ce5b82ab9 100644 --- a/lib/jt65_decode.f90 +++ b/lib/jt65_decode.f90 @@ -1,6 +1,6 @@ module jt65_decode - integer, parameter :: NSZ=3413, NZMAX=60*12000, NFFT=1000 + integer, parameter :: NSZ=3413, NZMAX=60*12000 type :: jt65_decoder procedure(jt65_decode_callback), pointer :: callback => null() @@ -13,7 +13,7 @@ module jt65_decode ! abstract interface subroutine jt65_decode_callback(this,utc,sync,snr,dt,freq,drift, & - decoded,ft,qual,nsmo,nsum,minsync,nsubmode,naggressive) + width,decoded,ft,qual,nsmo,nsum,minsync,nsubmode,naggressive) import jt65_decoder implicit none @@ -24,6 +24,7 @@ module jt65_decode real, intent(in) :: dt integer, intent(in) :: freq integer, intent(in) :: drift + real, intent(in) :: width character(len=22), intent(in) :: decoded integer, intent(in) :: ft integer, intent(in) :: qual @@ -142,7 +143,38 @@ contains ! If a candidate was found within +/- ntol of nfqso, move it into ca(1). call fqso_first(nfqso,ntol,ca,ncand) - if(single_decode) ncand=1 + 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 nvec=ntrials if(ncand.gt.75) then @@ -150,7 +182,6 @@ contains nvec=100 endif - df=12000.0/NFFT !df = 12000.0/8192 = 1.465 Hz mode65=2**nsubmode nflip=1 !### temporary ### nqd=0 @@ -164,9 +195,9 @@ contains endif do icand=1,ncand - freq=ca(icand)%freq - dtx=ca(icand)%dt sync1=ca(icand)%sync + dtx=ca(icand)%dt + freq=ca(icand)%freq if(ipass.eq.1) ntry65a=ntry65a + 1 if(ipass.eq.2) ntry65b=ntry65b + 1 call timer('decod65a',0) @@ -190,6 +221,8 @@ 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 ### nsnr=nint(s2db) if(nsnr.lt.-30) nsnr=-30 if(nsnr.gt.-1) nsnr=-1 @@ -208,7 +241,7 @@ contains if (associated(this%callback) .and. nsum.ge.2) then call this%callback(nutc,sync1,nsnr,dtx-1.0,nfreq,ndrift, & - avemsg,nftt,nqual,nsmo,nsum,minsync,nsubmode, & + width,avemsg,nftt,nqual,nsmo,nsum,minsync,nsubmode, & naggressive) prtavg=.true. cycle @@ -231,7 +264,8 @@ contains if(rtt.gt.r0(n)) cycle endif -5 if(decoded.eq.decoded0 .and. abs(freq-freq0).lt. 3.0 .and. & +5 continue + if(decoded.eq.decoded0 .and. abs(freq-freq0).lt. 3.0 .and. & minsync.ge.0) cycle !Don't display dupes if(decoded.ne.' ' .or. minsync.lt.0) then if( nsubtract .eq. 1 ) then @@ -259,7 +293,7 @@ contains nqual=min(qual,9999.0) if (associated(this%callback)) then call this%callback(nutc,sync1,nsnr,dtx-1.0,nfreq,ndrift, & - decoded,nft,nqual,nsmo,nsum,minsync,nsubmode, & + width,decoded,nft,nqual,nsmo,nsum,minsync,nsubmode, & naggressive) end if endif diff --git a/lib/jt65_mod.f90 b/lib/jt65_mod.f90 index 8a7f01a3a..9bf9e2617 100644 --- a/lib/jt65_mod.f90 +++ b/lib/jt65_mod.f90 @@ -8,5 +8,6 @@ module jt65_mod real s1(-255:256,126) real s3a(64,63) real pr(126) + real width end module jt65_mod diff --git a/lib/jt65_test.f90 b/lib/jt65_test.f90 index 00381d3cd..40fbc43ba 100644 --- a/lib/jt65_test.f90 +++ b/lib/jt65_test.f90 @@ -29,14 +29,14 @@ contains call timer('jt65a ',0) call my_decoder%decode(my_callback,dd,npts=52*12000,newdat=.true., & nutc=nutc,nf1=nflow,nf2=nfhigh,nfqso=nfqso,ntol=ntol, & - nsubmode=nsubmode, minsync=0,nagain=.false.,n2pass=n2pass, & + nsubmode=nsubmode, minsync=-1,nagain=.false.,n2pass=n2pass, & nrobust=nrobust,ntrials=ntrials,naggressive=naggressive, & ndepth=ndepth,nclearave=nclearave,mycall=mycall,hiscall=hiscall, & hisgrid=hisgrid,nexp_decode=nexp_decode) call timer('jt65a ',1) end subroutine test - subroutine my_callback (this, utc, sync, snr, dt, freq, drift, decoded & + subroutine my_callback (this,utc,sync,snr,dt,freq,drift,width,decoded & , ft, qual, smo, sum, minsync, submode, aggression) use jt65_decode implicit none @@ -48,6 +48,7 @@ contains real, intent(in) :: dt integer, intent(in) :: freq integer, intent(in) :: drift + real, intent(in) :: width character(len=22), intent(in) :: decoded integer, intent(in) :: ft integer, intent(in) :: qual @@ -56,11 +57,14 @@ contains integer, intent(in) :: minsync integer, intent(in) :: submode integer, intent(in) :: aggression + integer nwidth + nwidth=max(nint(width),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,float(freq),drift,decoded,ft -1012 format(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22,' JT65',i4) + write(13,1012) utc,nint(sync),snr,dt,freq,drift,nwidth, & + decoded,ft,sum,smo +1012 format(i4.4,i4,i5,f6.2,i5,i4,i3,1x,a22,' JT65',3i3) call flush(6) ! write(79,3001) utc,sync,snr,dt,freq,candidates, & ! hard_min,total_min,rtt,tries,ft,qual,decoded diff --git a/lib/jt65sim.f90 b/lib/jt65sim.f90 index bca8793fd..17ce50b61 100644 --- a/lib/jt65sim.f90 +++ b/lib/jt65sim.f90 @@ -213,7 +213,7 @@ program jt65sim a=0. if(x.lt.3.0) then !Cutoff beyond x=3 ! a=sqrt(exp(-x*x)) !Gaussian - a=sqrt(1.0/(1.0+x*x)) !Lorentzian + a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian call random_number(r1) phi1=twopi*r1 z=a*cmplx(cos(phi1),sin(phi1))