Experimental changes to JT65 decoder:

1. Change to rectangular FFT window for 2D sync spectrum (ss).
2. Move 2D sync spectrum array to common block.
3. Change to quarter-symbol steps for the ss array.
4. Allow up to 4 decoding passes.
5. Wire up Fast/Normal/Deep for 2, 3, or 4 decoding passes.
6. Make minsmo=0 (instead of 1) for minimally spread JT65B/C signals.

 


git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8178 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2017-10-22 00:09:01 +00:00
parent 7149d138f5
commit e25dd201a3
8 changed files with 84 additions and 49 deletions

View File

@ -97,8 +97,8 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials, &
minsmo=0 minsmo=0
maxsmo=0 maxsmo=0
if(mode65.ge.2 .and. mode65.ne.101) then if(mode65.ge.2 .and. mode65.ne.101) then
minsmo=nint(width/df) minsmo=nint(width/df)-1
maxsmo=2*minsmo maxsmo=2*nint(width/df)
endif endif
nn=0 nn=0
do ismo=minsmo,maxsmo do ismo=minsmo,maxsmo

View File

@ -11,7 +11,6 @@ subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
real*4 s3(64,63),afac1 real*4 s3(64,63),afac1
real*8 fs(64)
integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63) integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63)
if(nadd.eq.-999) return if(nadd.eq.-999) return
@ -29,7 +28,6 @@ subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
psum=0. psum=0.
do i=1,64 do i=1,64
x=min(afac*s3(i,j)/ave,50.d0) x=min(afac*s3(i,j)/ave,50.d0)
fs(i)=exp(x)
psum=psum+s3(i,j) psum=psum+s3(i,j)
if(s3(i,j).gt.s1) then if(s3(i,j).gt.s1) then
s1=s3(i,j) s1=s3(i,j)

View File

@ -57,7 +57,7 @@ contains
character(len=6), intent(in) :: hisgrid character(len=6), intent(in) :: hisgrid
real dd(NZMAX) real dd(NZMAX)
real ss(322,NSZ) real ss(552,NSZ)
real savg(NSZ) real savg(NSZ)
real a(5) real a(5)
character*22 decoded,decoded0,avemsg,deepave character*22 decoded,decoded0,avemsg,deepave
@ -81,6 +81,7 @@ contains
real r0(0:11) real r0(0:11)
common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano
common/steve/thresh0 common/steve/thresh0
common/sync/ss
! 0 1 2 3 4 5 6 7 8 9 10 11 ! 0 1 2 3 4 5 6 7 8 9 10 11
data h0/41,42,43,43,44,45,46,47,48,48,49,49/ data h0/41,42,43,43,44,45,46,47,48,48,49,49/
@ -119,23 +120,38 @@ contains
go to 900 go to 900
endif endif
do ipass=1,n2pass !Two-pass decoding loop ! do ipass=1,n2pass !Two-pass decoding loop
npass=1
if(n2pass .gt. 1) npass=ndepth+1 !**** TEMPORARY ****
do ipass=1,npass
first_time=.true. first_time=.true.
if(ipass.eq.1) then !First-pass parameters if(ipass.eq.1) then !First-pass parameters
thresh0=2.5 thresh0=2.5
nsubtract=1 nsubtract=1
nrob=0
elseif( ipass.eq.2 ) then !Second-pass parameters elseif( ipass.eq.2 ) then !Second-pass parameters
thresh0=2.5 thresh0=2.0
nsubtract=1
nrob=0
elseif( ipass.eq.3 ) then
thresh0=2.0
nsubtract=1
nrob=0
elseif( ipass.eq.4 ) then
thresh0=2.0
nsubtract=0 nsubtract=0
nrob=1
endif
if(npass.eq.1) then
nsubtract=0
thresh0=2.0
endif endif
if(n2pass.lt.2) nsubtract=0
! if(newdat) then
call timer('symsp65 ',0) call timer('symsp65 ',0)
ss=0. ss=0.
call symspec65(dd,npts,ss,nhsym,savg) !Get normalized symbol spectra ! call symspec65(dd,npts,ss,nqsym,savg) !Get normalized symbol spectra
call symspec65(dd,npts,nqsym,savg) !Get normalized symbol spectra
call timer('symsp65 ',1) call timer('symsp65 ',1)
! endif
nfa=nf1 nfa=nf1
nfb=nf2 nfb=nf2
single_decode=iand(nexp_decode,32).ne.0 .or. nagain single_decode=iand(nexp_decode,32).ne.0 .or. nagain
@ -161,7 +177,8 @@ contains
ncand=0 ncand=0
call timer('sync65 ',0) call timer('sync65 ',0)
call sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,0,bVHF) ! call sync65(ss,nfa,nfb,naggressive,ntol,nqsym,ca,ncand,0,bVHF)
call sync65(nfa,nfb,naggressive,ntol,nqsym,ca,ncand,nrob,bVHF)
call timer('sync65 ',1) call timer('sync65 ',1)
! If a candidate was found within +/- ntol of nfqso, move it into ca(1). ! If a candidate was found within +/- ntol of nfqso, move it into ca(1).
@ -171,9 +188,6 @@ contains
if(abs(ca(1)%freq - f0).gt.width) width=2*df !### ??? ### if(abs(ca(1)%freq - f0).gt.width) width=2*df !### ??? ###
endif endif
nvec=ntrials nvec=ntrials
if(ncand.gt.75) then
nvec=100
endif
mode65=2**nsubmode mode65=2**nsubmode
nflip=1 nflip=1
@ -195,11 +209,11 @@ contains
ca(ncand)%dt=2.5 ca(ncand)%dt=2.5
ca(ncand)%freq=nfqso ca(ncand)%freq=nfqso
endif endif
do icand=1,ncand do icand=1,ncand
sync1=ca(icand)%sync sync1=ca(icand)%sync
dtx=ca(icand)%dt dtx=ca(icand)%dt
freq=ca(icand)%freq freq=ca(icand)%freq
!write(*,*) icand,sync1,dtx,freq,ndepth,bVHF,mode65
if(bVHF) then if(bVHF) then
flip=ca(icand)%flip flip=ca(icand)%flip
nflip=flip nflip=flip
@ -314,9 +328,8 @@ contains
if(decoded0.eq.' ') decoded0='*' if(decoded0.eq.' ') decoded0='*'
endif endif
enddo !Candidate loop enddo !Candidate loop
if(ndecoded.lt.1) exit if(ipass.eq.2 .and. ndecoded.lt.1) exit
enddo !Two-pass loop enddo !Two-pass loop
900 return 900 return
end subroutine decode end subroutine decode

