diff --git a/lib/decode65a.f90 b/lib/decode65a.f90 index 5ad316771..c53d49dc6 100644 --- a/lib/decode65a.f90 +++ b/lib/decode65a.f90 @@ -97,8 +97,8 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials, & minsmo=0 maxsmo=0 if(mode65.ge.2 .and. mode65.ne.101) then - minsmo=nint(width/df) - maxsmo=2*minsmo + minsmo=nint(width/df)-1 + maxsmo=2*nint(width/df) endif nn=0 do ismo=minsmo,maxsmo diff --git a/lib/demod64a.f90 b/lib/demod64a.f90 index 4eaf0c1d2..b1176efb7 100644 --- a/lib/demod64a.f90 +++ b/lib/demod64a.f90 @@ -11,7 +11,6 @@ subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow) implicit real*8 (a-h,o-z) real*4 s3(64,63),afac1 - real*8 fs(64) integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63) if(nadd.eq.-999) return @@ -29,7 +28,6 @@ subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow) psum=0. do i=1,64 x=min(afac*s3(i,j)/ave,50.d0) - fs(i)=exp(x) psum=psum+s3(i,j) if(s3(i,j).gt.s1) then s1=s3(i,j) diff --git a/lib/ftrsd/ftrsd2.c b/lib/ftrsd/ftrsd2.c index 241312483..6da9b417b 100644 --- a/lib/ftrsd/ftrsd2.c +++ b/lib/ftrsd/ftrsd2.c @@ -38,7 +38,7 @@ void ftrsd2_(int mrsym[], int mrprob[], int mr2sym[], int mr2prob[], int nera_best=0; float pp,pp1,pp2; static unsigned int nseed; - + // Power-percentage symbol metrics - composite gnnf/hf int perr[8][8] = { { 4, 9, 11, 13, 14, 14, 15, 15}, diff --git a/lib/jt65_decode.f90 b/lib/jt65_decode.f90 index fe0a1edec..df136fde9 100644 --- a/lib/jt65_decode.f90 +++ b/lib/jt65_decode.f90 @@ -57,7 +57,7 @@ contains character(len=6), intent(in) :: hisgrid real dd(NZMAX) - real ss(322,NSZ) + real ss(552,NSZ) real savg(NSZ) real a(5) character*22 decoded,decoded0,avemsg,deepave @@ -81,6 +81,7 @@ contains real r0(0:11) common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano common/steve/thresh0 + common/sync/ss ! 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/ @@ -119,23 +120,38 @@ contains go to 900 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. if(ipass.eq.1) then !First-pass parameters thresh0=2.5 nsubtract=1 + nrob=0 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 + nrob=1 + endif + if(npass.eq.1) then + nsubtract=0 + thresh0=2.0 endif - if(n2pass.lt.2) nsubtract=0 - ! if(newdat) then call timer('symsp65 ',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) - ! endif nfa=nf1 nfb=nf2 single_decode=iand(nexp_decode,32).ne.0 .or. nagain @@ -161,7 +177,8 @@ contains ncand=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) ! 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 !### ??? ### endif nvec=ntrials - if(ncand.gt.75) then - nvec=100 - endif mode65=2**nsubmode nflip=1 @@ -195,11 +209,11 @@ contains ca(ncand)%dt=2.5 ca(ncand)%freq=nfqso endif - do icand=1,ncand sync1=ca(icand)%sync dtx=ca(icand)%dt freq=ca(icand)%freq +!write(*,*) icand,sync1,dtx,freq,ndepth,bVHF,mode65 if(bVHF) then flip=ca(icand)%flip nflip=flip @@ -314,9 +328,8 @@ contains if(decoded0.eq.' ') decoded0='*' endif enddo !Candidate loop - if(ndecoded.lt.1) exit + if(ipass.eq.2 .and. ndecoded.lt.1) exit enddo !Two-pass loop - 900 return end subroutine decode diff --git a/lib/slope.f90 b/lib/slope.f90 index 0b2b98da0..b374ee801 100644 --- a/lib/slope.f90 +++ b/lib/slope.f90 @@ -12,7 +12,7 @@ subroutine slope(y,npts,xpk) sumxy=0. sumy2=0. do i=1,npts - if(abs(i-xpk).gt.2.0) then + if(abs(i-xpk).gt.4.0) then sumw=sumw + 1.0 x=i sumx=sumx + x diff --git a/lib/symspec65.f90 b/lib/symspec65.f90 index e9879edf8..522fdadae 100644 --- a/lib/symspec65.f90 +++ b/lib/symspec65.f90 @@ -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 (NSZ=3413) !NFFT*5000/12000 parameter (MAXHSYM=322) + parameter (MAXQSYM=552) real*8 hstep 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 x(NFFT) real*4 w(NFFT) @@ -17,27 +20,33 @@ subroutine symspec65(dd,npts,ss,nhsym,savg) equivalence (x,c) data first/.true./ save /refspec/,first,w + common/sync/ss 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) df=12000.0/NFFT nhsym=(npts-NFFT)/hstep + nqsym=(npts-NFFT)/qstep savg=0. fac1=1.e-3 if(first) then ! Compute the FFT window - pi=4.0*atan(1.0) - width=0.25*nsps +! width=0.25*nsps do i=1,NFFT - z=(i-NFFT/2)/width - w(i)=exp(-z*z) +! z=(i-NFFT/2)/width + w(i)=1 + if(i.gt.4458) w(i)=0 +! w(i)=exp(-z*z) enddo first=.false. endif - do j=1,nhsym - i0=(j-1)*hstep +! do j=1,nhsym + do j=1,nqsym +! i0=(j-1)*hstep + i0=(j-1)*qstep x=fac1*w*dd(i0+1:i0+NFFT) call four2a(c,NFFT,1,-1,0) !r2c forward FFT do i=1,NSZ @@ -48,11 +57,12 @@ subroutine symspec65(dd,npts,ss,nhsym,savg) enddo 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() - savg=savg/ref - do j=1,nhsym +! do j=1,nhsym + do j=1,nqsym ss(j,1:NSZ)=ss(j,1:NSZ)/ref enddo diff --git a/lib/sync65.f90 b/lib/sync65.f90 index 794d5593f..c759d5481 100644 --- a/lib/sync65.f90 +++ b/lib/sync65.f90 @@ -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) parameter (NSZ=3413,NFFT=8192,MAXCAND=300) - real ss(322,NSZ) - real ccfblue(-11:540) !CCF with pseudorandom sequence +! real ss(322,NSZ) + real ss(552,NSZ) + real ccfblue(-32:82) !CCF with pseudorandom sequence real ccfred(NSZ) !Peak of ccfblue, as function of freq logical bVHF @@ -16,14 +19,20 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, & type(candidate) ca(MAXCAND) common/steve/thresh0 + common/sync/ss if(ntol.eq.-99) stop !Silence compiler warning call setup65 + df=12000.0/NFFT !df = 12000.0/8192 = 1.465 Hz ia=max(2,nint(nfa/df)) ib=min(NSZ-1,nint(nfb/df)) - lag1=-11 - lag2=59 +! lag1=-11 +! lag2=59 +! lag1=-22 +! lag2=118 + lag1=-32 + lag2=82 !may need to be extended for EME nsym=126 ncand=0 fdot=0. @@ -31,9 +40,9 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, & ccfblue=0. ccfmax=0. ipk=0 - 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 if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1, & lagpk0-lag1+1.0) @@ -64,8 +73,9 @@ subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, & endif endif if(itry.ne.0) then - call xcor(ss,i,nhsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot, & - nrobust) +! call xcor(ss,i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot, & +! nrobust) + call xcor(i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot,nrobust) if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1, & lagpk-lag1+1.0) 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) xlag=lagpk+dx2 endif - dtx=xlag*2048.0/11025.0 +! dtx=xlag*2048.0/11025.0 + dtx=xlag*1024.0/11025.0 ccfblue(lag1)=0. ccfblue(lag2)=0. ca(ncand)%freq=freq diff --git a/lib/xcor.f90 b/lib/xcor.f90 index 3f95f768a..83937736a 100644 --- a/lib/xcor.f90 +++ b/lib/xcor.f90 @@ -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 ! 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 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 a(NSMAX) - real ccf(-11:540) +! real ccf(-44:118) + real ccf(lag1:lag2) data lagmin/0/ !Silence g77 warning - save +! save + common/sync/ss df=12000.0/8192. - dtstep=0.5/df +! dtstep=0.5/df + dtstep=0.25/df fac=dtstep/(60.0*df) - do j=1,nsteps ii=nint((j-nsteps/2)*fdot*fac)+ipk 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 x=0. 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) enddo ccf(lag)=2*x !The 2 is for plotting scale