git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@337 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

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

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