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
This commit is contained in:
Joe Taylor 2016-03-17 13:28:57 +00:00
parent 8d21cdf1c7
commit 9f388b63dc
8 changed files with 99 additions and 42 deletions

View File

@ -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

View File

@ -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)

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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))