Transcript Spectral analysis: Foundations

Spectral analysis: Foundations







Orthogonal functions
Fourier Series
Discrete Fourier Series
Fourier Transform: properties
Chebyshev polynomials
Convolution
DFT and FFT
Scope: Understanding where the Fourier Transform comes
from. Moving from the continuous to the discrete world.
(Almost) everything we need to understand for filtering.
Spectral analysis: foundations
Computational Geophysics and Data Analysis
1
Fourier Series: one way to derive them
The Problem
we are trying to approximate a function f(x) by another function gn(x)
which consists of a sum over N orthogonal functions F(x) weighted by
some coefficients an.
N
f ( x)  g N ( x)   ai F i ( x)
i 0
Spectral analysis: foundations
Computational Geophysics and Data Analysis
2
The Problem
... and we are looking for optimal functions in a least squares (l2) sense ...


2
f ( x)  g N ( x)   f ( x)  g N ( x) dx
2
a

b


1/ 2
 Min!
... a good choice for the basis functions F(x) are orthogonal functions.
What are orthogonal functions? Two functions f and g are said to be
orthogonal in the interval [a,b] if
b
 f ( x) g ( x)dx  0
a
How is this related to the more conceivable concept of orthogonal
vectors? Let us look at the original definition of integrals:
Spectral analysis: foundations
Computational Geophysics and Data Analysis
3
Orthogonal Functions
N

  f i ( x) gi ( x)x 
a f ( x) g ( x)dx  Nlim

 i 1

b
... where x0=a and xN=b, and xi-xi-1=x ...
If we interpret f(xi) and g(xi) as the ith components of an N component
vector, then this sum corresponds directly to a scalar product of vectors.
The vanishing of the scalar product is the condition for orthogonality of
vectors (or functions).
fi
Spectral analysis: foundations
gi
f i  gi   f i gi  0
i
Computational Geophysics and Data Analysis
4
Periodic functions
Let us assume we have a piecewise continuous function of the form
f ( x  2 )  f ( x)
40
f ( x  2 )  f ( x)  x 2
30
20
10
0
-15
-10
-5
0
5
10
15
20
... we want to approximate this function with a linear combination of 2
periodic functions:
1, cos(x), sin(x), cos(2 x), sin(2 x),...,cos(nx), sin(nx)
N
1
 f ( x)  g N ( x)  a0   ak cos(kx)  bk sin(kx)
2
k 1
Spectral analysis: foundations
Computational Geophysics and Data Analysis
5
Orthogonality
... are these functions orthogonal ?
jk
 0

cos(
jx
)
cos(
kx
)
dx

2 j  k  0


jk 0

 0 j  k , j, k  0


sin(
jx
)
sin(
kx
)
dx




jk 0



 cos( jx) sin(kx)dx  0
j  0, k  0

... YES, and these relations are valid for any interval of length 2.
Now we know that this is an orthogonal basis, but how can we obtain the
coefficients for the basis functions?
from minimising f(x)-g(x)
Spectral analysis: foundations
Computational Geophysics and Data Analysis
6
Fourier coefficients
optimal functions g(x) are given if
g n ( x)  f ( x)
2
 Min!
or

ak
g
n
( x)  f ( x )
2
 0
... with the definition of g(x) we get ...

g n ( x)  f ( x)
ak
2
2
2

N

  1


   a0   ak cos(kx)  bk sin(kx) f ( x) dx
ak   2
k 1




leading to
N
1
g N ( x)  a0   ak cos(kx)  bk sin(kx)
2
k 1
ak 
bk 
Spectral analysis: foundations
1

1

wit h

 f ( x) cos(kx)dx,
k  0,1,..., N


 f ( x) sin(kx)dx,
k  1,2,..., N

Computational Geophysics and Data Analysis
7
Fourier approximation of |x|
... Example ...
  x  
f ( x)  x ,
leads to the Fourier Serie
g ( x) 
1
4  cos(x) cos(3x) cos(5 x)

 



...

2
  12
32
52

.. and for n<4 g(x) looks like
4
3
2
1
0
-20
-15
Spectral analysis: foundations
-10
-5
0
5
10
Computational Geophysics and Data Analysis
15
20
8
Fourier approximation of x2
... another Example ...
0  x  2
f ( x)  x 2 ,
leads to the Fourier Serie
N
4 2
4
4

