(b) The DFT of a real-time series has the symmetry F(−k) =
F (N−k) = F*(k).
(c) If the DFT of f(n) is F(k), the DFT of f(n − j) is W
kj F(k),
and the IDFT of F(k − m) is W
−mn f(n), where W = e
−2πi/ N .
11. As derived in Eqn 4.3.10, the depth h of an earthquake can be
estimated from the difference in arrival times δt between the direct
P wave and pP, the P wave reflected from the surface, using
δt = (2h cos i)/v where i and v are the incidence angle and velocity.
(a) Express the depth as a function of the parameters δt, v, i.
(b) Find the depth for a measured time difference of 2.7 s and
assumed velocity of 6.8 km/s and incidence angle of 24°.
(c) Use the propagation of errors relation to show how the
uncertainty in depth depends on the uncertainties of the
three parameters.
(d) Use the results of (c) to find the uncertainty in depth corresponding to uncertainties (one standard deviation) of 0.5 s
in time difference, 0.5 km/s in velocity, and 3° in incidence
angle. (Remember to convert to radians.)
Computer problems
C-1. Using the Fourier series coefficients for the step function, derived
in problem 1a, plot the first ten terms of the series and their sum.
Also plot the sum of the first 20 and 30 terms.
C-2. Write a subroutine to prepare a time series for taking the fast
Fourier transform and take it. The subroutine should call a set
of separate subroutines that extend the time series to a power of 2
as required, allow for a taper of a length which you input, take
the FFT using the subroutine (COOLB) provided (Box 6C-2)
or another, and plot the amplitude spectrum. The subroutine
should have the option to list the real and imaginary parts of
the spectrum, and the amplitude and phase spectra, at each
frequency.
C-3. (a) Write a subroutine to generate values of the function sin
2πt
T
from t = 0 to t = T max , where the time step ∆t, the period T,
and the total data length T max are inputs.
(b) Plot this function for ∆t = 0.25, T = 5, T max = 20.
(c) Use the results of C-2 to find the amplitude spectrum, with no
tapering and with 10% and 20% tapering.
(d) Do parts (b) and (c) for ∆t = 0.25, T = 8, T max = 50.
(e) Do parts (b) and (c) for the function
sin
( . ) sin
,
2
5
0 5
2
8
π
π
t
t
+
with ∆t = 0.25, T max = 256.
Problems 413
4
J=J-M
M=M/2
IF(M-2)5,3,3
5
J=J+M
MMAX=2
6
IF(MMAX-N)7,10,10
7
ISTEP=2*MMAX
THETA=SIGNI*6.2831831/FLOAT(MMAX)
SINTH=SIN(THETA/2.)
WSTPR=-2.0 *SINTH*SINTH
WSTPI=SIN(THETA)
WR=1.
WI=0.
DO 9 M=1,MMAX,2
DO 8 I=M,N,ISTEP
J=I+MMAX
TEMPR=WR*DATAI(J)-WI*DATAI(J+1)
TEMPI=WR*DATAI(J+1)+WI*DATAI(J)
DATAI(J)=DATAI(I)-TEMPR
DATAI(J+1)=DATAI(I+1)-TEMPI
DATAI(I)=DATAI(I)+TEMPR
8
DATAI(I+1)=DATAI(I+1)+TEMPI
TEMPR=WR
WR=WR*WSTPR-WI*WSTPI+WR
9
WI=WI*WSTPR+TEMPR*WSTPI+WI
MMAX=ISTEP
GO TO 6
10
RETURN
END
Box 6C-2 COOLB subroutine.
1
SUBROUTINE COOLB(NN,DATAI, SIGNI)
C CLASSIC - BUT USABLE - FFT PROGRAM
C DATAI IS DATA ARRAY, 2*NP REAL NUMBERS REPRESENTING
C NP COMPLEX POINTS, SO EACH PAIR OF POINTS ARE THE
C (REAL, IMAGINARY) PARTS OF A COMPLEX NUMBER.
C NN IS POWER OF TWO, CAN BE FOUND BY
C NN=(ALOG10(FLOAT(NP))/ALOG10(2.))+.99
C TRANSFORM DIRECTION CONTROLLED BY REAL VARIABLE
C SIGNI (SIGN OF EXPONENTIAL):-1. FORWARD, 1. TO
C INVERT.
C DIMENSIONS: IF TIME SERIES HAS TIME INCREMENT DT,
C TRANSFORM HAS DELTA FREQ=1/ (2**NN*DT)
C NOTE: AFTER TAKING INVERSE FFT DIVIDE OUTPUT BY 2**NN
INTEGER NN
REAL SIGNI
DIMENSION DATAI(1)
N=2**(NN+1)
J=1
DO 5 I=1,N,2
IF(I-J)1,2,2
1
TEMPR=DATAI(J)
TEMPI=DATAI(J+1)
DATAI(J)=DATAI(I)
DATAI(J+1)=DATAI(I+1)
DATAI(I)=TEMPR
DATAI(I+1)=TEMPI
2
M=N/2
3
IF(J-M)5,5,4
1 COOLB, written in 1960s vintage Fortran (note the arithmetic IF statements), has been left in original form to illustrate both the persistence of programs that
work and the advantages of subsequent developments in programming practice and documentation (Section A.8.2).
F (N−k) = F*(k).
(c) If the DFT of f(n) is F(k), the DFT of f(n − j) is W
kj F(k),
and the IDFT of F(k − m) is W
−mn f(n), where W = e
−2πi/ N .
11. As derived in Eqn 4.3.10, the depth h of an earthquake can be
estimated from the difference in arrival times δt between the direct
P wave and pP, the P wave reflected from the surface, using
δt = (2h cos i)/v where i and v are the incidence angle and velocity.
(a) Express the depth as a function of the parameters δt, v, i.
(b) Find the depth for a measured time difference of 2.7 s and
assumed velocity of 6.8 km/s and incidence angle of 24°.
(c) Use the propagation of errors relation to show how the
uncertainty in depth depends on the uncertainties of the
three parameters.
(d) Use the results of (c) to find the uncertainty in depth corresponding to uncertainties (one standard deviation) of 0.5 s
in time difference, 0.5 km/s in velocity, and 3° in incidence
angle. (Remember to convert to radians.)
Computer problems
C-1. Using the Fourier series coefficients for the step function, derived
in problem 1a, plot the first ten terms of the series and their sum.
Also plot the sum of the first 20 and 30 terms.
C-2. Write a subroutine to prepare a time series for taking the fast
Fourier transform and take it. The subroutine should call a set
of separate subroutines that extend the time series to a power of 2
as required, allow for a taper of a length which you input, take
the FFT using the subroutine (COOLB) provided (Box 6C-2)
or another, and plot the amplitude spectrum. The subroutine
should have the option to list the real and imaginary parts of
the spectrum, and the amplitude and phase spectra, at each
frequency.
C-3. (a) Write a subroutine to generate values of the function sin
2πt
T
from t = 0 to t = T max , where the time step ∆t, the period T,
and the total data length T max are inputs.
(b) Plot this function for ∆t = 0.25, T = 5, T max = 20.
(c) Use the results of C-2 to find the amplitude spectrum, with no
tapering and with 10% and 20% tapering.
(d) Do parts (b) and (c) for ∆t = 0.25, T = 8, T max = 50.
(e) Do parts (b) and (c) for the function
sin
( . ) sin
,
2
5
0 5
2
8
π
π
t
t
+
with ∆t = 0.25, T max = 256.
Problems 413
4
J=J-M
M=M/2
IF(M-2)5,3,3
5
J=J+M
MMAX=2
6
IF(MMAX-N)7,10,10
7
ISTEP=2*MMAX
THETA=SIGNI*6.2831831/FLOAT(MMAX)
SINTH=SIN(THETA/2.)
WSTPR=-2.0 *SINTH*SINTH
WSTPI=SIN(THETA)
WR=1.
WI=0.
DO 9 M=1,MMAX,2
DO 8 I=M,N,ISTEP
J=I+MMAX
TEMPR=WR*DATAI(J)-WI*DATAI(J+1)
TEMPI=WR*DATAI(J+1)+WI*DATAI(J)
DATAI(J)=DATAI(I)-TEMPR
DATAI(J+1)=DATAI(I+1)-TEMPI
DATAI(I)=DATAI(I)+TEMPR
8
DATAI(I+1)=DATAI(I+1)+TEMPI
TEMPR=WR
WR=WR*WSTPR-WI*WSTPI+WR
9
WI=WI*WSTPR+TEMPR*WSTPI+WI
MMAX=ISTEP
GO TO 6
10
RETURN
END
Box 6C-2 COOLB subroutine.
1
SUBROUTINE COOLB(NN,DATAI, SIGNI)
C CLASSIC - BUT USABLE - FFT PROGRAM
C DATAI IS DATA ARRAY, 2*NP REAL NUMBERS REPRESENTING
C NP COMPLEX POINTS, SO EACH PAIR OF POINTS ARE THE
C (REAL, IMAGINARY) PARTS OF A COMPLEX NUMBER.
C NN IS POWER OF TWO, CAN BE FOUND BY
C NN=(ALOG10(FLOAT(NP))/ALOG10(2.))+.99
C TRANSFORM DIRECTION CONTROLLED BY REAL VARIABLE
C SIGNI (SIGN OF EXPONENTIAL):-1. FORWARD, 1. TO
C INVERT.
C DIMENSIONS: IF TIME SERIES HAS TIME INCREMENT DT,
C TRANSFORM HAS DELTA FREQ=1/ (2**NN*DT)
C NOTE: AFTER TAKING INVERSE FFT DIVIDE OUTPUT BY 2**NN
INTEGER NN
REAL SIGNI
DIMENSION DATAI(1)
N=2**(NN+1)
J=1
DO 5 I=1,N,2
IF(I-J)1,2,2
1
TEMPR=DATAI(J)
TEMPI=DATAI(J+1)
DATAI(J)=DATAI(I)
DATAI(J+1)=DATAI(I+1)
DATAI(I)=TEMPR
DATAI(I+1)=TEMPI
2
M=N/2
3
IF(J-M)5,5,4
1 COOLB, written in 1960s vintage Fortran (note the arithmetic IF statements), has been left in original form to illustrate both the persistence of programs that
work and the advantages of subsequent developments in programming practice and documentation (Section A.8.2).