1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
subroutine sync4(dat,jz,ntol,nfqso,mode,mode4,minwidth,dtx,dfx,snrx, &
snrsync,flip,width)
! Synchronizes JT4 data, finding the best-fit DT and DF.
parameter (NFFTMAX=2520) !Max length of FFTs
parameter (NHMAX=NFFTMAX/2) !Max length of power spectra
parameter (NSMAX=525) !Max number of half-symbol steps
integer ntol !Range of DF search
real dat(jz)
real s2(NHMAX,NSMAX) !2d spectrum, stepped by half-symbols
real ccfblue(-5:540) !CCF with pseudorandom sequence
real ccfred(NHMAX) !Peak of ccfblue, as function of freq
real red(NHMAX) !Peak of ccfblue, as function of freq
integer ipk1(1)
integer nch(7)
logical savered
equivalence (ipk1,ipk1a)
data nch/1,2,4,9,18,36,72/
save
! Do FFTs of twice symbol length, stepped by half symbols. Note that
! we have already downsampled the data by factor of 2.
nsym=207
nfft=2520
nh=nfft/2
nq=nfft/4
nsteps=jz/nq - 1
df=0.5*11025.0/nfft
ftop=nfqso + 7*mode4*df
if(ftop.gt.11025.0/4.0) then
print*,'*** Rx Freq is set too high for this submode ***'
go to 900
endif
if(mode.eq.-999) width=0. !Silence compiler warning
do j=1,nsteps !Compute spectrum for each step, get average
k=(j-1)*nq + 1
call ps4(dat(k),nfft,s2(1,j))
enddo
! Set freq and lag ranges
ia=(nfqso-ntol)/df !Index of lowest tone, bottom of search range
ib=(nfqso+ntol)/df !Index of lowest tone, top of search range
iamin=nint(100.0/df)
if(ia.lt.iamin) ia=iamin
ibmax=nint(2700.0/df) - 6*mode4
if(ib.gt.ibmax) ib=ibmax
lag1=-5
lag2=59
syncbest=-1.e30
snrx=-26.0
ccfred=0.
red=0.
i0=nint(nfqso/df)
do ich=minwidth,7 !Find best width
kz=nch(ich)/2
savered=.false.
iaa=ia+kz
ibb=ib-kz
do i=iaa,ibb !Find best frequency channel for CCF
call xcor4(s2,i,nsteps,nsym,lag1,lag2,ich,mode4,ccfblue,ccf0, &
lagpk0,flip)
ccfred(i)=ccf0
! Find rms of the CCF, without main peak
call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0)
sync=abs(ccfblue(lagpk0))
! write(*,3000) ich,i,i*df,ccf0,sync,syncbest
!3000 format(2i5,4f12.3)
! Find best sync value
if(sync.gt.syncbest*1.03) then
ipk=i
lagpk=lagpk0
ichpk=ich
syncbest=sync
savered=.true.
endif
enddo
if(savered) red=ccfred
enddo
if(syncbest.lt.-1.e29) go to 900
ccfred=red
call pctile(ccfred(ia:ib),ib-ia+1,45,base)
ccfred=ccfred-base
dfx=ipk*df
! Peak up in time, at best whole-channel frequency
call xcor4(s2,ipk,nsteps,nsym,lag1,lag2,ichpk,mode4,ccfblue,ccfmax, &
lagpk,flip)
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
! Find rms of the CCF, without the main peak
call slope(ccfblue(lag1),lag2-lag1+1,xlag-lag1+1.0)
sq=0.
nsq=0
do lag=lag1,lag2
if(abs(lag-xlag).gt.2.0) then
sq=sq+ccfblue(lag)**2
nsq=nsq+1
endif
enddo
rms=sqrt(sq/nsq)
snrsync=max(0.0,db(abs(ccfblue(lagpk)/rms - 1.0)) - 4.5)
dt=2.0/11025.0
istart=xlag*nq
dtx=istart*dt
ipk1=maxloc(ccfred)
ccf10=0.5*maxval(ccfred)
do i=ipk1a,ia,-1
if(ccfred(i).le.ccf10) exit
enddo
i1=i
do i=ipk1a,ib
if(ccfred(i).le.ccf10) exit
enddo
nw=i-i1
width=nw*df
sq=0.
ns=0
iaa=max(ipk1a-10*nw,ia)
ibb=min(ipk1a+10*nw,ib)
jmax=2*mode4/3
do i=iaa,ibb
j=abs(i-ipk1a)
if(j.gt.nw .and. j.lt.jmax) then
sq=sq + ccfred(j)*ccfred(j)
ns=ns+1
endif
enddo
rms=sqrt(sq/ns)
snrx=10.0*log10(ccfred(ipk1a)/rms) - 41.2
900 return
end subroutine sync4
|
