Geology Reference
In-Depth Information
34 CONTINUE
C Divide by segment length to get right-hand side of Parseval's relation,
C then divide by left side of Parseval's relation to convert SVL to
C Parseval's ratio.
SVL=SVL/(T*SUMS)
C Construct approximation, GP, to original data set G, from DFT SV.
DO 36 J=1,N
GP(J)=(0.D0,0.D0)
DO 37 I=1,M
AAI=DFLOAT(I-MP-1)
ARG=OMEGA*AAI*DAY(J)/T
GP(J)=GP(J)+SV(I)*DCMPLX(DCOS(ARG),DSIN(ARG))
37 CONTINUE
C Normalize reconstructed data set.
GP(J)=GP(J)/T
36 CONTINUE
C Write out original and reconstructed data sets and calculate mean square
C error in reconstruction ERROR.
ERROR=0.D0
DO 38 J=1,N
ERRORC=G(J)-GP(J)
ERROR=ERROR+ERRORC*DCONJG(ERRORC)
WRITE(7,39)DAY(J),G(J),GP(J)
39 FORMAT(1X,F12.4,4D15.7)
38 CONTINUE
C Find mean square error.
ERROR=ERROR/AN
C Find relative reconstruction error.
ERROR=DSQRT(ERROR)/RMS
C Write out number of singular values eliminated (NSVE),
C relative reconstruction error (RRE), and Parseval's ratio (PR).
WRITE(2,40)NSVE,ERROR,SVL
WRITE(6,40)NSVE,ERROR,SVL
40 FORMAT(1X,I10,2D15.8)
C Decide whether to try new value of NSVE or to end programme.
WRITE(6,41)
41 FORMAT(1X,'Enter 1 for new NSVE or 0 to end.')
READ(5,*)ISWITCH
IF(ISWITCH.NE.0)WRITE(6,42)
IF(ISWITCH.EQ.0)GO TO 43
42 FORMAT(1X,'Enter new value of NSVE.')
READ(5,*)NSVE
GO TO 30
43 CONTINUE
C Write out recovered DFT.
DO 44 J=1,M
AJ=DFLOAT(J-MP-1)
FREQ=AJ/T
WRITE(3,45)FREQ,SV(J)
45 FORMAT(1X,D15.5,2D26.15)
44 CONTINUE
END
The series of 3560 pole positions and their standard errors, running from
3 August 1979 to 22 October 2009, in the file ppy2009a.dat, spans a record length
of 30.218 years. The initial part of the record is very sparsely sampled and we take
the record to be analysed to begin at 4.0 years into the record at point number 48.
The end of the record is set at 30.25 years, the last sample being at point 3560. The
record then has 3513 time points and we use the programme SVDDFT to find the
Search WWH ::




Custom Search