SVN r8568

This commit is contained in:
Jordan Sherer
2018-03-17 12:56:24 -04:00
parent 587950f372
commit 45cc6416c1
633 changed files with 186 additions and 366401 deletions
@@ -1,73 +0,0 @@
subroutine analytic(d,npts,nfft,c,pc,beq)
! Convert real data to analytic signal
parameter (NFFTMAX=1024*1024)
real d(npts) ! passband signal
real h(NFFTMAX/2) ! real BPF magnitude
real*8 pc(5),pclast(5) ! static phase coeffs
real*8 ac(5),aclast(5) ! amp coeffs
real*8 fp
complex corr(NFFTMAX/2) ! complex frequency-dependent correction
complex c(NFFTMAX) ! analytic signal
logical*1 beq ! boolean static equalizer flag
data nfft0/0/
data aclast/1.0,0.0,0.0,0.0,0.0/
data ac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000
save corr,nfft0,h,ac,aclast,pclast,pi,t,beta
df=12000.0/nfft
nh=nfft/2
if( nfft.ne.nfft0 ) then
pi=4.0*atan(1.0)
t=1.0/2000.0
beta=0.1
do i=1,nh+1
ff=(i-1)*df
f=ff-1500.0
h(i)=1.0
if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then
h(i)=h(i)*0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t))))
elseif( abs(f) .gt. (1+beta)/(2*t) ) then
h(i)=0.0
endif
enddo
nfft0=nfft
endif
if( any(aclast .ne. ac) .or. any(pclast .ne. pc) ) then
aclast=ac
pclast=pc
write(*,3001) pc
3001 format('Phase coeffs:',5f12.6)
do i=1,nh+1
ff=(i-1)*df
f=ff-1500.0
fp=f/1000.0
corr(i)=ac(1)+fp*(ac(2)+fp*(ac(3)+fp*(ac(4)+fp*ac(5))))
pd=fp*fp*(pc(3)+fp*(pc(4)+fp*pc(5))) ! ignore 1st two terms
corr(i)=corr(i)*cmplx(cos(pd),sin(pd))
enddo
endif
fac=2.0/nfft
c(1:npts)=fac*d(1:npts)
c(npts+1:nfft)=0.
call four2a(c,nfft,1,-1,1) !Forward c2c FFT
if( beq ) then
c(1:nh+1)=h(1:nh+1)*corr(1:nh+1)*c(1:nh+1)
else
c(1:nh+1)=h(1:nh+1)*c(1:nh+1)
endif
c(1)=0.5*c(1) !Half of DC term
c(nh+2:nfft)=0. !Zero the negative frequencies
call four2a(c,nfft,1,1,1) !Inverse c2c FFT
return
end subroutine analytic
@@ -1,92 +0,0 @@
subroutine sync65(ss,nfa,nfb,naggressive,ntol,nhsym,ca,ncand,nrobust, &
bVHF)
parameter (NSZ=3413,NFFT=8192,MAXCAND=300)
real ss(322,NSZ)
real ccfblue(-11:540) !CCF with pseudorandom sequence
real ccfred(NSZ) !Peak of ccfblue, as function of freq
logical bVHF
type candidate
real freq
real dt
real sync
real flip
end type candidate
type(candidate) ca(MAXCAND)
common/steve/thresh0
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
nsym=126
ncand=0
fdot=0.
ccfred=0.
ccfblue=0.
ccfmax=0.
ipk=0
do i=ia,ib
call xcor(ss,i,nhsym,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)
ccfred(i)=ccfblue(lagpk0)
if(ccfred(i).gt.ccfmax) then
ccfmax=ccfred(i)
ipk=i
endif
enddo
call pctile(ccfred(ia:ib),ib-ia+1,35,xmed)
ccfred(ia:ib)=ccfred(ia:ib)-xmed
ccfred(ia-1)=ccfred(ia)
ccfred(ib+1)=ccfred(ib)
do i=ia,ib
freq=i*df
itry=0
! if(naggressive.gt.0 .and. ntol.lt.1000 .and. ccfmax.ge.thresh0) then
if(naggressive.gt.0 .and. ccfmax.ge.thresh0) then
if(i.ne.ipk) cycle
itry=1
ncand=ncand+1
else
if(ccfred(i).ge.thresh0 .and. ccfred(i).gt.ccfred(i-1) .and. &
ccfred(i).gt.ccfred(i+1)) then
itry=1
ncand=ncand+1
endif
endif
if(itry.ne.0) then
call xcor(ss,i,nhsym,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
if(lagpk.gt.lag1 .and. lagpk.lt.lag2) then
call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2)
xlag=lagpk+dx2
endif
dtx=xlag*2048.0/11025.0
ccfblue(lag1)=0.
ccfblue(lag2)=0.
ca(ncand)%freq=freq
ca(ncand)%dt=dtx
ca(ncand)%flip=flip
if(bVHF) then
ca(ncand)%sync=db(ccfred(i)) - 16.0
else
ca(ncand)%sync=ccfred(i)
endif
endif
if(ncand.eq.MAXCAND) exit
enddo
return
end subroutine sync65