Transcript Fourier Transform: Applications - uni

Fourier Transform:
Applications
• Seismograms
• Eigenmodes of the Earth
• Time derivatives of
seismograms
• The pseudo-spectral
method for acoustic wave
propagation
Fourier: Applications
Modern Seismology – Data processing and inversion
1
Fourier: Space and Time
Space
x
space variable
L
spatial wavelength
k=2p/l spatial wavenumber
F(k)
wavenumber spectrum
t
T
f
w=2pf
Time
Time variable
period
frequency
angular frequency
Fourier integrals
With the complex representation of sinusoidal functions eikx (or
(eiwt) the Fourier transformation can be written as:
1
f ( x) 
2p
1
F (k ) 
2p
Fourier: Applications

ikx
F
(
k
)
e
dx




f ( x)eikxdx

Modern Seismology – Data processing and inversion
2
The Fourier Transform
discrete vs. continuous
Whatever we do on the
computer with data will
be based on the discrete
Fourier transform
discrete
1
Fk 
N
N 1

j 0
continuous
1
f ( x) 
2p
1
F (k ) 
2p

ikx
F
(
k
)
e
dx




f ( x)eikxdx

f j e  2pikj / N , k  0,1,..., N  1
N 1
f k   F j e 2pikj / N , k  0,1,..., N  1
j 0
Fourier: Applications
Modern Seismology – Data processing and inversion
3
The Fast Fourier Transform
... the latter approach became interesting with the introduction of the
Fast Fourier Transform (FFT). What’s so fast about it ?
The FFT originates from a paper by Cooley and Tukey (1965, Math.
Comp. vol 19 297-301) which revolutionised all fields where Fourier
transforms where essential to progress.
The discrete Fourier Transform can be written as
1
uˆk 
N
N 1
 2pikj / N
u
e
, k  0,1,..., N  1
 j
j 0
N 1
uk   uˆ j e 2pikj / N , k  0,1,..., N  1
j 0
Fourier: Applications
Modern Seismology – Data processing and inversion
4
The Fast Fourier Transform
... this can be written as matrix-vector products ...
for example the inverse transform yields ...
1
1
1
w

1 w 2






N 1
1 w
1
1
w2
w4
w3
w6


  uˆ0   u0 
 w N 1   uˆ1   u1 
 w 2 N  2   uˆ 2   u2 



      
      
 

( N 1) 2  
 w
 uˆ N 1  u N 1 

1
.. where ...
w  e 2pi / N
Fourier: Applications
Modern Seismology – Data processing and inversion
5
The Fast Fourier Transform
... the FAST bit is recognising that the full matrix - vector multiplication
can be written as a few sparse matrix - vector multiplications
(for details see for example Bracewell, the Fourier Transform and its
applications, MacGraw-Hill) with the effect that:
Number of multiplications
full matrix
FFT
N2
2Nlog2N
this has enormous implications for large scale problems.
Note: the factorisation becomes particularly simple and effective
when N is a highly composite number (power of 2).
Fourier: Applications
Modern Seismology – Data processing and inversion
6
The Fast Fourier Transform
Number of multiplications
Problem
1D (nx=512)
1D (nx=2096)
1D (nx=8384)
full matrix
2.6x105
FFT
Ratio full/FFT
9.2x103
28.4
94.98
312.6
.. the right column can be regarded as the speedup of an algorithm
when the FFT is used instead of the full system.
Fourier: Applications
Modern Seismology – Data processing and inversion
7
Spectral synthesis
The red trace is the sum of all blue traces!
Fourier: Applications
Modern Seismology – Data processing and inversion
8
Phase and amplitude spectrum
The spectrum consists of two real-valued functions of angular
frequency, the amplitude spectrum mod (F(w)) and the phase
spectrum f(w)
i (w )
F (w)  F (w) e
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)
Fourier: Applications
Modern Seismology – Data processing and inversion
9
… remember …
z*  a  ib  r (cosf  i sin f )
 r cos f  ri sin(f )  r
 if