View File

@ -12,7 +12,7 @@ subroutine slope(y,npts,xpk)
sumxy=0. sumxy=0.
sumy2=0. sumy2=0.
do i=1,npts do i=1,npts
if(abs(i-xpk).gt.2.0) then if(abs(i-xpk).gt.4.0) then
sumw=sumw + 1.0 sumw=sumw + 1.0
x=i x=i
sumx=sumx + x sumx=sumx + x

View File

@ -1,13 +1,16 @@
subroutine symspec65(dd,npts,ss,nhsym,savg) !subroutine symspec65(dd,npts,ss,nqsym,savg)
subroutine symspec65(dd,npts,nqsym,savg)
! Compute JT65 symbol spectra at half-symbol steps ! Compute JT65 symbol spectra at quarter-symbol steps
parameter (NFFT=8192) parameter (NFFT=8192)
parameter (NSZ=3413) !NFFT*5000/12000 parameter (NSZ=3413) !NFFT*5000/12000
parameter (MAXHSYM=322) parameter (MAXHSYM=322)
parameter (MAXQSYM=552)
real*8 hstep real*8 hstep
real*4 dd(npts) real*4 dd(npts)
real*4 ss(MAXHSYM,NSZ) ! real*4 ss(MAXHSYM,NSZ)
real*4 ss(MAXQSYM,NSZ)
real*4 savg(NSZ) real*4 savg(NSZ)
real*4 x(NFFT) real*4 x(NFFT)
real*4 w(NFFT) real*4 w(NFFT)
@ -17,27 +20,33 @@ subroutine symspec65(dd,npts,ss,nhsym,savg)
equivalence (x,c) equivalence (x,c)
data first/.true./ data first/.true./
save /refspec/,first,w save /refspec/,first,w
common/sync/ss
hstep=2048.d0*12000.d0/11025.d0 !half-symbol = 2229.116 samples hstep=2048.d0*12000.d0/11025.d0 !half-symbol = 2229.116 samples
qstep=hstep/2.0 !quarter-symbol = 1114.558 samples
nsps=nint(2*hstep) nsps=nint(2*hstep)
df=12000.0/NFFT df=12000.0/NFFT
nhsym=(npts-NFFT)/hstep nhsym=(npts-NFFT)/hstep
nqsym=(npts-NFFT)/qstep
savg=0. savg=0.
fac1=1.e-3 fac1=1.e-3
if(first) then if(first) then
! Compute the FFT window ! Compute the FFT window
pi=4.0*atan(1.0) ! width=0.25*nsps
width=0.25*nsps
do i=1,NFFT do i=1,NFFT
z=(i-NFFT/2)/width ! z=(i-NFFT/2)/width
w(i)=exp(-z*z) w(i)=1
if(i.gt.4458) w(i)=0
! w(i)=exp(-z*z)
enddo enddo
first=.false. first=.false.
endif endif
do j=1,nhsym ! do j=1,nhsym
i0=(j-1)*hstep do j=1,nqsym
! i0=(j-1)*hstep
i0=(j-1)*qstep
x=fac1*w*dd(i0+1:i0+NFFT) x=fac1*w*dd(i0+1:i0+NFFT)
call four2a(c,NFFT,1,-1,0) !r2c forward FFT call four2a(c,NFFT,1,-1,0) !r2c forward FFT
do i=1,NSZ do i=1,NSZ
@ -48,11 +57,12 @@ subroutine symspec65(dd,npts,ss,nhsym,savg)
enddo enddo
savg=savg/nhsym savg=savg/nhsym
call flat65(ss,nhsym,MAXHSYM,NSZ,ref) !Flatten the 2d spectrum, saving ! call flat65(ss,nhsym,MAXQSYM,NSZ,ref) !Flatten the 2d spectrum, saving
call flat65(ss,nqsym,MAXQSYM,NSZ,ref) !Flatten the 2d spectrum, saving
dfref=df ! the reference spectrum ref() dfref=df ! the reference spectrum ref()
savg=savg/ref savg=savg/ref
do j=1,nhsym ! do j=1,nhsym
do j=1,nqsym
ss(j,1:NSZ)=ss(j,1:NSZ)/ref ss(j,1:NSZ)=ss(j,1:NSZ)/ref
enddo enddo

