390 Seismograms as Signals
Commonly, the same computer program is used for both forward and inverse FFTs, except that the sign of the exponential
must be changed and the 1/N normalization remembered (the
last being a traditional bane of students).
In using the FFT to transform data as part of a filtering
operation, the factor of 1/N may be included at any step in the
process. Often, however, we use the FFT to obtain the Fourier
transform of a time series, and compare this to a result derived
in the frequency domain, such as an analytic expression for a
synthetic seismogram as a function of ω. In this case, we have to
consider the units of both the forward and the inverse DFT.
The forward DFT is an approximate way of evaluating the
Fourier transform integral (Eqn 13), in which the differential
dt is replaced by the difference ∆t. Thus, the FFT results are
multiplied by ∆t. Similarly, the IDFT approximates the inverse
transform integral (Eqn 21), with the differential dω replaced
by the difference ∆ω. Hence the results from inverting the FFT
are multiplied by ∆ω/(2π). The product of these two factors is
∆ω∆t/2π = 1/N, as expected.
This discussion assumes that the series length N is a power of
2. If this is not the case, a number of zeroes necessary to obtain
a power of 2 can be added to the end of the time series. Such
zero padding has the effect of sampling the spectrum more
densely, because the sample interval is unchanged, but the
frequency interval ∆ω = 2π/(N∆t) decreases. Despite the denser
sampling, the real resolution in frequency is not increased
beyond that resulting from the real (nonzero) data length.
Instead, smooth interpolation is done within the range of actual resolution ∆ω real = 2π /T nonzero .
Finally, it is worth distinguishing between the DFT and the
FFT. The DFT is the discrete approximation to the Fourier
transform which has the periodic properties we have discussed.
The FFT is a clever method for computing the DFT with many
fewer operations.
6.4.5 Digital convolution
As discussed in Section 6.3.2, the convolution is used in many
seismological applications. This operation has some special
features when carried out with discretized time series and their
transforms.
Given two discrete time series with unit sample period,
x(m) with M points x(0), x(1), . . . , x(M − 1) and f(n) with
N points f(0), f(1), . . . , f(N − 1), the convolution in the
time domain is written, by analogy to the integral definition, as
y(t) = x(t) * f(t) =
m
M
=
−
∑
0
1
x(m)f(t − m).
(40)
We evaluate the summation for each value of t that yields a
nonzero value. Because f(n) is zero for n outside the range (0,
N − 1) and x(m) is zero for m outside the range (0, M − 1), there
are N + M − 1 terms in the convolution, and y(t) is defined for
t = 0, 1, . . . , N + M − 2. For example, if N = 3 and M = 4, the
3 + 4 − 1 = 6 terms are
B(k) =
n
N
=
−
∑
0
2 1
/
b(n)e −4πikn/N ,
(33)
where k goes from 0 to N/2 − 1, and the factor of 4 comes from
the fact that the subseries lengths are N/2.
The DFT of the original series can be written in terms of the
DFTs of the subseries,
F(k) =
n
N
=
−
∑
0
1
f(n)e −2π ikn / N
=
n
N
=
−
∑
0
2 1
/
[a(n)e −2πik(2n)/N + b(n)e −2πik(2n+1)/ N ]
= A(k) + e −2πik/N B(k) for k = 0, 1, . . . , N/2 − 1,
(34)
giving the first N/2 points of F(k). The second N/2 points come
from replacing k by k + N/2,
F(k + N/2) = A(k + N/2) + e −2πi(k+N/2)/N B(k + N/2),
(35)
and noting that, because the DFTs of the subseries are periodic
with a period equal to their length, N/2,
A(k + N/2) = A(k) and B(k + N/2) = B(k).
(36)
Because the exponential can be written as
e
−2π i(k+N/2)/N = e −πi e −2πik/N = −e −2πik/N ,
(37)
the second half of the transform can be found from the first,
using
F(k + N/2) = A(k) − e −2πik/N B(k).
(38)
In terms of W = e −2πi/N , the expressions for the two parts of the
transform (Eqns 34 and 38) have the simple form of
F(k) = A(k) + W k B(k) and F(k + N/2) = A(k) − W k B(k).
(39)
This method is called doubling a finding the transform of
an N-point series from the transforms of its two N/2-point
subseries. Doubling can be applied recursively, because we
can find the transform of each N/2-point series from that of two
N/4-point series, etc. Ultimately, a series of length N = 2
n can be
evaluated via n = log 2 N such stages. In the final stage, the transform of each 2-point series is found from two 1-point series,
but the transform of a 1-point series is itself. Various methods
can be used to further speed up operations.
Thus, to obtain the FFT of a time series, we treat the data
points as N 1-point series, use doubling to form (N/2) 2-point
series, and so on until the final N-point transform. The same
FFT algorithm can also be used to take the inverse transform.
Commonly, the same computer program is used for both forward and inverse FFTs, except that the sign of the exponential
must be changed and the 1/N normalization remembered (the
last being a traditional bane of students).
In using the FFT to transform data as part of a filtering
operation, the factor of 1/N may be included at any step in the
process. Often, however, we use the FFT to obtain the Fourier
transform of a time series, and compare this to a result derived
in the frequency domain, such as an analytic expression for a
synthetic seismogram as a function of ω. In this case, we have to
consider the units of both the forward and the inverse DFT.
The forward DFT is an approximate way of evaluating the
Fourier transform integral (Eqn 13), in which the differential
dt is replaced by the difference ∆t. Thus, the FFT results are
multiplied by ∆t. Similarly, the IDFT approximates the inverse
transform integral (Eqn 21), with the differential dω replaced
by the difference ∆ω. Hence the results from inverting the FFT
are multiplied by ∆ω/(2π). The product of these two factors is
∆ω∆t/2π = 1/N, as expected.
This discussion assumes that the series length N is a power of
2. If this is not the case, a number of zeroes necessary to obtain
a power of 2 can be added to the end of the time series. Such
zero padding has the effect of sampling the spectrum more
densely, because the sample interval is unchanged, but the
frequency interval ∆ω = 2π/(N∆t) decreases. Despite the denser
sampling, the real resolution in frequency is not increased
beyond that resulting from the real (nonzero) data length.
Instead, smooth interpolation is done within the range of actual resolution ∆ω real = 2π /T nonzero .
Finally, it is worth distinguishing between the DFT and the
FFT. The DFT is the discrete approximation to the Fourier
transform which has the periodic properties we have discussed.
The FFT is a clever method for computing the DFT with many
fewer operations.
6.4.5 Digital convolution
As discussed in Section 6.3.2, the convolution is used in many
seismological applications. This operation has some special
features when carried out with discretized time series and their
transforms.
Given two discrete time series with unit sample period,
x(m) with M points x(0), x(1), . . . , x(M − 1) and f(n) with
N points f(0), f(1), . . . , f(N − 1), the convolution in the
time domain is written, by analogy to the integral definition, as
y(t) = x(t) * f(t) =
m
M
=
−
∑
0
1
x(m)f(t − m).
(40)
We evaluate the summation for each value of t that yields a
nonzero value. Because f(n) is zero for n outside the range (0,
N − 1) and x(m) is zero for m outside the range (0, M − 1), there
are N + M − 1 terms in the convolution, and y(t) is defined for
t = 0, 1, . . . , N + M − 2. For example, if N = 3 and M = 4, the
3 + 4 − 1 = 6 terms are
B(k) =
n
N
=
−
∑
0
2 1
/
b(n)e −4πikn/N ,
(33)
where k goes from 0 to N/2 − 1, and the factor of 4 comes from
the fact that the subseries lengths are N/2.
The DFT of the original series can be written in terms of the
DFTs of the subseries,
F(k) =
n
N
=
−
∑
0
1
f(n)e −2π ikn / N
=
n
N
=
−
∑
0
2 1
/
[a(n)e −2πik(2n)/N + b(n)e −2πik(2n+1)/ N ]
= A(k) + e −2πik/N B(k) for k = 0, 1, . . . , N/2 − 1,
(34)
giving the first N/2 points of F(k). The second N/2 points come
from replacing k by k + N/2,
F(k + N/2) = A(k + N/2) + e −2πi(k+N/2)/N B(k + N/2),
(35)
and noting that, because the DFTs of the subseries are periodic
with a period equal to their length, N/2,
A(k + N/2) = A(k) and B(k + N/2) = B(k).
(36)
Because the exponential can be written as
e
−2π i(k+N/2)/N = e −πi e −2πik/N = −e −2πik/N ,
(37)
the second half of the transform can be found from the first,
using
F(k + N/2) = A(k) − e −2πik/N B(k).
(38)
In terms of W = e −2πi/N , the expressions for the two parts of the
transform (Eqns 34 and 38) have the simple form of
F(k) = A(k) + W k B(k) and F(k + N/2) = A(k) − W k B(k).
(39)
This method is called doubling a finding the transform of
an N-point series from the transforms of its two N/2-point
subseries. Doubling can be applied recursively, because we
can find the transform of each N/2-point series from that of two
N/4-point series, etc. Ultimately, a series of length N = 2
n can be
evaluated via n = log 2 N such stages. In the final stage, the transform of each 2-point series is found from two 1-point series,
but the transform of a 1-point series is itself. Various methods
can be used to further speed up operations.
Thus, to obtain the FFT of a time series, we treat the data
points as N 1-point series, use doubling to form (N/2) 2-point
series, and so on until the final N-point transform. The same
FFT algorithm can also be used to take the inverse transform.