subroutine sync4d(cd0,i0,ctwk,itwk,sync) ! Compute sync power for a complex, downsampled FT4 signal. include 'ft4_params.f90' parameter(NP=NMAX/NDOWN,NSS=NSPS/NDOWN) complex cd0(0:NP-1) complex csynca(2*NSS),csyncb(2*NSS),csyncc(2*NSS),csyncd(2*NSS) complex csync2(2*NSS) complex ctwk(2*NSS) complex z1,z2,z3,z4 complex zz1,zz2,zz3,zz4 logical first integer icos4a(0:3),icos4b(0:3),icos4c(0:3),icos4d(0:3) data icos4a/0,1,3,2/ data icos4b/1,0,2,3/ data icos4c/2,3,1,0/ data icos4d/3,2,0,1/ data first/.true./ save first,twopi,csynca,csyncb,csyncc,csyncd,fac p(z1)=real(z1*fac)**2 + aimag(z1*fac)**2 !Statement function for power if( first ) then twopi=8.0*atan(1.0) k=1 phia=0.0 phib=0.0 phic=0.0 phid=0.0 do i=0,3 dphia=2*twopi*icos4a(i)/real(NSS) dphib=2*twopi*icos4b(i)/real(NSS) dphic=2*twopi*icos4c(i)/real(NSS) dphid=2*twopi*icos4d(i)/real(NSS) do j=1,NSS/2 csynca(k)=cmplx(cos(phia),sin(phia)) csyncb(k)=cmplx(cos(phib),sin(phib)) csyncc(k)=cmplx(cos(phic),sin(phic)) csyncd(k)=cmplx(cos(phid),sin(phid)) phia=mod(phia+dphia,twopi) phib=mod(phib+dphib,twopi) phic=mod(phic+dphic,twopi) phid=mod(phid+dphid,twopi) k=k+1 enddo enddo first=.false. fac=1.0/(2.0*NSS) endif i1=i0 !four Costas arrays i2=i0+33*NSS i3=i0+66*NSS i4=i0+99*NSS z1=0. z2=0. z3=0. z4=0. if(itwk.eq.1) csync2=ctwk*csynca !Tweak the frequency if(i1.ge.0 .and. i1+4*NSS-1.le.NP-1) then z1=sum(cd0(i1:i1+4*NSS-1:2)*conjg(csync2)) elseif( i1.lt.0 ) then npts=(i1+4*NSS-1)/2 if(npts.le.8) then z1=0. else z1=sum(cd0(0:i1+4*NSS-1:2)*conjg(csync2(2*NSS-npts:))) endif endif if(itwk.eq.1) csync2=ctwk*csyncb !Tweak the frequency if(i2.ge.0 .and. i2+4*NSS-1.le.NP-1) z2=sum(cd0(i2:i2+4*NSS-1:2)*conjg(csync2)) if(itwk.eq.1) csync2=ctwk*csyncc !Tweak the frequency if(i3.ge.0 .and. i3+4*NSS-1.le.NP-1) z3=sum(cd0(i3:i3+4*NSS-1:2)*conjg(csync2)) if(itwk.eq.1) csync2=ctwk*csyncd !Tweak the frequency if(i4.ge.0 .and. i4+4*NSS-1.le.NP-1) then z4=sum(cd0(i4:i4+4*NSS-1:2)*conjg(csync2)) elseif( i4+4*NSS-1.gt.NP-1 ) then npts=(NP-1-i4+1)/2 if(npts.le.8) then z4=0. else z4=sum(cd0(i4:i4+2*npts-1:2)*conjg(csync2(1:npts))) endif endif sync = p(z1) + p(z2) + p(z3) + p(z4) return end subroutine sync4d