Assignment 6 (15 points)
Considerananalog280-Hzsinusoidalsignal.Ifitwassampledata4-ms interval, what is the alias frequency?
Sketchthespectrumyouwouldexpectfromseismogramsrecordedunder the following circumstances. Label the axes.A noise study reveals a generally white background of−6
StudyLectures17&15 TheamplitudesinQ3haveabroadrange.Theverticalaxisdoesnot
need to draw to scale. See L15 for “How to determine the frequency band of a signal”. Below is another example.
amplitude3x10period range 4-11s. Teleseismic body waves have typical amplitude
, plus microseismic noise of amplitude 10 3×10−4 ms−1 and frequency 0.5-10 Hz; surface waves have typical
Note:
−s
msamplitude 5×10−2 ms−1 and period 10-100 s.
−3 −1msin the
Geophys 6001: Advanced Geophysical Data Analysis
Assignment 6 (15 points)
1. Consider an analog 280-Hz sinusoidal signal. If it was sampled at a 4-ms
interval, what is the alias frequency?
2. At what frequency would a 210 Hz signal be recorded by a digital
recording system with a sampling rate of 130 Hz?
3. Sketch the spectrum you would expect from seismograms recorded under
the following circumstances. Label the axes.
A noise study reveals a generally white background of
amplitude 3 x10 −6 ms − s , plus microseismic noise of amplitude 10 −3 ms −1 in the
period range 4-11s. Teleseismic body waves have typical amplitude
3 x10 −4 ms −1 and frequency 0.5-10 Hz; surface waves have typical
amplitude 5 x10 −2 ms −1 and period 10-100 s.
Note:
Study Lectures 17 & 15
The amplitudes in Q3 have a broad range. The vertical axis does not
need to draw to scale. See L15 for “How to determine the frequency band
of a signal”. Below is another example.
A0
If A0=2
0.707xA0
Then 0.701xA0=1.414
fL
fH
2
GEOPHYS 6001
Advanced Geophysical Data Analysis
Lecture 15
Contents
FFT
Kelly Liu
Digital filter
Spectral leakage, windowing, and taper
Moving average
Procedure to digitize a continuous time series and produce a spectrum
Reference: Gubbins, “Time series analysis and inverse theory for geophysicists”
Liu’s Geophys 6001, Summer 2020
1
FFT
Fast Fourier Transform (FFT)
1
An =
N
N −1
∑ ak e
− 2πink / N
k =0
DFT
n=0,1,2,3,…N-1
The above equation has N items. Each item requires multiplying N complex
sinusoids. The total operations required are therefore O(N2) multiplications.
1 operation
1 sample
A0=a0 x exponential
4 operations 2 samples A0=a0 x exponential + a1 x exponential
A1=a0 x exponential + a1 x exponential
9 operations 3 samples A0=a0 x exponential + a1 x exponential + a2 x exponential
A1=a0 x exponential + a1 x exponential + a2 x exponential
A2=a0 x exponential + a1 x exponential + a2 x exponential
Liu’s Geophys 6001, Summer 2020
2
FFT
Fast Fourier Transform (FFT)
A clever algorithm, called the fast Fourier transform or FFT, accomplished
the same task with substantially less computation.
1
An =
N
N −1
− 2πink / N
a
e
∑ k
k =0
n=0,1,2,3,…N-1
To illustrate the principle of the FFT, suppose N is divisible by 2 and split the sum
into two parts:
N −1
N / 2 −1
N −1
− 2πink / N
− 2πink / N
− 2πink / N
n
k
k
k
k =0
k =0
k=N / 2
NA = ∑ a e
=
N / 2 −1
∑ (a
k
+ ak + N / 2 e
=
−πin
∑a e
)e
− 2πink / N
+
∑a e
=
Divide and conquer algorithm
k =0
1. Forming the combinations (N/2 operations)
2. Performing the DFT the sum of N/2 terms.
3. The total number of operations for all frequency terms is now 2*(N/2)2=N2/2
->computing time has been reduced in half.
Liu’s Geophys 6001, Summer 2020
3
FFT
Fast Fourier Transform (FFT)
Divided
log 2 N
times
Divide and Conquer Algorithm.
Liu’s Geophys 6001, Summer 2020
Smith, “The Scientist and Engineer’s Guide to Digital Signal Processing”
4
FFT
Fast Fourier Transform (FFT)
log 2 N
Cooley–Tukey FFT algorithm, 1965
Divide and conquer algorithm
Liu’s Geophys 6001, Summer 2020
5
FFT
Exercise
Consider a case
N = 10 6 ≈ 2 20
A computer performing a million operations per second.
Calculate the computing time for the conventional DFT method.
Calculate the computing time of the FFT method.
Liu’s Geophys 6001, Summer 2020
6
FFT
The speed advantage of the FFT is enormous!
Liu’s Geophys 6001, Summer 2020
7
FFT
For convolution of two sequences, we transform both sequences with the
FFT, multiple the transforms together, then transform back using the FFT.
It is much faster than doing convolution in time domain.
Liu’s Geophys 6001, Summer 2020
8
Filter
A digital filter processing is defined as a computational process by which an input
time sequence is transformed into a different sequence as an output.
Filtering of a time sequence is convolution with a second, usually shorter,
sequence.
Many important linear processes take the form of a convolution.
For example: the seismometer convolves the ground motion with its own
impulse response, and can therefore be thought of as a filtering
operations.
Liu’s Geophys 6001, Summer 2020
9
Filter
H(f)
Low-pass filters
eliminate all
frequencies
above a certain
value.
Pass band: the range
of frequencies
allowed through.
H(f)
A
H(f)
B
H(f)
D
C
𝜔𝜔c: cut-off frequency
or corner frequency
Or Band Reject Filter or Notch Filter
An ideal filter should have an amplitude spectrum that is equal to zero
outside the pass band and equal to unity inside it.
Liu’s Geophys 6001, Summer 2020
10
Butterworth low-pass filter
1
| F1 (ω ) | =
2n
1 + (ω / ω c )
2
Butterworth filter function;
n = pole
The sharp discontinuity cannot be represented well by a finite number of terms
in the time because of Gibbs’s phenomenon.
A compromise is needed between the sharpness of the cut-off and the length
of the sequence in the time domain.
The compromise is similar to that in tapering windows for spectra.
Tapering windows is to reduce spectral leakage in Fourier spectra.
Liu’s Geophys 6001, Summer 2020
11
Butterworth low-pass filter
1
| F1 (ω ) | =
2n
1 + (ω / ω c )
2
Liu’s Geophys 6001, Summer 2020
12
Filter
Pole=1
Butterworth
Pole=2
Pole=8
Liu’s Geophys 6001, Summer 2020
13
Gubbins “Time Series Analysis and Inverse Theory for Geophysicists”
Filter
High-pass
High-pass filter amplitude spectrum Fh ( w) = 1 − Fl ( w)
Fl (ω ) low-pass filter
Band-pass
| Fb (ω ) | 2 =
1
1 + [(ω − ω b ) / ω c ]
2n
Shifting the lowpass filter function
along the frequency axis to centre ω b
Notch
Fn (ω ) = 1 − Fb (ω )
Liu’s Geophys 6001, Summer 2020
14
Windowing
Ideally our seismogram would be a continuous function of time extending
forever; the spectrum would then be the integral Fourier transform.
In practice, we have a discretized series of finite length, which introduces
some fundamental limitations on the derived spectrum.
Selecting an interval from a longer time sequence is called windowing, and it
can be achieved mathematically by multiplying the ideal sequence by a
boxcar function that is zero everywhere outside the range (0, T) and equal to
unity inside the range.
T
t
Liu’s Geophys 6001, Summer 2020
15
Gubbins “Time Series Analysis and Inverse Theory for Geophysicists”
Spectral leakage
In time-domain, we multiple
the time sequence with a
Boxcar function.
In frequency domain, the
boxcar function contains a
peak with many spurious
secondary peaks. This is
called spectral leakage.
side lobes
Liu’s Geophys 6001, Summer 2020
16
Spectral leakage & Taper
Corresponding to a boxcar window
The peaks have been
broadened by the finite length
of record.
Prominent side lobes appear
as the fuzzy in-fill between the
major peaks; they overlap
throughout the low frequency
(long period) part of the
spectrum.
Ideal, calculated theoretically
would be obtained from a record of infinite length
A synthetic normal mode spectrum from an
explosive source at a distance of 950.
Figure. 3.3. Synthetic boxcar window power spectra
for an explosive source recorded at a distance of 950.
This is a fundamental limitation placed on the data by the limited recording time.
Liu’s Geophys 6001, Summer 2020
17
Gubbins “Time Series Analysis and Inverse Theory for Geophysicists”
Spectral leakage & Taper
To reduce the side lobes, we can taper the time series using different window
functions, rather than the Boxcar window.
Many tapering functions
Bartlett window (triangle)
Gaussian window
Hanning window
Liu’s Geophys 6001, Summer 2020
18
Spectral Leakage & Taper
c
Effect of Hanning and Boxcar
b
a
Ideal
Figure. 3.3. Synthetic boxcar window power spectra
for an explosive source recorded at a distance of 950.
Side lobes are reduced by smoothing the sharp edges of
Liu’s Geophys 6001, Summer 2020
the rectangular boxcar window. This
is called tapering.
19
Gubbins “Time Series Analysis and Inverse Theory for Geophysicists”
Filter
aftershock
A high-pass filter does a good job of removing surface waves from
the record, leaving the body waves from the aftershock unchanged.
Liu’s Geophys 6001, Summer 2020
20
Gubbins, “Time series analysis and inverse theory for geophysicists
Frequency band
How to determine the frequency band of a signal?
Apeak
√2/2 (~0.707) Apeak
Liu’s Geophys 6001, Summer 2020
21
Seismic spectrum
Liu’s Geophys 6001, Summer 2020
22
Kearey, Brooks, and Hill, “An Introduction to Geophysical Exploration”
Moving average
x1 + x2 + x3 + x4 + x5 x2 + x3 + x4 + x5 + x6
,……,
,
X5 = (
5
5
xn−4 + xn−3 + xn−2 + xn−1 + xn
)
5
Moving average for 5 samples
In stock market, moving average can smooth out the “noise” of shorter-term price
fluctuations, for better identifying and define significant underlying trends.
Liu’s Geophys 6001, Summer 2020
23
Moving average
Convolution of the original time
sequence with a boxcar function
Filtering operation
Liu’s Geophys 6001, Summer 2020
24
Smith, “The Scientist and Engineer’s Guide to Digital Signal Processing”
Procedure
To digitize a continuous time series and produce a spectrum
Filter the analogue record to remove all energy near and above the Nyquist
frequency. The signal is now band-limited. will learn in future
Digitize using a sampling interval such that the Nyquist frequency lies above the
highest frequency in the original data. The sampling theorem ensures that no
information is lost.
Window the data with a finite-length window. In practice this may be imposed by the
duration of the experiment, or you may be choosing a suitable piece of a continuous
recording.
Detrend (e.g., by removing a best-fitting straight line). Show variations across the
entire time sequence contain significant energy at all frequencies and can mask the
high frequency oscillations of interest.
Taper to smooth the ends of the time series to zero. This reduces spectral leakage
but broadens any central peaks. It also has the effect of equalizing end points of the
data to zero, thus avoiding any discontinuities seen by the DFT from periodic
repetition. Tapering is equivalent to convolution in the frequency domain; it smoothes
the spectrum.
Pad with zeroes, usually to make the number of samples a power of 2 for the FFT.
This reduces the spacing in the frequency domain and smoothes the spectrum by
interpolation.
Liu’s Geophys 6001, Summer 2020
25
GEOPHYS 6001
Advanced Geophysical Data Analysis
Lecture 17
Contents
Sampling theorem
Frequency aliasing
Nyquist frequency
Anti-alias filter
Kelly Liu
Reference: Yilmaz, “Seismic data analysis”
Liu’s Geophys 6001, Summer 2020
1
Sampling Theorem
To digitize a signal, we take
samples of the analog signal
every certain second.
Apparently the reconstructed
signal lacks the details
presented in the original
analog signal. The high
frequency components were
lost by sampling.
Liu’s Geophys 6001, Summer 2020
2
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Sampling period (or sampling rate Δt): the time interval between
successive samples
Sampling frequency (Hz): the reciprocal of sampling period (1/Δt=Δf)
SPS = samples per second
Given an analog signal, how should we select the sampling rate Δt , or
the sampling frequency Δ f ?
Know major information about the characteristics of the signal to be sampled.
know the maximum frequency content of the general class of signals, we can specify
the sampling rate necessary to convert the analog signals to digital signals.
Liu’s Geophys 6001, Summer 2020
3
Sampling Theorem
∆f > 2 f max
f N = 2 f max or
fN
1
fN =
2∆t
Nyquist frequency
is half of the
sampling rate.
—-Nyquist frequency
Need to know the frequency content of the signal
Liu’s Geophys 6001, Summer 2020
4
Sampling Theorem
Assume that any analog signal can be represented as a sum of sinusoids
of different amplitudes, frequencies, and phases.
N
x(t ) = ∑ Ai cos(2πf i t + θ i )
Exercise
i =1
x (t ) = 3 cos 50πt + 10 sin 300πt − cos 100πt
Q: Nyquist frequency for this signal?
Liu’s Geophys 6001, Summer 2020
5
Sampling Theorem
N
x(t ) = ∑ Ai cos(2πf i t + θ i )
i =1
Exercise Solution
x (t ) = 3 cos 50πt + 10 sin 300πt − cos 100πt
Q: Nyquist frequency for this signal?
The frequencies present in the signal above are
f1 = 25 Hz
Thus
f 2 = 150 Hz
f 3 = 50 Hz
f max = 150 Hz
Therefore the Nyquist frequency is
and ∆f > 2 f max = 300 Hz
f N = 2 f max = 300 Hz
The sampling frequency must be greater than 300 Hz.
Liu’s Geophys 6001, Summer 2020
6
Sampling Theorem
Discussion
x (t ) = 3 cos 50πt + 10 sin 300πt − cos 100πt
It should be observed that the signal component
If we sample it at the Nyquist rate 300 Hz,
It results in the samples
10 sin 300πt
f N = 2 f max = 300 Hz
10 sin πt
which is identically zero. In other words, we are sampling the analog sinusoid at its
zero-crossing points, and hence we miss this signal component completely. This
situation does not occur if the sine function is offset in phase by some amount.
A simple way to avoid this potentially troublesome situation is to sample the
analog signal at a rate higher than the Nyquist rate.
Liu’s Geophys 6001, Summer 2020
7
Sampling Theorem
1 ms = 1 millisecond = 10-3 s
Exercise
fN = ?
fN = ?
fN = ?
1
fN =
2∆t
The original time
sequence was resampled with 4- and 8-ms
sampling intervals.
What are the
corresponding Nyquist
frequencies?
The coarser the sampling
interval, the smoother the
series. What make them
smoother? Loss of high
frequencies!
Liu’s Geophys 6001, Summer 2020
8
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
f N =250
? Hz
f N =125
? Hz
125 and
250 Hz are
missing!
f N =62.5
? Hz
62.5 and
250 Hz are
missing!
Liu’s Geophys 6001, Summer 2020
9
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
In the last panel, frequency components between 62.5 and 250 Hz are missing
from the series re-sampled to 8 ms.
Can these frequencies be recovered? No! Once a continuous signal is digitized,
the high frequency can not be restored.
We can interpolate 4 ms and 8 ms signal to 2 ms sampling interval. But
interpolation does not recover the frequencies lost by sampling; it only generates
extra samples.
Liu’s Geophys 6001, Summer 2020
10
Yilmaz, “Seismic Data Analysis”
Frequency Aliasing
f a =| 2 f N − f s |
fa
fN
—-Alias frequency
—-Nyquist frequency
f s —Signal frequency
Liu’s Geophys 6001, Summer 2020
11
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
A 25 Hz sine function. This
signal is re-sampled to 4
and 8 ms.
The signals didn’t change
after re-sampling it to a
coarser sampling interval.
fN=62.5 Hz
Liu’s Geophys 6001, Summer 2020
12
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
fN=62.5 Hz
No frequency aliasing
Liu’s Geophys 6001, Summer 2020
13
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
fN=250 Hz
A 75 Hz frequency.
fN=125 Hz
The signal is sampled as 2
originally.
The signal looks much smoother
after resampled using a coarser
sampling interval, 8 ms.
fN=62.5 Hz
f a =| 2 f − f s |=| 2 x62.5 − 75 |= 50 Hz
Liu’s Geophys 6001, Summer 2020
14
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
The Nyquist
frequency for
an 8 ms
sampling
interval is 62.5
Hz. The true
signal
frequency is 75
Hz.
fN=125 Hz
As a result of
resampling, the
signal with 75
Hz frequency
folded back
onto the
spectrum and
appeared at its
alias frequency
of 50 Hz.
fN=62.5 Hz
folded back
Liu’s Geophys 6001, Summer 2020
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
A 150 Hz frequency.
The signal is sampled as 2
originally.
The signal looks much smoother
after resampled using a coarser
sampling intervals 4 & 8 ms.
f a =| 2 f − f s |=| 2 x125 − 150 |= 100 Hz
fN=62.5 Hz
f a =| 2 f − f s |=| 2 x62.5 − 150 |= 25 Hz
Liu’s Geophys 6001, Summer 2020
16
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
fN=62.5 Hz
Liu’s Geophys 6001, Summer 2020
17
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
fN=62.5 Hz
By testing using sine function, you can see the frequencies above the Nyquist frequency
Geophysbelow
6001,the
Summer
really are not lost after sampling, but appear Liu’s
at frequencies
Nyquist.
2020
18
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
fN=62.5 Hz
Superposition of two sine
functions with frequencies of
12.5 and 75 Hz
Digitized at 8 ms, the 12.5 Hz
component is not affected, because
8 ms sampling still is sufficient to
sample this low-frequency
component.
75 Hz shows as 50 Hz.
Liu’s Geophys 6001, Summer 2020
19
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
fN=250 Hz
fN=125 Hz
fN=62.5 Hz
Folded back
Liu’s Geophys 6001, Summer 2020
20
Yilmaz, “Seismic Data Analysis”
Sampling Theorem
Exercise
1
fN =
2∆t
∆t = 4ms
f s = 75 Hz
fN=? fN=125 Hz
fa=?
f s = 175 Hz
fa=|2fN-fs|=|2×125-175|=75 Hz
f s = 250 Hz
fa=|2fN-fs|=|2×125-250|=0 Hz
Liu’s Geophys 6001, Summer 2020
21
Sheriff and Geldart, “Exploration Seismology”
Sampling Theorem
Exercise
∆t = 2ms
fs=500 Hz
fN=250 Hz
Liu’s Geophys 6001, Summer 2020
22
Gubbins “Time Series Analysis and Inverse Theory for Geophysicists”
Sampling Theorem
∆t = 2ms
fs=500 Hz
fN=250 Hz
fa=|2fN-fs|=|2×250-417|=83 Hz Liu’s Geophys 6001, Summer 2020
23
Gubbins “Time Series Analysis and Inverse Theory for Geophysicists”
Sampling Theorem
Undersampling has two effects
Band limiting the spectrum of the continuous signal, with the
maximum frequency being the Nyquist frequency.
Contamination of the digital signal spectrum by high frequencies
beyond the Nyquist, which may have been present in the continuous
signal. Very high frequencies can enter at very low frequency and
completely mask the signal of interest.
Liu’s Geophys 6001, Summer 2020
24
Sampling Theorem
Nothing can be done about the first problem.
To keep the recoverable frequency back between zero and the Nyquist
frequency free from aliased frequencies, a high-cut anti-aliasing filter
can be applied before analog-to-digital conversion of signals.
This filter eliminates those frequency components that would have
been aliased during sampling. This is usually done by an anti-alias
filter applied to the analogue signal.
—–Anti-alias filter
Liu’s Geophys 6001, Summer 2020
25
Procedure
To digitize a continuous time series and produce a spectrum
Filter the analogue record to remove all energy near and above the Nyquist frequency. The
signal is now band-limited.
Digitize using a sampling interval such that the Nyquist frequency lies above the highest
frequency in the original data. The sampling theorem ensures that no information is lost.
Window the data with a finite-length window.
Detrend (e.g., by removing a best-fitting straight line). Show variations across the entire time
sequence contain significant energy at all frequencies and can mask the high frequency
oscillations of interest.
Taper to smooth the ends of the time series to zero. This reduces spectral leakage but
broadens any central peaks. Tapering is equivalent to convolution in the frequency domain; it
smoothes the spectrum.
Pad with zeroes usually to make the number of samples a power of 2 for the FFT. This
reduces the spacing in the frequency domain and smoothes the spectrum by interpolation.
Liu’s Geophys 6001, Summer 2020
26
Summary
Sampling theorem
Frequency aliasing
Nyquist frequency
Anti-alias filter
Liu’s Geophys 6001, Summer 2020
27