Fourier Transform: Applications in seismology
• Estimation of spectra
– windowing – resampling
• Seismograms – frequency content
• Eigenmodes of the Earth
• „Seismo-weather“ with FFTs
• Derivative using FFTs – pseudospectral method for wave propagation
• Convolutional operators: finite-difference operators using FFTs
Scope: Understand how to calculate the spectrum from time series and interpret both phase and amplitude part. Learn other
applications of the FFT in seismology
Fourier: Space and Time
Space
x space variable
L spatial wavelength k=2p/ spatial wavenumber F(k) wavenumber
spectrum
Space
x space variable
L spatial wavelength k=2p/ spatial wavenumber F(k) wavenumber
spectrum
Time
t Time variable
T period
f frequency
=2f angular frequency Time
t Time variable
T period
f frequency
=2f angular frequency Fourier integrals
Fourier integrals
With the complex representation of sinusoidal functions eikx (or (eit) the Fourier transformation can be written as:
With the complex representation of sinusoidal functions eikx (or (eit) the Fourier transformation can be written as:
dx e
x f k
F
dx e
k F x
f
ikx ikx
) 2 (
) 1 (
) 2 (
) 1 (
The Fourier Transform
discrete vs. continuous
dx e
x f k
F
dx e
k F x
f
ikx ikx
) 2 (
) 1 (
) 2 (
) 1 (
1 ,...,
1 , 0 ,
1 ,...,
1 , 0 1 ,
/ 1 2
/ 1 2
0
N k
e F f
N k
e N f
F
N N ikj
j k
N N ikj
j
j k
discrete
continuous
Whatever we do on the computer with data will be based on the discrete Fourier transform
Whatever we do on the computer with data will be based on the discrete Fourier transform
Phase and amplitude spectrum
)
)
(( )
( F e
i F
The spectrum consists of two real-valued functions of angular frequency, the amplitude spectrum mod (F()) and the phase spectrum
In many cases the amplitude spectrum is the most important part to be considered. However there are cases where the phase spectrum plays an important role (-> resonance, seismometer)
Spectral synthesis
The spectrum
Amplitude spectrum Phase spectrum
Fourier space Physical space
The Fast Fourier Transform (FFT)
Most processing tools (e.g. octave, Matlab, Mathematica,
Fortran, etc) have intrinsic functions for FFTs
Most processing tools (e.g. octave, Matlab, Mathematica,
Fortran, etc) have intrinsic functions for FFTs
>> help fft
FFT Discrete Fourier transform.
FFT(X) is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is applied to each column. For N-D arrays, the FFT operation operates on the first non-singleton dimension.
FFT(X,N) is the N-point FFT, padded with zeros if X has less than N points and truncated if it has more.
FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the dimension DIM.
For length N input vector x, the DFT is a length N vector X, with elements
N
X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
n=1
The inverse DFT (computed by IFFT) is given by N
x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
k=1
See also IFFT, FFT2, IFFT2, FFTSHIFT.
>> help fft
FFT Discrete Fourier transform.
FFT(X) is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is applied to each column. For N-D arrays, the FFT operation operates on the first non-singleton dimension.
FFT(X,N) is the N-point FFT, padded with zeros if X has less than N points and truncated if it has more.
FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the dimension DIM.
For length N input vector x, the DFT is a length N vector X, with elements
N
X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
n=1
The inverse DFT (computed by IFFT) is given by N
x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
k=1
See also IFFT, FFT2, IFFT2, FFTSHIFT.
Matlab FFT
Good practice – for estimating spectra
1. Filter the analogue record to avoid aliasing 2. Digitise such that the Nyquist lies above the
highest frequency in the original data 3. Window to appropriate length
4. Detrend (e.g., by removing a best-fitting line) 5. Taper to smooth ends to avoid Gibbs
6. Pad with zeroes (to smooth spectrum or to use
2
npoints for FFT)
Resampling (Decimating)
• Often it is useful to down-sample a time series (e.g., from 100Hz to 1Hz, when looking at
surface waves).
• In this case the time series has to be preprocessed to avoid aliasing effects
• All frequencies above twice the new sampling
interval have to be filtered out
Spectral leakage, windowing, tapering
Care must be taken when extracting time windows when estimating spectra:
– as the FFT assumes periodicity, both ends must have the same value
– this can be achieved by „tapering“
– It is useful to remove drifts as to avoid any discontinuities in the time series -> Gibbs phenonemon
-> practicals
Spectral leakage
Fourier Spectra: Main Cases
random signals
Random signals may contain all frequencies. A spectrum with constant contribution of all frequencies is called a white spectrum
Random signals may contain all frequencies. A spectrum with constant contribution of all frequencies is called a white spectrum
Fourier Spectra: Main Cases
Gaussian signals
The spectrum of a Gaussian function will itself be a Gaussian function. How does the spectrum change, if I make the Gaussian
narrower and narrower?
The spectrum of a Gaussian function will itself be a Gaussian function. How does the spectrum change, if I make the Gaussian
narrower and narrower?
Fourier Spectra: Main Cases
Transient waveform
A transient wave form is a wave form limited in time (or space) in comparison with a harmonic wave form that is infinite
A transient wave form is a wave form limited in time (or space) in comparison with a harmonic wave form that is infinite
Puls-width and Frequency Bandwidth
time (space) spectrum
Narrowing physical signal Widening frequency band
Frequencies in seismograms
Amplitude spectrum
Eigenfrequencies
Sound of an instrument
a‘ - 440Hz
Instrument Earth
26.-29.12.2004 (FFB )
0
S
2– Earth‘s gravest tone T=3233.5s =53.9min
Theoretical eigenfrequencies
The 2009 earthquake „swarm“
Rotations …
First observations of eigenmodes with ring laser!
16 hr time window (Samoa)
36 hr time window (Chile quake)
Time – frequency analysis
Spectral analysis: windowed spectra
24 hour ground motion, do you see any signal?
Seismo-“weather“
Running spectrum of the same data
Shear wave speed variationsGorbatov & Kennett (2003)
?
Western Samoa, 1993/5/16, D=34.97°
Western Samoa, 1993/5/16, D=34.97°
Vz [m/s]Vz [m/s]
t [s]
t [s]
Going beyond ray tomography …
Kristeková et al., 2006 : ... the most complete and informative characterization of a signal can be obtained by its decomposition in the time-frequency plane ... .
u( )g(τ t)e dτ 2π
G(u) 1 :
t) ,
û( i
u ( )g(τ t)e dτ 2π
) 1 G(u :
t) , (
û0 0 0 i
[ t-w representation of synthetics, u(t) ]
[ t-w representation of data, u0(t) ]
Time-frequency representation of seismograms
comparing data with synthetics
… an elegant way of separating phase (i.e., travel time) and amplitude (envelope) information!
Time-frequency misfits
... individually weighted …
Some properties of FT
• FT is linear
signals can be treated as the sum of several signals, the transform will be the sum of their transforms -> stacking is possible!
• FT of a real signals
has symmetry properties
the negative frequencies can be obtained from symmetry
properties
• Shifting corresponds to changing the phase (shift theorem)
• Derivative
) (
* )
( F
F
) ( )
(
) ( )
(
t f e
a F
F e
a t f
a i
a i
