Fourier Transform:

(1)

Fourier Transform:

Applications

• Seismograms

• Eigenmodes of the Earth

• Time derivatives of seismograms

• The pseudo-spectral

method for acoustic wave propagation

(2)

Modern Seismology – Data processing and inversion

2

Fourier: Applications

Fourier: Space and Time

Space

x space variable

L spatial wavelength k=2π/λ spatial wavenumber F(k) wavenumber spectrum

Space

x space variable

L spatial wavelength k=2π/λ spatial wavenumber F(k) wavenumber spectrum

t Time variableTime

T period

f frequency

ω=2πf angular frequency t Time variableTime

T period

f frequency

ω=2πf angular frequency Fourier integrals

Fourier integrals

With the complex representation of sinusoidal functions e^ikx (or (e^iwt) 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 (

π π

(3)

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 ,

/ 2 1

0

/ 2 1

0

−

=

−

=

∑

−

=

− −

=

N k

e F f

N k

e N f

F

N ikj N

j

j k

N ikj N

j

j k

π

discrete

continuous

Whatever we do on the computer with data will be based on the discrete Fourier transform

(4)

4

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 ,...,

1 , 0 ˆ ,

1 ,...,

1 , 0 1 ,

ˆ

/ 2 1

0

/ 2 1

0

−

=

−

=

∑

−

=

− −

=

N k

e u u

N k

e N u

u

N ikj N

j

j k

N ikj N

j

j k

π

(5)

... this can be written as matrix-vector products ...

for example the inverse transform yields ...

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

=

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

−

− −

−

1 2 1 0

) 1 ( 1

2 2 6

4 2

1 3

2

ˆ ˆ

1 1 1

1 1

1

2

N N

u u

M M M

M L

L L

M M

M

M M

M

K K K

ω ω

.. where ...

N

e²^πi ^/

ω =

(6)

6

... 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 Number of multiplications

full matrix FFT

N² 2Nlog₂N

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).

(7)

.. the right column can be regarded as the speedup of an algorithm when the FFT is used instead of the full system.

Number of multiplications Number of multiplications

Problem full matrix FFT Ratio full/FFT 1D (nx=512) 2.6x10⁵ 9.2x10³ 28.4

1D (nx=2096) 94.98 1D (nx=8384) 312.6

(8)

8

Spectral synthesis

The red trace is the sum of all blue traces!

(9)

Phase and amplitude spectrum

)

(

( )

( ω = F ω e

ⁱ^Φ ^ω

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)

(10)

10

… remember …

2

* ( )( )

) sin(

cos

) sin

(cos

*

r ib

a ib

a zz

z

r ri

r

i r

ib a

z

i

=

− +

=

−

=

−

=

−

=

− φ

φ φ

φ

(11)

The spectrum

Amplitude spectrum

Amplitude spectrum Phase spectrumPhase spectrum

FourierspaceFourierspace PhysicalspacePhysicalspace

(12)

12

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

X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.N The inverse DFT (computed by IFFT) is given byn=1

x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.N k=1

Frequencies in seismograms

(14)

14

Amplitude spectrum

Eigenfrequencies

(15)

Sound of an instrument

a‘ - 440Hz

(16)

16

Instrument Earth

26.-29.12.2004 (FFB )

0 S₂– Earth‘s gravest tone T=3233.5s =53.9min

Theoretical eigenfrequencies

(17)

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

(18)

18

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?

(19)

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

(20)

20

Puls-width and Frequency Bandwidth

time (space) spectrum

Narrowingphysicalsignal Wideningfrequencyband

(21)

Spectral analysis: an Example

24 hour ground motion, do you see any signal?

(22)

22

Seismo-Weather

Running spectrum of the same data

(23)

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

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

ω ω

ω

−

→

−

→

−

) ( )

( )

(t iω F ω

dt f

dⁿ → _n

(24)

24

Fourier Derivatives

∫

∞

−

∞

−

=

⎟⎟⎠

⎞

⎜⎜⎝

∂ ⎛

=

∂

dk e

k ikF

dk e

k F x

f

ikx ikx x

x

) (

) ( )

(

.. let us recall the definition of the derivative using Fourier integrals ...

... 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

(25)

Acoustic Wave Equation - Fourier Method

let us take the acoustic wave equation with variable density

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ ∂

∂

=

∂ p p

c ^t ^x ρ ^x

ρ

1

1 ₂

2

the left hand side will be expressed with our standard centered finite-difference approach

[ ] _⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛ ∂

∂

=

− +

−

+ dt p t p t dt p

t dt p

c ^x ρ ^x

ρ

) 1 (

) ( 2 )

1 (

2 2

... leading to the extrapolation scheme ...

(26)

26

Acoustic Wave Equation - Fourier Method

where the space derivatives will be calculated using the Fourier Method.

The highlighted term will be calculated as follows:

) (

) ( 1 2

)

(t dt c ² dt ² p p t p t dt

p _x _x ⎟⎟ + − −

⎠

⎜⎜ ⎞

⎝

⎛ ∂

∂

=

+ ρ ρ

n j x n

n n

j P ik P P

P → FFT → ˆ_υ → _υ ˆ_υ → FFT ⁻¹ → ∂

multiply by 1/ρ

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ ∂

∂

→

⎟⎟ →

⎠

⎜⎜ ⎞

⎝

⎛ ∂

⎟⎟ →

⎠

⎜⎜ ⎞

⎝

⎛ ∂

→

∂ ⁻ _x _x _jⁿ

n x

n j

xP P ik P P

ρ ρ

ρ

ρ _υ ^υ _υ

FFT 1 1 ˆ

1 ˆ

1 FFT ₁

... then extrapolate ...

(27)

... and the first derivative using FFTs ...

function df=sder1d(f,dx)

% SDER1D(f,dx) spectral derivative of vector nx=max(size(f));

% initialize k kmax=pi/dx;

dk=kmax/(nx/2);

for i=1:nx/2, k(i)=(i)*dk; k(nx/2+i)=-kmax+(i)*dk; end k=sqrt(-1)*k;

% FFT and IFFT

ff=fft(f); ff=k.*ff; df=real(ifft(ff));

.. simple and elegant ...

(28)

28

Fourier Method - Comparison with FD - Table

0 20 40 60 80 100 120 140 160

5 Hz 10 Hz 20 Hz

3 point 5 point Fourier

Difference (%) between numerical and analytical solution as a function of propagating frequency

Simulation time 5.4s

7.8s 33.0s

(29)

500 1000 1500 0

1 2 3 4 5 6 7 8 9 10

500 100 0

1 2 3 4 5 6 7 8 9 10

500 1000 1500 0

1 2 3 4 5 6 7 8 9 10

Numerical solutions and Green’s Functions

3 point operator 5 point operator Fourier Method

Frequency increases Impulse response (analytical) concolvedwith source Impulse response (numerical convolved with source