View File

@ -1,9 +1,12 @@
subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, & !subroutine sync65(ss,nfa,nfb,naggressive,ntol,nqsym,ca,ncand,nrobust, &
! bVHF)
subroutine sync65(nfa,nfb,naggressive,ntol,nqsym,ca,ncand,nrobust, &
bVHF) bVHF)
parameter (NSZ=3413,NFFT=8192,MAXCAND=300) parameter (NSZ=3413,NFFT=8192,MAXCAND=300)
real ss(322,NSZ) ! real ss(322,NSZ)
real ccfblue(-11:540) !CCF with pseudorandom sequence real ss(552,NSZ)
real ccfblue(-32:82) !CCF with pseudorandom sequence
real ccfred(NSZ) !Peak of ccfblue, as function of freq real ccfred(NSZ) !Peak of ccfblue, as function of freq
logical bVHF logical bVHF
@ -16,14 +19,20 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, &
type(candidate) ca(MAXCAND) type(candidate) ca(MAXCAND)
common/steve/thresh0 common/steve/thresh0
common/sync/ss
if(ntol.eq.-99) stop !Silence compiler warning if(ntol.eq.-99) stop !Silence compiler warning
call setup65 call setup65
df=12000.0/NFFT !df = 12000.0/8192 = 1.465 Hz df=12000.0/NFFT !df = 12000.0/8192 = 1.465 Hz
ia=max(2,nint(nfa/df)) ia=max(2,nint(nfa/df))
ib=min(NSZ-1,nint(nfb/df)) ib=min(NSZ-1,nint(nfb/df))
lag1=-11 ! lag1=-11
lag2=59 ! lag2=59
! lag1=-22
! lag2=118
lag1=-32
lag2=82 !may need to be extended for EME
nsym=126 nsym=126
ncand=0 ncand=0
fdot=0. fdot=0.
@ -31,9 +40,9 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, &
ccfblue=0. ccfblue=0.
ccfmax=0. ccfmax=0.
ipk=0 ipk=0
do i=ia,ib do i=ia,ib
call xcor(ss,i,nhsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk0,flip,fdot,nrobust) ! call xcor(ss,i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk0,flip,fdot,nrobust)
call xcor(i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk0,flip,fdot,nrobust)
! Remove best-fit slope from ccfblue and normalize so baseline rms=1.0 ! Remove best-fit slope from ccfblue and normalize so baseline rms=1.0
if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1, & if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1, &
lagpk0-lag1+1.0) lagpk0-lag1+1.0)
@ -64,8 +73,9 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, &
endif endif
endif endif
if(itry.ne.0) then if(itry.ne.0) then
call xcor(ss,i,nhsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot, & ! call xcor(ss,i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot, &
nrobust) ! nrobust)
call xcor(i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot,nrobust)
if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1, & if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1, &
lagpk-lag1+1.0) lagpk-lag1+1.0)
xlag=lagpk xlag=lagpk
@ -73,7 +83,8 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, &
call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2) call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2)
xlag=lagpk+dx2 xlag=lagpk+dx2
endif endif
dtx=xlag*2048.0/11025.0 ! dtx=xlag*2048.0/11025.0
dtx=xlag*1024.0/11025.0
ccfblue(lag1)=0. ccfblue(lag1)=0.
ccfblue(lag2)=0. ccfblue(lag2)=0.
ca(ncand)%freq=freq ca(ncand)%freq=freq

