Ch2 Fundamentals 2013
Ch2 Fundamentals 2013
Ch2 Fundamentals 2013
Chapter 2
Fundamentals
The Fourier transform (FT)
- Consider a mono-frequency sinusoidal function, g(t), given by:
g(t) = |A
0
| cos(2t f
0
t - |
0
). (2.1)
- g(t) can be completely described in terms of the following parameters:
The peak amplitude |A
0
|, which is the highest amplitude of g(t).
The frequency f
0
, which is the number of cycles per second.
The phase |
0
, which is a time-shift applied to g(t) with respect to t = 0.
- A time-dependent periodic signal, h(t) with a fundamental (dominant) period T
0
, can
be analyzed to an infinite series of sinusoids, each having its own peak-amplitude,
frequency, and phase; using the Fourier series:
=
=
0
) 2 cos( ) (
n
n n n
t f A t h | t
; (2.2)
where A
n
, f
n
, and |
n
are the peak amplitude, frequency, and phase of the n
th
sinusoidal
component and f
0
= 1/T
0
is the fundamental (dominant) frequency of h(t). (Examples
1 and 2; Yilmaz, 2001)
- We use the forward Fourier transform to compute the peak amplitude A
n
and phase |
n
at every frequency f
n
.
- The forward Fourier transform H(f) of the time-function h(t) is given by:
}
+
= dt e t h f H
ft i t 2
) ( ) (
. (2.3)
2
- H(f) is generally a complex function that can be represented as:
)] ( / ) ( [ tan ) (
) ( ) (
)]} ( sin[ )] ( {cos[ ) (
) ( ) (
1
) (
f H f H f
f iH f H
f i f f H
e f H f H
r i
i r
f i
=
+ =
+ =
=
|
| |
|
, (2.4)
where H
r
(f), H
i
(f), |H(f)|, and |(f) are the real part, imaginary part, amplitude
spectrum, and phase spectrum of H(f), respectively.
- Practically, when we carry out the forward Fourier transform, it returns the amplitude
spectrum |H(f)| and phase spectrum |(f) as functions of frequency f.
- Conversely, given |H(f)| and |(f), we can synthesize h(t) using the inverse Fourier
transform:
}
+
= df e f H t h
ft i t 2
) ( ) (
, (2.5)
where H(f) is calculated from |H(f)| and |(f) using equation (2.4).
- The algorithm used commonly to calculate the Fourier transform numerically is the
fast Fourier transform (FFT).
- The Fourier transform will transform any function p(a) from its domain (a) to a
function P(b) in another domain (b) provided that: b = 1/a (e.g., time to frequency).
- In this course, we will denote a Fourier transform pair as: h(t) |H(f)| and we will
only calculate the amplitude spectrum of the function. The phase spectrum will only
be calculated if needed.
- More Internet resources are on: http://www.answers.com/topic/fourier-transform.
3
The impulse (delta) function
- The impulse function o(t) is defined practically as:
)
`
=
=
=
0 0
0
) (
t
t
t o . (2.6a)
- o(t) can be thought of as a function that has a very large magnitude at t = 0 and
infinitesimally small duration such that the area under the curve is one.
- For sampling purposes, o(t) can be treated as a function that has a magnitude of
one at t = 0 and zero elsewhere (Figure, Wikipedia, 2013).
- Generally, for constant k and t
0
, the function k o(t-t
0
) has a magnitude of k at t = t
0
and zero elsewhere.
- o(t) 1, - < f < . That is, the FT of o(t) is constant for all frequencies
(Figure, www.physforum.com, 2013).
- Another function that is related to o(t) is the sampling (Dirac) function D(t) given
as:
)
`
A A A =
A A A =
= A =
=
t N t t t
t N t t t
t n t t D
N
n
,..., 2 , , 0 0
,..., 2 , , 0 1
) ( ) (
0
o
, (2.6b)
where At is the sampling interval and N is the number of samples (Figure,
Wikipedia, 2013).
- D(t) is used for sampling continuous functions by extracting the function values only
at intervals of At and discarding all other values.
- More Internet resources are on: http://www.answers.com/topic/dirac-delta-function.
4
The sinc function
- The sinc function is defined as:
t f
t f
t
0
0
2
) 2 sin(
) sinc(
t
t
=
, (2.7)
where f
0
is the frequency of the sine function.
- The sinc(t) function is basically a mono-frequency sinusoidal function scaled by its
own time.
- sinc(t) 1, |f| < f
0
/2. That is, the FT of sinc(t) has a magnitude of one for
frequencies in the open interval (-f
0
/2, +f
0
/2) and zero elsewhere. It is a square
(boxcar) function (Figure, www.physforum.com, 2013).
- The sinc(t) function is commonly used in the frequency domain for filtering (i.e.,
keeping frequencies in an interval (f
1
, f
2
) and discarding frequencies outside this
range) as we will discuss in frequency filtering.
- More Internet resources are on: http://www.answers.com/sinc%20function.
Convolution
- The convolution of two time-dependent functions x(t) and h(t) yields the time-
dependent function y(t) given by:
}
= = t t t d t h x t h t x t y ) ( ) ( ) ( * ) ( ) (
, (2.8)
where * denotes the convolution operation. Examples (Ex1, Ex2, Brigham, 1988).
- To convolve two sampled functions x(t) and h(t), we:
1. Fold h(t)
2. Shift it one sample to the right
5
3. Multiply overlapping samples of x(t) and h(t)
4. Add the result
5. Repeat steps 1-4 until there is no more overlapping samples of x(t) and h(t).
- For sampled functions x(t) and h(t) that have number of samples N
x
and N
h
,
respectively, the convolution y(t) has a number of samples N
y
given by:
N
y
= N
x
+ N
h
1. (2.9)
- The convolution of a function x(t) with o(t t
0
) is evaluated by simply shifting x(t) so
that its vertical axis is positioned at t = t
0
. That is: x(t) * o(t t
0
) = x(t t
0
).
Examples (Ex1, Ex2, Brigham, 1988).
- The convolution theorem states that convolution of two functions in the time-domain
is equivalent to multiplying their Fourier transforms in the frequency domain:
y(t) = x(t) * h(t) Y(f) = X(f) H(f), (2.10)
where:
Y(f) = |A
y
(f)| e
i|
y
(f)
, |A
y
(f)| and |
y
(f) are amplitude and phase spectra of y(t);
|A
y
(f)| = |A
x
(f)| |A
h
(f)|, |A
x
(f)| and |A
h
(f)| are amplitude spectra of x(t) and h(t);
|
y
(f) = |
x
(f) + |
h
(f), |
x
(f) and |
h
(f) are phase spectra of x(t) and h(t) (Ex1, Ex2,
Brigham, 1988).
- The inverse of the convolution theorem is also true (Example, Brigham, 1988):
x(t) h(t) X(f) * H(f). (2.11)
- Convolution is commutative. That is: x(t) * h(t) = h(t) * x(t). This means that it
doesnt matter if you fold and shift x(t) or h(t). Examples (Ex1, Ex2, Brigham, 1988).
- It is usually faster to perform convolution in the frequency domain for functions of
long durations (i. e., N
x
& N
h
> 32) if the FFT is used.
6
- Numerical example of convolution.
- More Internet resources are on: http://www.answers.com/topic/convolution.
Correlation
- The correlation of two time-dependent functions x(t) and h(t) yields the time-
dependent function y(t) given by:
}
+ = = t t t d t h x t h t x t y ) ( ) ( ) ( ) ( ) (
, (2.12)
where denotes the correlation operation. Examples (Ex1, Ex2, Brigham, 1988).
- To correlate two sampled functions x(t) and h(t), we:
1. Shift h(t) all the way to the left until its last sample overlaps with x(t) first sample
2. Shift h(t) one sample to the right
3. Multiply overlapping samples of x(t) and h(t)
4. Add the result
5. Repeat steps 2-4 until there is no more overlapping samples of x(t) and h(t).
- For sampled functions x(t) and h(t) that have number of samples N
x
and N
h
,
respectively, the correlation y(t) has a number of samples N
y
given by:
N
y
= N
x
+ N
h
1. (2.13)
- The correlation of a function with another different function is called crosscorrelation
while the correlation of a function with itself is called autocorrelation.
- Applying the convolution theorem on the correlation yields:
y(t) = x(t) h(t) = x(t) * h(-t) Y(f) = X(f) H
*
(f) (2.14)
7
where H
*
(f) = H
r
(f) i H
i
(f) = |A
h
(f)|e
-iuh(f)
is the complex conjugate of H(f) = H
r
(f) +
i H
i
(f) = |A
h
(f)|e
iuh(f)
; Y(f) = |A
y
(f)| e
i|
y
(f)
; |A
y
(f)| = |A
x
(f)| |A
h
(f)|; |
y
(f) = |
x
(f) |
h
(f).
- Note that correlation is NOT commutative. That is: x(t) h(t) h(t) x(t). This
means that, given x(t) h(t), you must shift h(t).
- Numerical examples of cross-correlation and auto-correlation.
- More Internet resources are on: http://www.answers.com/topic/cross-correlation and
http://www.answers.com/topic/autocorrelation-1.
The z-transform
- The z-transform of a sampled function x(t) = (x
0
, x
1
, , x
N
) is defined as:
N
N
N
n
n
n
z x z x z x x z x z X + + + + = =
=
... ) (
2
2 1 0
0
. (2.15)
- For example, the z-transform of x(t) = (1, -1/2, 2, -1) is X(z) = 1 ()z + 2z
2
z
3
.
- The inverse z-transform is found by collecting the coefficients of the z-transform
polynomial.
- For example, the inverse z-transform of Y(z) = 2 + z - 4z
2
+ z
3
is y(t) = (2, 1, -4, 1).
- The convolution theorem applies to the z-transform: x(t) * h(t) = X(z) H(z).
- Another example of Z-transform.
- More Internet resources are on: http://www.answers.com/topic/z-transform-1.
Phase considerations
- A wavelet is a time-domain signal that has a start and end times.
8
- A minimum-phase wavelet has its energy concentrated at its start time.
- A maximum-phase wavelet has its energy concentrated at its end time.
- A mixed-phase wavelet has its energy concentrated (or distributed) between its start
and end times.
- A zero-phase wavelet is a non-realizable wavelet that starts from negative time and
ends at positive time with its energy symmetric about t = 0. A non-realizable wavelet
is one that cannot be produced directly by real sources.
- Examples.
- The wavelet shape and position with respect to t = 0 can be modified by changing the
phase spectrum. (Example)
- For processing, it is desirable to have a minimum-phase wavelet because
deconvolution requires it; but the wavelet peak will be shifted further in time from the
actual reflection time.
- For interpretation, it is desirable to have a zero-phase wavelet since the wavelet peak
will occur exactly at the actual reflection time.
- It is common to do processing using minimum-phase wavelets and transform them to
zero-phase wavelets before handing the data to the interpreter.
Frequency filtering
- Frequency filtering of a function x(t) means that we apply an operator g(t) to the
function so that the output is a function y(t) that contains only frequencies that lie
within our desired frequency range [f
i
, f
f
].
9
- Filtering can be done in the time-domain by convolving x(t) with g(t) or in the
frequency domain by multiplying their Fourier transforms (the convolution theorem).
- Since we want to change only the amplitude spectrum of x(t), we must select a zero-
phase g(t).
- The desired amplitude spectrum of G(f) is a square function that has a value of one in
the range [f
i
, f
f
], and zero elsewhere.
- Because the FT exists only for continuous functions and G(f) is discontinuous at f
i
and f
f
, we must multiply G(f) by a taper (smoothing) function S(f) to eliminate this
problem. Examples of taper functions include Hanning and Hamming.
- Furthermore, the sharp slopes in G(f) at f
i
and f
f
causes the output y(t) and Y(f) to be
oscillatory because of Gibbs phenomenon. Therefore, relatively gentle slopes are
used at the edges of the square function. (Example)
- With these modifications, the final shape of G(f) will be a trapezoid with smooth
edges.
- The trapezoid G(f) will have a gentle slope from f
i
Af
i
to f
i
, one between f
i
and f
f
, a
gentle slope from f
f
to f
f
+ Af
f
, and zero elsewhere.
- The slope at the high-frequency end should be gentler than that at the low-frequency
end (i. e., Af
f
> Af
i
) because Gibbs phenomenon increases slightly with frequency.
- Shifting the band-pass range [f
i
, f
f
] to higher frequencies will not improve the vertical
resolution of the whole seismic trace. Instead, we have to do this by increasing the
width of the pass-band frequency filter. (Examples 1 and 2)
- Because of absorption of high frequencies by the Earth, a constant pass-band filter
might not be sufficient and a time-variant filter (TVF) might be needed.
10
- TVF means that we segment the seismic record to 3 or 4 segments and design a filter
for each segment separately.
- Real-data and numerical examples.
- More internet resources on tapering functions.
Frequency aliasing
- When sampling a continuous time-dependent signal, we must choose a sampling rate
(sampling interval), At.
- When we sample a function using a sampling rate At, the highest frequency that can
be retrieved is the Nyquist frequency: f
N
= 1/(2At).
- Common sampling rates used in seismic exploration are: 2, 4, and 8 ms. The
corresponding Nyquist frequencies are: 250, 125, and 62.5 Hz, respectively.
(Example)
- Sampling in the time domain means that we multiply the continuous function h(t) by
the time sampling function D(t), which is a series of impulses spaced at At with each
impulse having an amplitude of one.
- Using the convolution theorem, this means that we convolve the amplitude spectra of
h(t) and D(t) namely |H(f)| and D(f), where D(f) is also a series of impulses spaced at
Af = 1/At = 2f
N
with each impulse having an amplitude of 1/At.
- Convolution of |H(f)| and D(f) will replicate |H(f)| at an interval of 2f
N
.
- |H(f)| naturally has a maximum usable frequency f
h
beyond which there is no signal.
- If At is small such that f
N
> f
h
; then, |H(f)| replicas will be totally separated and both
f
N
and f
h
can be retrieved. This is the case with no frequency aliasing. (Examples 1
and 2)
11
- If At is large such that f
N
< f
h
; then, |H(f)| replicas will interfere and neither f
N
nor f
h
can be retrieved. This is the case with frequency aliasing. (Example)
- When frequency aliasing occurs, the maximum retrievable frequency is the alias
frequency given by: f
a
= |2f
N
f
h
|. (Examples Ex2, Ex3)
- For example, if a Vibroseis signal with a frequency range of 10-160 Hz was sampled
at a sampling rate of 4 ms; then:
o f
N
= 1/(2*0.004) = 125 Hz
o f
h
= 160 Hz.
o Since f
N
< f
h
, aliasing will occur.
o The alias frequency is: f
a
= |2125 160| = 90 Hz.
o This means that we lost all the frequencies above 90 Hz. Note that we could not
retrieve neither f
N
nor f
h
because of aliasing.
- On the other hand, if we sampled the same signal at 2 ms sampling interval, f
N
will be
250 Hz while f
h
is still 160 Hz. Since f
N
> f
h
, no aliasing will occur and we can easily
retrieve both f
h
and f
N
.
- Practically, we want to sample the data such that we avoid aliasing. In the field, we
use an anti-aliasing analog filter before sampling to insure that f
N
> f
h
. The anti-
aliasing filter makes f
h
of the analog (continuous) data equal to f
N
/2 or 3f
N
/4.
Gain Applications
- Example.
- Gain is a time-variant scaling of the amplitudes of the data.
- Gain is done to account for spherical divergence and absorption effects.
12
- There are two main types of gain methods:
(1) Physical gain methods, which depend on subsurface velocities.
(2) Display gain methods, which are done purely for display purposes without
requiring velocity information.
(1) Physical Gain Methods:
(a) Single layer:
A seismic source produces a spherical wavefront whose energy decays as 1/r
2
and amplitude as 1/r, where r is the distance from the source.
Therefore, in order to restore the amplitudes to their original values (at the
source), we have to multiply them by the distance traveled by the wave.
Therefore, the gained amplitude is:
G(t) = A(t).V.t, (2.16a)
where A(t) is the amplitude at time t before gain application, V is the layer
velocity, and t is the event time (TWTT).
(b) Multiple layers:
Refraction (ray bending) at layer interfaces increases the decay of amplitudes.
If we know the RMS velocities (V
rms
); then, the gained amplitude is:
G(t) = A(t).V
rms
2
.t. (2.16b)
If we do not know the velocities; then, the gained amplitude is:
G(t) = t
m
, (2.16c)
where m is a number between 1 and 2.
- If absorption effects are considerable, an exponential correction can be used such
that the gained amplitude due to absorption only is:
13
G(t) = A(t).e
(o t)
, (2.17)
where o = 0.01 1 sec
-1
.
- Claerbout (1985) suggests using m = 2 in equation (2.16c) to account for both
spherical divergence and absorption.
(2) Display Gain Methods:
(a) Automatic Gain Control (AGC):
The most common type of AGC is the instantaneous AGC that is performed
as follows:
(1) The absolute mean ( ( )
+
=
=
N k
k i
k
Ai N A / 1 ) of trace samples within a time
window containing samples k to k+N is computed, where k=1,,L-N,
where L is number of samples in the trace.
(2) The trace samples within the time window starting with the k
th
sample and
ending with the (k+N)
th
sample are divided by
k
A .
(3) The time window is shifted down by one sample and steps (1) and (2) are
repeated until the trace end (sample L-N) is reached.
In practice, 200 500 ms AGC time windows are commonly chosen.
This method attempts to scale the amplitudes in the trace to one.
(a) Trace balancing (trace equalization):
Single-trace balancing:
o When the relative amplitudes along a trace is important, then trace
balancing is required for that trace.
14
o Trace balancing is performed by computing the RMS value using the
amplitudes in a single window of the trace and scaling all amplitudes
of the trace by this RMS value.
Relative-trace balancing:
o When the relative amplitudes across a group of traces is important
(such as in AVO studies), then trace balancing is required for that
group of traces (e. g., CMP gather).
o Trace balancing is performed by computing the RMS value using the
amplitudes in a single window of a single trace and scaling all
amplitudes of all traces in the group by this RMS value.
- AGC is a time-variant gain function (i.e., it involves a sliding window) while
trace balancing is time-invariant (i.e., it does not involve a sliding window).
- AGC and trace balancing should only be used to display the data. Whenever
possible, use physical gain methods.
- Example.
15
Appendix A
Numerical convolution, cross-, and auto-correlations
x(t) = (2, 1, -2)
y(t) = (1, -1)
Convolution
f(t) = x(t) * y(t) = (2, -1, -3, 2) = y(t) * x(t) (Prove!)
2 1 -2
-1 1 f(0) = (2)(1) = 2
-1 1 f(1) = (2)(-1) + (1)(1) = -2 + 1 = -1
-1 1 f(2) = (1)(-1) + (-2)(1) = -1 - 2 = -3
-1 1 f(3) = (-2)(-1) = 2
Back
Cross-correlation
f(t) = x(t) y(t) = (-2, 1, 3, -2) = y(t) x(t) (Prove!)
2 1 -2
1 -1 f(-1) = (2)(-1) = -2
1 -1 f(0) = (2)(1) + (1)(-1) = 2 - 1 = 1
1 -1 f(1) = (1)(1) + (-2)(-1) = 1 + 2 = 3
1 -1 f(2) = (-2)(1) = -2
Back
16
Auto-correlation
f(t) = x(t) x(t) = (-4, 0, 9, 0, -4) [Note that: f(-t) = f(t)]
2 1 -2
2 1 -2 f(-2) = (2)(-2) = -4
2 1 -2 f(-1) = (2)(1) + (1)(-2) = 2 - 2 = 0
2 1 -2 f(0) = (2)(2) + (1)(1) + (-2)(-2)= 4 + 1 + 4=9
2 1 -2 f(1) = (1)(2) + (-2)(1) = 2 - 2 = 0
2 1 -2 f(2) = (-2)(2) = -4
Back
Z-Transform
X(z) = (2)z
0
+ (1)z
1
+ (-2)z
2
= 2 + z - 2 z
2
Y(z) = (1)z
0
+ (-1)z
1
= 1 z
Back