Filtering Geophysical Data: Be careful!

(1)

1

 Filtering: basic concepts

 Seismogram examples, high-low-bandpass filters

 The crux with causality

 Windowing seismic signals

 Various window functions

 Multitaper approach

 Wavelets (principle)

Scope: Understand the effects of filtering on time series (seismograms). Get to know frequently used windowing functions.

(2)

Why filtering

1. Get rid of unwanted frequencies

2. Highlight signals of certain frequencies 3. Identify harmonic signals in the data

4. Correcting for phase or amplitude characteristics of instruments

5. Prepare for down-sampling 6. Avoid aliasing effects

(3)

3

A seismogram

Time (s)

Frequency (Hz) Amplitude Spectral amplitude

(4)

Digital Filtering

Often a recorded signal contains a lot of information that we are not interested in (noise). To get rid of this noise we can apply a filter in the frequency domain.

The most important filters are:

• High pass: cuts out low frequencies

• Low pass: cuts out high frequencies

• Band pass: cuts out both high and low frequencies and leaves a band of frequencies

• Band reject: cuts out certain frequency band and leaves all other frequencies

(5)

5

Cutoff frequency

(6)

Cut-off and slopes in spectra

(7)

7

Digital Filtering

(8)

Low-pass filtering

(9)

9

Lowpass filtering

(10)

High-pass filter

(11)

11

Band-pass filter

(12)

The simplemost filter

The simplemost filter gets rid of all frequencies above a certain cut-off frequency (low-pass), „boxcar“











 

0 0

) 1

(  



 

if H_L if

-1000 -50 0 50 100

0.2 0.4 0.6 0.8 1

Frequency (Hz)

Filter Amplitude

(13)

13

The simplemost filter

… and its brother

… (high-pass)











 





0 0

) 1 ( 1

)

(  



 

 if

H if

H_H _L

-1000 -50 0 50 100

0.2 0.4 0.6 0.8 1

Frequency (Hz)

Filter Amplitude

(14)

… let‘s look at the consequencse







 F k e dt

t

f ^ik^

 ⁽ ⁾ 2

) 1 (

-1000 -50 0 50 100

0.2 0.4 0.6 0.8 1

Frequency (Hz)

Filter Amplitude







 H F e dt

t

f ^filt   ^ik^

 ⁽ ⁾ ⁽ ⁾ 2

) 1 (

) ( H

… but what does H() look like in the time domain … remember the

convolution theorem?

(15)

15

… surprise …







 H e dt

t

h  ^ik^

 ⁽ ⁾

2 ) 1

(

-1000 -50 0 50 100

0.2 0.4 0.6 0.8 1

Frequency (Hz)

Filter Amplitude )(H

-40 -30 -20 -10 0 10 20 30 40

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2

(16)

Zero phase and causal filters

Zero phase filters can be realised by

 Convolve first with a chosen filter

 Time reverse the original filter and convolve again

 First operation multiplies by F(), the 2nd operation is a multiplication by F*()

 The net multiplication is thus | F(w)|²

 These are also called two-pass filters

(17)

17

The Butterworth Filter (Low-pass, 0-phase)

n c

FL ₂

) /

( 1

) 1

(  

 

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

1 Butterworth n=1, f0=20 Hz

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

(18)

… effect on a spike …

0 1 2 3 4 5

0 0.2 0.4 0.6 0.8

1 Original function

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0.05 0.1 0.15 0.2

Filtered with n=1, f0=20 Hz

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15 0.2

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15

0.2 Filtered with n=9, f0=20 Hz

Time (s)

Amplitude

(19)

19

… on a seismogram …

… varying the order …

40 45 50 55 60 65

-5 0 5

x 10^-6 Original function

Time (s)

Amplitude

40 45 50 55 60 65

-1 0 1

x 10^-6 Filtered with n=4, f0=1 Hz

Time (s)

Amplitude

40 45 50 55 60 65

-2 -1 0 1

2x 10^-6 Filtered with n=9, f0=1 Hz

Time (s)

Amplitude

40 45 50 55 60 65

-2 -1 0 1 2

Time (s)

Amplitude

(20)

… varying the cut-off frequency…

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

40 45 50 55 60 65

-1 0 1

Time (s)

Amplitude

40 45 50 55 60 65

-1 -0.5 0 0.5 1

x 10^-6 Filtered with n=4, f0=0.666667 Hz

Time (s)

Amplitude

40 45 50 55 60 65

-1 -0.5 0 0.5 1

Time (s)

Amplitude

(21)

21

The Butterworth Filter (High-Pass)

n c

FH ₂

) /