View File

@ -1,4 +1,5 @@
subroutine xcor(ss,ipk,nsteps,nsym,lag1,lag2,ccf,ccf0,lagpk,flip,fdot,nrobust) !subroutine xcor(ss,ipk,nsteps,nsym,lag1,lag2,ccf,ccf0,lagpk,flip,fdot,nrobust)
subroutine xcor(ipk,nsteps,nsym,lag1,lag2,ccf,ccf0,lagpk,flip,fdot,nrobust)
! Computes ccf of a row of ss and the pseudo-random array pr. Returns ! Computes ccf of a row of ss and the pseudo-random array pr. Returns
! peak of the CCF and the lag at which peak occurs. For JT65, the ! peak of the CCF and the lag at which peak occurs. For JT65, the
@ -7,17 +8,19 @@ subroutine xcor(ss,ipk,nsteps,nsym,lag1,lag2,ccf,ccf0,lagpk,flip,fdot,nrobust)
use jt65_mod use jt65_mod
parameter (NHMAX=3413) !Max length of power spectra parameter (NHMAX=3413) !Max length of power spectra
parameter (NSMAX=322) !Max number of half-symbol steps parameter (NSMAX=552) !Max number of quarter-symbol steps
real ss(NSMAX,NHMAX) !2d spectrum, stepped by half-symbols real ss(NSMAX,NHMAX) !2d spectrum, stepped by half-symbols
real a(NSMAX) real a(NSMAX)
real ccf(-11:540) ! real ccf(-44:118)
real ccf(lag1:lag2)
data lagmin/0/ !Silence g77 warning data lagmin/0/ !Silence g77 warning
save ! save
common/sync/ss
df=12000.0/8192. df=12000.0/8192.
dtstep=0.5/df ! dtstep=0.5/df
dtstep=0.25/df
fac=dtstep/(60.0*df) fac=dtstep/(60.0*df)
do j=1,nsteps do j=1,nsteps
ii=nint((j-nsteps/2)*fdot*fac)+ipk ii=nint((j-nsteps/2)*fdot*fac)+ipk
if( (ii.ge.1) .and. (ii.le.NHMAX) ) then if( (ii.ge.1) .and. (ii.le.NHMAX) ) then
@ -43,7 +46,7 @@ subroutine xcor(ss,ipk,nsteps,nsym,lag1,lag2,ccf,ccf0,lagpk,flip,fdot,nrobust)
do lag=lag1,lag2 do lag=lag1,lag2
x=0. x=0.
do i=1,nsym do i=1,nsym
j=2*i-1+lag j=4*i-3+lag
if(j.ge.1 .and. j.le.nsteps) x=x+a(j)*pr(i) if(j.ge.1 .and. j.le.nsteps) x=x+a(j)*pr(i)
enddo enddo
ccf(lag)=2*x !The 2 is for plotting scale ccf(lag)=2*x !The 2 is for plotting scale