diff --git a/ccf65.f b/ccf65.f new file mode 100644 index 000000000..ec705dce5 --- /dev/null +++ b/ccf65.f @@ -0,0 +1,114 @@ + subroutine ccf65(ss,sync1,ipol1,dt1,flipk,syncshort, + + snr2,ipol2,dt2) + + parameter (NFFT=512,NH=NFFT/2) + real ss(4,322) !Input: half-symbol powers, 4 pol'ns + real s(NFFT) !CCF = ss*pr + complex cs(0:NH) !Complex FT of s + real s2(NFFT) !CCF = ss*pr2 + complex cs2(0:NH) !Complex FT of s2 + real pr(NFFT) !JT65 pseudo-random sync pattern + complex cpr(0:NH) !Complex FT of pr + real pr2(NFFT) !JT65 shorthand pattern + complex cpr2(0:NH) !Complex FT of pr2 + real tmp1(322) + real tmp2(322) + real ccf(-27:27,4) + logical first + integer npr(126) + data first/.true./ + equivalence (s,cs),(pr,cpr),(s2,cs2),(pr2,cpr2) + save + +C The JT65 pseudo-random sync pattern: + data npr/ + + 1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, + + 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, + + 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, + + 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, + + 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, + + 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, + + 1,1,1,1,1,1/ + + if(first) then +C Initialize pr, pr2; compute cpr, cpr2. + fac=1.0/NFFT + do i=1,NFFT + pr(i)=0. + k=2*mod((i-1)/8,2)-1 + pr2(i)=fac*k + enddo + do i=1,126 + j=2*i + pr(j)=fac*(2*npr(i)-1) + enddo + call four2a(pr,NFFT,1,-1,0) + call four2a(pr2,NFFT,1,-1,0) + first=.false. + endif + +C Look for JT65 sync pattern and shorthand square-wave pattern. + ccfbest=0. + ccfbest2=0. + do ip=1,4 !Do all four pol'ns + do i=1,321 + s(i)=min(4.0,ss(ip,i)+ss(ip,i+1)) + enddo + do i=322,NFFT + s(i)=0. + enddo + call four2a(s,NFFT,1,-1,0) !Real-to-complex FFT + do i=0,NH + cs2(i)=cs(i)*conjg(cpr2(i)) !Mult by complex FFT of pr2 + cs(i)=cs(i)*conjg(cpr(i)) !Mult by complex FFT of pr + enddo + call four2a(cs,NFFT,1,1,-1) !Complex-to-real inv-FFT + call four2a(cs2,NFFT,1,1,-1) !Complex-to-real inv-FFT + + do lag=-27,27 !Check for best JT65 sync + ccf(lag,ip)=s(lag+28) + if(abs(ccf(lag,ip)).gt.ccfbest) then + ccfbest=abs(ccf(lag,ip)) + lagpk=lag + ipol1=ip + flipk=1.0 + if(ccf(lag,ip).lt.0.0) flipk=-1.0 + endif + enddo + + do lag=-8,7 !Check for best shorthand + ccf2=s2(lag+28) + if(ccf2.gt.ccfbest2) then + ccfbest2=ccf2 + lagpk2=lag + ipol2=ip + endif + enddo + + enddo + +C Find rms level on baseline of "ccfblue", for normalization. + sum=0. + do lag=-26,26 + if(abs(lag-lagpk).gt.1) sum=sum + ccf(lag,ipol1) + enddo + base=sum/50.0 + sq=0. + do lag=-26,26 + if(abs(lag-lagpk).gt.1) sq=sq + (ccf(lag,ipol1)-base)**2 + enddo + rms=sqrt(sq/49.0) + sync1=ccfbest/rms - 4.0 + dt1=2.5 + lagpk*(2048.0/11025.0) + +C Find base level for normalizing snr2. + do i=1,322 + tmp1(i)=ss(ipol2,i) + enddo + call pctile(tmp1,tmp2,322,40,base) + snr2=0.398107*ccfbest2/base !### empirical + syncshort=0.5*ccfbest2/rms - 4.0 !### better normalizer than rms? + dt2=2.5 + lagpk2*(2048.0/11025.0) + + return + end diff --git a/fchisq.f b/fchisq.f new file mode 100644 index 000000000..f6e051ecc --- /dev/null +++ b/fchisq.f @@ -0,0 +1,72 @@ + real function fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax) + + parameter (NMAX=60*96000) !Samples per 60 s + complex cx(npts),cy(npts) + real a(5) + complex w,wstep,za,zb,z + real ss(2600) + complex csx(0:NMAX/64),csy(0:NMAX/64) + data twopi/6.283185307/a1,a2,a3/99.,99.,99./ + save + + baud=11025.0/4096.0 + if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then + a1=a(1) + a2=a(2) + a3=a(3) + +C Mix and integrate the complex X and Y signals + csx(0)=0. + csy(0)=0. + w=1.0 + x0=0.5*(npts+1) + s=2.0/npts + do i=1,npts + x=s*(i-x0) + if(mod(i,100).eq.1) then + p2=1.5*x*x - 0.5 +! p3=2.5*(x**3) - 1.5*x +! p4=4.375*(x**4) - 3.75*(x**2) + 0.375 + dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample) + wstep=cmplx(cos(dphi),sin(dphi)) + endif + w=w*wstep + csx(i)=csx(i-1) + w*cx(i) + csy(i)=csy(i-1) + w*cy(i) + enddo + endif + +C Compute 1/2-symbol powers at 1/16-symbol steps. + fac=1.e-4 + pol=a(4)/57.2957795 + aa=cos(pol) + bb=sin(pol) + nsps=nint(fsample/baud) !Samples per symbol + nsph=nsps/2 !Samples per half-symbol + + ndiv=16 !Output ss() steps per symbol + nout=ndiv*npts/nsps + dtstep=1.0/(ndiv*baud) !Time per output step + + do i=1,nout + j=i*nsps/ndiv + k=j-nsph + ss(i)=0. + if(k.ge.1) then + za=csx(j)-csx(k) + zb=csy(j)-csy(k) + z=aa*za + bb*zb + ss(i)=fac*(real(z)**2 + aimag(z)**2) + endif + enddo + + ccfmax=0. + call ccf2(ss,nout,nflip,ccf,lagpk) + if(ccf.gt.ccfmax) then + ccfmax=ccf + dtmax=lagpk*dtstep + endif + fchisq=-ccfmax + + return + end