( 1 1 1

)

(  

 



0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

(22)

0 1 2 3 4 5

0 0.2 0.4 0.6 0.8

1 Original function

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6

Time (s)

Amplitude

(23)

23

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

(24)

50 52 54 56 58 60

-5 0 5

Time (s)

Amplitude

50 52 54 56 58 60

-4 -2 0 2 4

Time (s)

Amplitude

50 52 54 56 58 60

-4 -2 0 2

4x 10^-6 Filtered with n=4, f0=1.5 Hz

Time (s)

Amplitude

50 52 54 56 58 60

-4 -2 0 2

4x 10^-6 Filtered with n=4, f0=2 Hz

Time (s)

Amplitude

(25)

25

The Butterworth Filter (Band-Pass)

 _b  ⁿ

FBP ₂

/ ) (

1 1 1 )

(   





 



0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8

Frequency (Hz)

Filter amplitude

Hz

 5



(26)

0 1 2 3 4 5

0 0.2 0.4 0.6 0.8

1 Original function

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

-0.05 0 0.05 0.1

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

-0.05 0 0.05

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

-0.05 0 0.05

Time (s)

Amplitude

(27)

27

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

40 45 50 55 60 65

-2 -1 0 1 2

Time (s)

Amplitude

40 45 50 55 60 65

-2 -1 0 1 2

Time (s)

Amplitude

40 45 50 55 60 65

-2 0 2

Time (s)

Amplitude

(28)

50 52 54 56 58 60

-5 0 5

Time (s)

Amplitude

50 52 54 56 58 60

-5 0 5

Time (s)

Amplitude

50 52 54 56 58 60

-5 0 5

Time (s)

Amplitude

50 52 54 56 58 60

-5 0 5

Time (s)

Amplitude

(29)

29

Zero phase and causal filters

0 1 2 3 4 5

0 0.2 0.4 0.6 0.8

1 Original function

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0.05 0.1 0.15 0.2

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15 0.2

Time (s)

Amplitude

0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15

Time (s)

Amplitude

When the phase of a filter is set to zero (and simply the amplitude spectrum is inverted) we obtain a zero-phase filter. It means a peak will not be shifted.

Such a filter is acausal. Why?

(30)

Butterworth Low-pass (20 Hz) on spike

0 50 100

0 0.5

Butterworth n=1, f0=20 Hz 1

Frequency (Hz)

Filter amplitude

0 50 100

0 0.5

Frequency (Hz)

Filter amplitude

0 50 100

0 0.5

Frequency (Hz)

Filter amplitude

0 50 100

0 0.5

Frequency (Hz)

Filter amplitude

(31)

31

(causal) Butterworth Low-pass (20 Hz) on spike

0 5

0 0.5

1 Original function

Time (s)

Amplitude

0 0.5 1

0 0.2

Time (s)

Amplitude

0 0.5 1

0 0.1 0.2

Time (s)

Amplitude

0 0.5 1

0 0.1

Filtered with n=9, f0=20 Hz 0.2

Time (s)

Amplitude

(32)

Butterworth Low-pass (20 Hz) on data

40 45 50 55 60 65

-5 0 5

Time (s)

Amplitude

40 45 50 55 60 65

-1 -0.5 0 0.5 1

Time (s)

Amplitude

40 45 50 55 60 65

-2 -1 0 1 2

Time (s)

Amplitude

40 45 50 55 60 65

-2 0 2

Time (s)

Amplitude

(33)

33

Other windowing functions

 So far we only used the Butterworth filtering window

 In general if we want to extract time windows from (permanent) recordings we have other options in the time domain.

 The key issues are

Do you want to preserve the main maxima at the expense of side maxima?

Do you want to have as little side lobes as posible?

(34)

Example

(35)

35

Possible windows

Plain box car (arrow stands for Fourier transform):

fM M fM

f M W

t

M t t

w_R _R



 2

2 2 sin

) , (

0 , ) 1

(  











 

Bartlett

2 2

) sin (

, 0

, ) 1

( 



 



 

 















 

fM M fM

f W

M t

M M t

t t

w_R _B



(36)

Possible windows

Hanning

The spectral representations of the boxcar, Bartlett (and Parzen) functions are:

)2

