Geology Reference
In-Depth Information
DO 13 I=1,N
C Find index.
IN=I-NP-1
AI=DFLOAT(IN)
C Calculate frequency.
FREQ(I)=AI/T
C Test for negative index, replace with positive index one period ahead.
IF(IN.LT.0)IN=IN+N
J=IN+1
C Apply phase shift to find DFT from FFT.
ARG=OMEGA/AN
ARG=ARG*AI*NP
GF(I)=GT(J)*DCMPLX(DCOS(ARG),DSIN(ARG))
C Normalize DFT.
GF(I)=DT*GF(I)
C Calculate spectral density.
SD(I)=GF(I)*DCONJG(GF(I))
SD(I)=SD(I)*560.D0/(151.D0*T)
C Write out spectral density.
WRITE(2,14)FREQ(I),SD(I)
14
FORMAT(1X,D15.5,D26.15)
13
CONTINUE
END
Both the spectral density estimate of the VLBI pole path found by singular value
decomposition (SVD), shown in Figure 4.4, and that found for the interpolated
pole path, shown in Figure 4.6, have only the Chandler resonance and the annual
line rising above the 95% confidence level. The true spectrum is convolved, or
averaged, with the square of the Parzen frequency window, as given in expression
(2.343). The half-power points of the Parzen window, in the present case, are separ-
ated by a distance of 1.82/26.25 cycles per year, or approximately 0.07 cpy. Thus,
the Chandler and annual features are barely separated, suggesting the application
of a super-resolution technique, such as the maximum entropy spectral analysis
considered in Section 2.6.
4.2.4 Maximum entropy spectral analysis of the VLBI pole path
The maximum entropy spectral density, given by expression (2.470), is propor-
tional to the inverse of the squared magnitude of the discrete Fourier transform
(DFT) of the prediction error filter. In turn, the prediction error filter is found using
the Burg algorithm described in Section 2.6.3, and implemented by the subroutine
BPEC.
The programme MEMSPEC finds the spectral density by the maximum entropy
method (MEM), with calls to the subroutine BPEC, and gives it as a function of
frequency in cycles per year (cpy) in the output file memspec.dat. The input file
is the pole path interpolated by natural splines, ppintp.dat, as generated by the
programme PPINTP. The length of the prediction error filter is entered as input on
request.
Search WWH ::




Custom Search