g N ( x) 
   2 cos(kx) 
sin(kx)
3
k

k 1  k
.. and for N<11, g(x) looks like
40
30
20
10
0
-10
-10
Spectral analysis: foundations
-5
0
5
Computational Geophysics and Data Analysis
10
15
9
Fourier - discrete functions
... what happens if we know our function f(x) only at the points
2
xi 
i
N
it turns out that in this particular case the coefficients are given by
2
ak 
N
*
2
bk 
N
*
N
 f (x
j 1
j
) cos(kx j ),
k  0,1,2,...
j
) sin(kx j ),
k  1,2,3,...
N
 f (x
j 1
.. the so-defined Fourier polynomial is the unique interpolating function to
the function f(xj) with N=2m


1 * m1 *
1
g ( x)  a 0   a k cos(kx)  b*k sin(kx)  am* cos(kx)
2
2
k 1
*
m
Spectral analysis: foundations
Computational Geophysics and Data Analysis
10
Fourier - collocation points
... with the important property that ...
gm* ( xi )  f ( xi )
... in our previous examples ...
3.5
3
2.5
2
1.5
1
0.5
0
-10
-5
0
5
10
f(x)=|x| => f(x) - blue ; g(x) - red; xi - ‘+’
Spectral analysis: foundations
Computational Geophysics and Data Analysis
11
Fourier series - convergence
f(x)=x2 => f(x) - blue ; g(x) - red; xi - ‘+’
Spectral analysis: foundations
Computational Geophysics and Data Analysis
12
Fourier series - convergence
f(x)=x2 => f(x) - blue ; g(x) - red; xi - ‘+’
Spectral analysis: foundations
Computational Geophysics and Data Analysis
13
Gibb’s phenomenon
f(x)=x2 => f(x) - blue ; g(x) - red; xi - ‘+’
N = 16
N = 64
N = 32
6
6
6
4
4
4
2
2
2
0
0
0
-2
-2
-2
-4
-4
-4
-6
0
0.5
1
N = 128
1.5
-6
0
6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6
0
0.5
1
Spectral analysis: foundations
1.5
-6
0.5
1
N = 256
1.5
-6
0
0.5
1
1.5
The overshoot for equispaced Fourier
interpolations is 14% of
the step height.
0
0.5
1
1.5
Computational Geophysics and Data Analysis
14
Chebyshev polynomials
We have seen that Fourier series are excellent for interpolating
(and differentiating) periodic functions defined on a regularly
spaced grid. In many circumstances physical phenomena which
are not periodic (in space) and occur in a limited area. This quest
leads to the use of Chebyshev polynomials.
We depart by observing that cos(n) can be expressed by a
polynomial in cos():
cos(2 )  2 cos2   1
cos(3 )  4 cos3   3 cos
cos(4 )  8 cos4   8 cos2   1
... which leads us to the definition:
Spectral analysis: foundations
Computational Geophysics and Data Analysis
15
Chebyshev polynomials - definition
cos(n )  Tn (cos( ))  Tn ( x),
x  cos( ),
x [1,1],
n N
... for the Chebyshev polynomials Tn(x). Note that because of
x=cos() they are defined in the interval [-1,1] (which - however can be extended to ). The first polynomials are
T0 ( x)  1
T1 ( x)  x
T2 ( x)  2 x 2  1
T3 ( x)  4 x 3  3 x
T4 ( x)  8 x 4  8 x 2  1
Tn ( x)  1
Spectral analysis: foundations
where
for x  [1,1]
and n  N 0
Computational Geophysics and Data Analysis
16
Chebyshev polynomials - Graphical
The first ten polynomials look like [0, -1]
1
T_n(x)
0.5
0
-0.5
-1
0
0.2
0.6
0.4
0.8
1
x
The n-th polynomial has extrema with values 1 or -1 at
x
Spectral analysis: foundations
( ext)
k
k
 cos(
),
n
k  0,1,2,3,...,n
Computational Geophysics and Data Analysis
17
Chebyshev collocation points
These extrema are not equidistant (like the Fourier extrema)
k
x(k)
x
Spectral analysis: foundations
( ext)
k
k
 cos(
),
n
k  0,1,2,3,...,n
Computational Geophysics and Data Analysis
18
Chebyshev polynomials - orthogonality
... are the Chebyshev polynomials orthogonal?
Chebyshev polynomials are an orthogonal set of functions in the
interval [-1,1] with respect to the weight function 1/ 1 x2
such that
0

dx

T
(
x
)
T
(
x
)

 / 2
1 k j
2
1 x
 

1
k j


for k  j  0,
for k  j  0 

for
k, j  N0
... this can be easily verified noting that
x  cos , dx   sin d
Tk ( x)  cos(k ), T j ( x)  cos( j )
Spectral analysis: foundations
Computational Geophysics and Data Analysis
19
Chebyshev polynomials - interpolation
... we are now faced with the same problem as with the Fourier
series. We want to approximate a function f(x), this time not a
periodical function but a function which is defined between [-1,1].
We are looking for gn(x)
n
1
f ( x)  g n ( x)  c0T0 ( x)   ck Tk ( x)
2
k 1
... and we are faced with the problem, how we can determine the
coefficients ck. Again we obtain this by finding the extremum
(minimum)
1
 
dx 
2
  g n ( x)  f ( x)
0
2
ck  1
1 x 
Spectral analysis: foundations
Computational Geophysics and Data Analysis
20
Chebyshev polynomials - interpolation
... to obtain ...
ck 
2

1

f ( x)Tk ( x)
1
dx
1 x
2
,
k  0,1,2,...,n
... surprisingly these coefficients can be calculated with FFT
techniques, noting that
ck 
2


f (cos ) cos kd ,
k  0,1,2,...,n
0
... and the fact that f(cos) is a 2-periodic function ...
ck 
1



f (cos ) cos kd ,
k  0,1,2,...,n

... which means that the coefficients ck are the Fourier coefficients
ak of the periodic function F()=f(cos )!
Spectral analysis: foundations
Computational Geophysics and Data Analysis
21
Chebyshev - discrete functions
... what happens if we know our function f(x) only at the points
xi  cos

N
i
in this particular case the coefficients are given by
2
ck 
N
*
N
 f (cos
j 1
j
) cos(k j ),
k  0,1,2,...N / 2
... leading to the polynomial ...
m
1 *
g ( x)  c0T 0   ck*Tk ( x)
2
k 1
*
m
... with the property
gm* ( x)  f ( x)
Spectral analysis: foundations
at
x j  cos(j/N)
j  0,1,2,...,N
Computational Geophysics and Data Analysis
22
Chebyshev - collocation points - |x|
f(x)=|x| => f(x) - blue ; gn(x) - red; xi - ‘+’
N=8
1
0.8
8 points
0.6
0.4
0.2
0
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
N = 16
1
0.8
0.6
16 points
0.4
0.2
0
-1
Spectral analysis: foundations
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Computational Geophysics and Data Analysis
0.8
1
23
Chebyshev - collocation points - |x|
f(x)=|x| => f(x) - blue ; gn(x) - red; xi - ‘+’
N = 32
1
0.8
0.6
32 points
0.4
0.2
0
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.2
0.4
0.6
0.8
1
N = 128
1
0.8
128 points
0.6
0.4
0.2
0
-1
Spectral analysis: foundations
-0.8
-0.6
-0.4
-0.2
0
0.6
Computational Geophysics and Data Analysis
0.8
1
24
Chebyshev - collocation points - x2
f(x)=x2 => f(x) - blue ; gn(x) - red; xi - ‘+’
N=8
1.2
1
0.8
8 points
0.6
0.4
0.2
0
-1
-0.8
-0.6
-0.4
-0.2
0
N = 64
0.2
0.4
0.6
0.8
1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1
The interpolating
function gn(x) was
shifted by a small
amount to be
visible at all!
0.8
0.6
64 points
0.4
0.2
0
-1
Spectral analysis: foundations
Computational Geophysics and Data Analysis
25
Chebyshev vs. Fourier - numerical
Chebyshev
Fourier
N = 16
N = 16
1
35
0.8
30
25
0.6
20
15
0.4
10
5
0.2
0
0
0
0.2
0.4
0.6
0.8
1
-5
0
2
4
6
f(x)=x2 => f(x) - blue ; gN(x) - red; xi - ‘+’
This graph speaks for itself ! Gibb’s phenomenon with Chebyshev?
Spectral analysis: foundations
Computational Geophysics and Data Analysis
26
Chebyshev vs. Fourier - Gibb’s
Chebyshev
Fourier
N = 16
N = 16
1.5
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5
-1
-0.5
0
0.5
1
-1.5
0
2
4
6
f(x)=sign(x-) => f(x) - blue ; gN(x) - red; xi - ‘+’
Gibb’s phenomenon with Chebyshev? YES!
Spectral analysis: foundations
Computational Geophysics and Data Analysis
27
Chebyshev vs. Fourier - Gibb’s
Chebyshev
Fourier
N = 64
N = 64
1.5
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5
-1
-0.5
0
0.5
1
-1.5
0
2
4
6
8
f(x)=sign(x-) => f(x) - blue ; gN(x) - red; xi - ‘+’
Spectral analysis: foundations
Computational Geophysics and Data Analysis
28
Fourier vs. Chebyshev
Chebyshev
Fourier
2
xi 
i
N
periodic functions
domain
1 *
a0
2
m 1


  a *k cos(kx)  b*k sin(kx)
interpolating
function

N
i
limited area [-1,1]
Tn ( x)  cos(n ),
basis functions
cos(nx), sin(nx)
g m* ( x ) 
xi  cos
collocation points
x  cos
m
1 *
g ( x)  c0T 0   ck*Tk ( x)
2
k 1
*
m
k 1

1 *
am cos(kx)
2
Spectral analysis: foundations
Computational Geophysics and Data Analysis
29
Fourier vs. Chebyshev (cont’d)
Chebyshev
Fourier
2
ak 
N
*
2
bk 
N
*
N
 f (x
j 1
j
) cos(kx j )
coefficients
N
 f (x
j 1
j
2
ck 
N
*
N
 f (cos
j 1
j
) cos(k j )
) sin(kx j )
• Gibb’s phenomenon for
discontinuous functions
• Efficient calculation via FFT
• infinite domain through
periodicity
• limited area calculations
some properties
• grid densification at boundaries
• coefficients via FFT
• excellent convergence at
boundaries
• Gibb’s phenomenon
Spectral analysis: foundations
Computational Geophysics and Data Analysis
30
The Fourier Transform Pair
1
F ( ) 
2

 f (t )e
i t
dt
Forward transform


f (t ) 
 i t
F
(

)
e
d

Inverse transform

Note the conventions concerning the sign of the exponents and the factor.
Spectral analysis: foundations
Computational Geophysics and Data Analysis
31
The Fourier Transform Pair
F ( )  R( )  iI ( )  A( )e
iF (  )
A( )  F ( )  R 2 ( )  I 2 ( )
I ( )
F ( )  arg F ( )  arctan
R( )
A( )
F( )
Amplitude spectrum
Phase spectrum
In most application it is the amplitude (or the power) spectrum that is of interest.
Spectral analysis: foundations
Computational Geophysics and Data Analysis
32
The Fourier Transform: when does it work?
Conditions that the integral transforms work:
 f(t) has a finite number of jumps and the limits exist from
both sides

 f(t) is integrable, i.e.

f (t ) dt  G  

Properties of the Fourier transform for special functions:
Function f(t)
Fouriertransform F()
even
even
odd
odd
real
hermitian
imaginary
antihermitian
hermitian
real
Spectral analysis: foundations
Computational Geophysics and Data Analysis
33
… graphically …
Spectral analysis: foundations
Computational Geophysics and Data Analysis
34
Some properties of the Fourier Transform
Defining as the FT:
af1 (t )  bf2 (t )  aF1 ()  bF2 ()
 Linearity
f (t )  2F ( )
 Symmetry
f (t  t )  eit F ()
 Time shifting
 Time differentiation
Spectral analysis: foundations
f (t )  F ( )
 n f (t )
n

(

i

)
F ( )
n
t
Computational Geophysics and Data Analysis
35
Differentiation theorem
 Time differentiation
Spectral analysis: foundations
 n f (t )
n

(

i

)
F ( )
n
t
Computational Geophysics and Data Analysis
36
Convolution
The convolution operation is at the heart of linear systems.
Definition:
f (t )  g (t ) 
Properties:




 f (t ' ) g (t  t ' )dt'   f (t  t ' ) g (t ' )dt'
f (t )  g (t )  g (t )  f (t )
f (t )  (t )  f (t )
f (t )  H (t )   f (t )dt
H(t) is the Heaviside function:
Spectral analysis: foundations
Computational Geophysics and Data Analysis
37
The convolution theorem
A convolution in the time domain corresponds to a
multiplication in the frequency domain.
… and vice versa …
a convolution in the frequency domain corresponds to a
multiplication in the time domain
f (t )  g (t )  F ( )G( )
f (t ) g (t )  F ( )  G( )
The first relation is of tremendous practical implication!
Spectral analysis: foundations
Computational Geophysics and Data Analysis
38
The convolution theorem
From Bracewell (Fourier transforms)
Spectral analysis: foundations
Computational Geophysics and Data Analysis
39
Discrete Convolution
Convolution is the mathematical description of the change of
waveform shape after passage through a filter (system).
There is a special mathematical symbol for convolution (*):
y(t )  g (t )  f (t )
Here the impulse response function g is convolved with the
input signal f. g is also named the „Green‘s function“
yk   g i f k i
gi
i  0,1,2,....,m
k  0,1,2,, m  n
fj
j  0,1,2,....,n
m
i 0
Spectral analysis: foundations
Computational Geophysics and Data Analysis
40
Convolution Example(Matlab)
Impulse response
>> x
x=
0
0
1
0
>> y
System input
y=
1
2
1
>> conv(x,y)
ans =
0
Spectral analysis: foundations
0
1
2
1
0
Computational Geophysics and Data Analysis
System output
41
Convolution Example (pictorial)
„Faltung“
x
0
0
1
0
2
1
1
0
0
1
2
1
0
1
0
0
1
2
1
0
2
0
2
Spectral analysis: foundations
0
2
0
1
0
x*y
1
0
1
1
y
1
1
1
1
0
0
0
0
1
0
0
1
2
1
0
1
Computational Geophysics and Data Analysis
42
The digital world
Spectral analysis: foundations
Computational Geophysics and Data Analysis
43
The digital world

g s (t )  g (t )   (t  jdt)
j  
gs is the digitized version of g and the sum is called the comb function.
Defining the Nyquist frequency fNy as
f Ny 
1
2dt
after a few operations the spectrum can be written as



1

Gs ( f )  G( f )   G( f  2nfNy )  G( f  2nfNy ) 
dt 
n 1

… with very important consequences …
Spectral analysis: foundations
Computational Geophysics and Data Analysis
44
The sampling theorem
The implications are that for the calculation of the
spectrum at frequency f there are also contributions
of frequencies f±2nfNy, n=1,2,3,…
f Ny
That means dt has to be chosen such that fN is the
largest frequency contained in the signal.
Spectral analysis: foundations
Computational Geophysics and Data Analysis
1

2dt
45
The Fast Fourier Transform FFT
... spectral analysis became interesting for computing 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
N 1
Fk   f j e  2ikj / N , k  0,1,..., N  1
j 0
1
fk 
N
Spectral analysis: foundations
N 1
2ikj / N
F
e
, k  0,1,..., N  1
 j
j 0
Computational Geophysics and Data Analysis
46
The Fast Fourier Transform FFT
... this can be written as matrix-vector products ...
for example the inverse transform yields ...
1
1
1


1  2






N 1
1 
  F0   f 0 
  N 1   F1   f1 
  2 N  2   F2   f 2 



      
      
 

( N 1) 2  
 
  FN 1   f N 1 

1
1
2
4
3
6


1
.. where ...
  e 2i / N
Spectral analysis: foundations
Computational Geophysics and Data Analysis
47
FFT
... 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).
Spectral analysis: foundations
Computational Geophysics and Data Analysis
48
FFT
Number of multiplications
Problem
1D (nx=512)
1D (nx=2096)
1D (nx=8384)
full matrix
2.6x105
FFT
9.2x103
Ratio full/FFT
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.
Spectral analysis: foundations
Computational Geophysics and Data Analysis
49
Summary
 The Fourier Transform can be derived from the problem of approximating
an arbitrary function.
 A regular set of points allows exact interpolation (or derivation) of arbitrary
functions
 There are other basis functions (e.g., Chebyshev polynomials) with similar
properties
 The discretization of signals has tremendous impact on the estimation of
spectra: aliasing effect
 The FFT is at the heart of spectral analysis
Spectral analysis: foundations
Computational Geophysics and Data Analysis
50