Geology Reference
In-Depth Information
The programme SVDDFT finds the DFT of an unequally spaced time sequence.
The name of the sequence file is entered, as is the name of the output file for the
DFT. The first and last point numbers of the data segment are entered, as are the
length and midpoint of the segment. The data are reduced to zero mean and Parzen
windowed to suppress finite record e
ff
ects by a call to the subroutine WINDOW.
The singular value decomposition is accomplished by the subroutine SVD with
calls to the subroutines BIDIAG and HHOLDER. The number of singular values
eliminated (NSVE) to satisfy Parseval's theorem as closely as possible is found
by iteration. The programme is quite general and could be used to find the DFT
of any unequally spaced sequence with suitable adjustment of the input FORMAT
statement.
PROGRAMME SVDDFT.FOR
C
C SVDDFT.FOR finds the DFT of an unequally spaced time sequence of VLBI
C polar motion or nutation observations using singular value decomposition
C (SVD) combined with Parseval's theorem to determine the number of
C singular values to eliminate (NSVE). The sequence is reduced to zero mean,
C and it is then Parzen windowed by a call to the subroutine WINDOW
C to suppress finite record effects. The SVD is found by a call to
C the subroutine SVD which in turn calls the subroutines HHOLDER and BIDIAG.
C The programme allows iteration to find the optimum value of NSVE
C to satisfy Parseval's theorem.
C
IMPLICIT DOUBLE COMPLEX(A-H,O-Z)
DOUBLE PRECISION SIGS(4000),DAY(4000),S(4000),XWOB,YWOB,EXWOB,
1 EYWOB,SVL,T,T0,DT,AM,AN,AJ1,AJ2,ARG,ARG1,ARG2,RMS,SUMX,SUMY,
2 AAI,ERROR,AJ,FREQ,W,PI,SUMS,OMEGA,TD
DIMENSION G(4000),GP(4000),A(4000,4000),SV(4000),U(4000,4000),
1 US(4000,4000),VH(4000,4000),B(4000)
CHARACTER*20 INSEQ,DFTOUT
C Type in name of input sequence file.
WRITE(6,10)
10 FORMAT(1X,'Type in name of input sequence file.')
READ(5,11)INSEQ
11 FORMAT(A20)
C Type in name of output DFT file.
WRITE(6,12)
12 FORMAT(1X,'Type in name of output DFT file.')
READ(5,11)DFTOUT
OPEN(UNIT=1,FILE=INSEQ,STATUS='OLD')
OPEN(UNIT=2,FILE='err.dat',STATUS='UNKNOWN')
OPEN(UNIT=3,FILE=DFTOUT,STATUS='UNKNOWN')
OPEN(UNIT=4,FILE='sv.dat',STATUS='UNKNOWN')
OPEN(UNIT=7,FILE='rec.dat',STATUS='UNKNOWN')
C Set values of pi and 2pi.
PI=3.141592653589793D0
OMEGA=2.D0*PI
C Enter segment length and time at segment midpoint, from beginning of record.
WRITE(6,13)
13 FORMAT(1X,'Enter segment length and time at segment midpoint.')
READ(5,*)T,T0
C Enter first and last point numbers of segment data set.
WRITE(6,14)
14
FORMAT(1X,'Enter first and last point numbers of segment.')
READ(5,*)IB,IE
Search WWH ::
Custom Search