2 ( 1

1 2

2 ) sin

( ,

0

), cos

1 2 ( 1 )

( fM fM

M fM f

W M

t

M M t

t t

w_H _H



 

 

 















 

2 /

; 4 , 2 , 2 1

2 ) sin

( n M T

fM f fM

W

n



 



 



 



(37)

37

Examples

(38)

Examples

(39)

39

The Gabor transform: t-f misfits

^

 



 u( )g(τ t)e dτ 2π

G(u) 1 :

t) ,

û(  ^i

^

 



 u ( )g(τ t)e dτ 2π

) 1 G(u :

t) , (

û₀  ₀ ₀  ⁱ^

[ t- representation of synthetics, u(t) ]

[ t- representation of data, u₀(t) ] phase information:

• can be measured reliably

• ± linearly related to Earth structure

• physically interpretable

amplitude information:

• hard to measure (earthquake magnitude often unknown)

• non-linearly related to structure

(40)

The Gabor time window

The Gaussian time windows is given by

(41)

41

Example

(42)

Multitaper

Goal: „obtaining a spectrum with little or no bias and small uncertainties“.

problem comes down to finding the right tapering to reduce the bias (i.e, spectral leakage).

In principle we seek:

This section follows Prieto eet al., GJI, 2007. Ideas go back to a paper by Thomson (1982).

(43)

43

Multi-taper Principle

• Data sequence x is multiplied by a set of orthgonal sequences (tapers)

• We get several single periodograms (spectra) that are then averaged

• The averaging is not even, various weights apply

• Tapers are constructed to optimize resistance to spectral leakage

• Weighting designed to generate smooth estimate with less variance than with single tapers

(44)

Spectrum estimates

We start with

with

To maintain total power.

(45)

45

Condition for optimal tapers

N is the number of points, W is the resolution bandwith (frequency increment)

One seeks to maximize  the fraction of energy in the interval (–W,W). From this equation one finds a‘s by an eigenvalue problem -> Slepian function

(46)

Slepian functions

The tapers (Slepian functions) in time and frequency domains

(47)

47

Final assembly

Final averaging of spectra Slepian sequences (tapers)

(48)

Example

0 5 10 15 20 25 30 35 40 45 50

-6 -4 -2 0 2 4 6

8 Signal Corrupted with Zero-Mean Random Noise

time (milliseconds)

(49)

49

Classical Periodogram

0 50 100 150 200 250 300 350 400 450 500

0 0.2 0.4 0.6 0.8 1 1.2

1.4 Single-Sided Amplitude Spectrum of y(t)

Frequency (Hz)

|Y(f)|

(50)

… and its power …

0 50 100 150 200 250 300 350 400 450 500

10^-6 10^-5 10^-4 10^-3 10^-2 10^-1 10⁰ 10¹

Single-Sided Power Spectrum of y(t)

Frequency (Hz)

|Y(f)|

(51)

51

… multitaper spectrum …

0 50 100 150 200 250 300 350 400 450 500

-160 -140 -120 -100 -80 -60 -40 -20 0

Frequency (Hz)

Power/frequency (dB/Hz)

Power Spectral Density

(52)

Wavelets – the principle

Motivation:

Time-frequency analysis

Multi-scale approach

„when do we hear what frequency?“

(53)

53

Continuous vs. local basis functions

(54)

(55)

55

Some maths

A wavelet can be defined as



 



 



 ^

a b a t

b t

a, ( ) 1/2







 



 



 

 dt

a b t x

a f b

a f

W 1 ( )

) , )(

(

dadb a

t f

C t

f ( ) ^ , ^a^,^b ^a^,^b( ) ^²







   

  

With the transform pair:

(56)

(57)

57

Resulting wavelet representation

(58)

Shifting and scaling

(59)

59

(60)

Application to seismograms

http://users.math.uni-potsdam.de/~hols/DFG1114/projectseis.html

(61)

61

Graphical comparison

(62)

Summary

 Filtering is not necessarily straight forward, even the fundamental operations (LP, HP, BP, etc) require some thinking before application to data.

 The form of the filter decides upon the changes to the waveforms of the time series you are filtering

 For seismological applications filtering might drastically influence observables such as travel times or amplitudes

 „Windowing“ the signals in the right way is fundamental to obtain the desired filtered sequence

(63)

63