z  zz*  (a  ib)(a  ib)  r
2
Fourier: Applications
Modern Seismology – Data processing and inversion
2
10
The spectrum
Phase spectrum
Fourier space
Amplitude spectrum
Fourier: Applications
Modern Seismology – Data processing and inversion
11
The Fast Fourier Transform
(FFT)
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.
Matlab FFT
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.
Fourier: Applications
Modern Seismology – Data processing and inversion
12
Frequencies in seismograms
Fourier: Applications
Modern Seismology – Data processing and inversion
13
Amplitude spectrum
Eigenfrequencies
Fourier: Applications
Modern Seismology – Data processing and inversion
14
Sound of an instrument
a‘ - 440Hz
Fourier: Applications
Modern Seismology – Data processing and inversion
15
Instrument Earth
26.-29.12.2004 (FFB )
0S2 –
Earth‘s gravest tone
T=3233.5s =53.9min
Theoretical eigenfrequencies
Fourier: Applications
Modern Seismology – Data processing and inversion
16
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
Fourier: Applications
Modern Seismology – Data processing and inversion
17
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?
Fourier: Applications
Modern Seismology – Data processing and inversion
18
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
Fourier: Applications
Modern Seismology – Data processing and inversion
19
Puls-width and Frequency Bandwidth
spectrum
Narrowing physical signal
Widening frequency band
time (space)
Fourier: Applications
Modern Seismology – Data processing and inversion
20
Spectral analysis: an Example
24 hour ground motion, do you see any signal?
Fourier: Applications
Modern Seismology – Data processing and inversion
21
Seismo-Weather
Running spectrum of the same data
Fourier: Applications
Modern Seismology – Data processing and inversion
22
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
• FT of a real signals
has symmetry properties
F (w )  F * (w )
the negative frequencies can be
obtained from symmetry
properties
Fourier: Applications
• Shifting corresponds to
changing the phase (shift
theorem)
f (t  a)  e iwa F (w )
F (w  a)  e iwa f (t )
• Derivative
dn
f (t )  (iw ) n F (w )
dt
Modern Seismology – Data processing and inversion
23
Fourier Derivatives
.. let us recall the definition of the derivative using Fourier integrals ...


ikx
 x f ( x)   x   F (k )e dk 
 


   ikF (k )e ikxdk

... we could either ...
1) perform this calculation in the space domain by convolution
2) actually transform the function f(x) in the k-domain and back
Fourier: Applications
Modern Seismology – Data processing and inversion
24
Acoustic Wave Equation - Fourier Method
let us take the acoustic wave equation with variable density
1

1
2

t p   x   x p 
2

c


the left hand side will be expressed with our
standard centered finite-difference approach
1

1

 p(t  dt)  2 p(t )  p(t  dt)   x   x p 
2
2
c dt


... leading to the extrapolation scheme ...
Fourier: Applications
Modern Seismology – Data processing and inversion
25
Acoustic Wave Equation - Fourier Method
1


p(t  dt)  c dt  x   x p 
 2 p(t )  p(t  dt)



2
2
where the space derivatives will be calculated using the Fourier Method.
The highlighted term will be calculated as follows:
ˆ n  ik P
ˆ n  FFT1   P n
Pjn  FFT  P

 
x j
multiply by 1/
n
n
 1 ˆ
 1 ˆ
1
1
n
1
n
 x Pj  FFT    x P   ik   x P   FFT   x   x Pj 







... then extrapolate ...
Fourier: Applications
Modern Seismology – Data processing and inversion
26
... and the first derivative using FFTs ...
function
df=
%
SDER1D(f,d
nx=max(size(
%
initialize
kmax=pi/dx;
dk=kmax/(nx/
for
i=1:nx/2
k=sqrt(-1)*k
%
FFT
and
IF
ff=fft(f);
f
.. simple and elegant ...
Fourier: Applications
Modern Seismology – Data processing and inversion
27
Fourier Method - Comparison with FD - Table
Difference (%) between numerical and analytical solution
as a function of propagating frequency
160
140
120
Simulation time
5.4s
3 point
7.8s
5 point
33.0s
Fourier
100
80
60
40
20
0
5 Hz
Fourier: Applications
10 Hz
20 Hz
Modern Seismology – Data processing and inversion
28
Numerical solutions and Green’s Functions
10
10
9
9
9
8
8
8
7
7
7
6
6
6
5
5
5
4
4
4
3
3
3
2
2
2
1
1
1
1000
1500
0
500
1000
0
1500 500
1000
Modern Seismology – Data processing and inversion
1500
Impulse response (numerical convolved with source
Fourier Method
10
0
500
Fourier: Applications
5 point operator
Impulse response (analytical) concolved with source
Frequency increases
3 point operator
29