[PPT] - We have developed numerical approaches to the following: Fourier PowerPoint Presentation

SLIDE 1

SIGNAL SMOOTHING AND ROOT FINDING

Lesson 9

SLIDE 2

We have developed numerical approaches to the following:
Fourier series on the periodic interval
Taylor series on the circle
Chebyshev series on the interval
Consider two applications:
Smoothing a noisy signal
Finding roots of a function

SLIDE 3

WHITE NOISE

SLIDE 4

Return to the case of sampling a function at evenly spaced points on the periodic interval
3
2
1

1 2 3

f(θ1) f(θm)

…

Suppose each time the function is evaluated it returns an identically distributed random number, say, in the interval

(–1,1)

I.e., the evaluation vector is a random vector:
3
2
1

1 2 3

… r1 rm

SLIDE 5

Instances of white noise

Function values
3
2
1

1 2 3

1.0
0.5

0.5 1.0

m = 5

SLIDE 6

Instances of white noise

Function values
3
2
1

1 2 3

1.0
0.5

0.5 1.0

m = 10

SLIDE 7

Instances of white noise

Function values
3
2
1

1 2 3

1.0
0.5

0.5 1.0

m = 20

SLIDE 8

Instances of white noise

Function values
3
2
1

1 2 3

1.0
0.5

0.5 1.0

m = 40

SLIDE 9

3
2
1

1 2 3

1.0
0.5

0.5 1.0

m = 100 Instances of white noise

Function values

SLIDE 10

White noise DFT coefficients

1

1 2 0.05 0.10 0.20 0.50 1.00

m = 5

Computed Fourier coefficients

SLIDE 11

White noise DFT coefficients

4
2

2 4 0.10 1.00 0.50 0.20 0.30 0.15 0.70

m = 10

Computed Fourier coefficients

SLIDE 12

White noise DFT coefficients

5

5 0.02 0.05 0.10 0.20 0.50 1.00

m = 20

Computed Fourier coefficients

SLIDE 13

White noise DFT coefficients

10

10 0.01 0.02 0.05 0.10 0.20 0.50 1.00

m = 40

Computed Fourier coefficients

SLIDE 14

White noise DFT coefficients

40
20

20 40 0.01 0.02 0.05 0.10 0.20 0.50 1.00

m = 100

Computed Fourier coefficients

SLIDE 15

White noise DFT coefficients

400
200

200 400 0.005 0.010 0.050 0.100 0.500 1.000

m = 1000

Computed Fourier coefficients

SLIDE 16

4000
2000

2000 4000 0.001 0.005 0.010 0.050 0.100 0.500 1.000

m = 10000

White noise DFT coefficients Computed Fourier coefficients

SLIDE 17

SIGNAL DENOISING

SLIDE 18

The Fourier coefficients are the same magnitude!
But unlike the values, they tend to zero as the number of samples increases
Thus, by computing the Fourier coefficients, we can remove white noise from a signal

SLIDE 19

3
2
1

1 2 3

1

1 2 3

m = 50

Original Signal Signal with noise

3
2
1

1 2 3

1

1 2 3

m = 50

SLIDE 20

Noisy signal Original signal DFT coefficients

20
10

10 20 10-14 10-11 10-8 10-5 0.01

m = 50

SLIDE 21

Noisy signal Original signal DFT coefficients

20
10

10 20 10-15 10-12 10-9 10-6 0.001 1

m = 100

SLIDE 22

Noisy signal Original signal DFT coefficients

20
10

10 20 10-15 10-12 10-9 10-6 0.001 1

m = 1000

SLIDE 23

Noisy signal Original signal DFT coefficients

20
10

10 20 10-16 10-12 10-8 10-4 1

m = 10000

SLIDE 24

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients

3
2
1

1 2 3 0.5 1.0 1.5 2.0 2.5

m = 100

SLIDE 25

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients

3
2
1

1 2 3 0.5 1.0 1.5 2.0 2.5

m = 1000

SLIDE 26

3
2
1

1 2 3 0.5 1.0 1.5 2.0 2.5

m = 10000

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients

SLIDE 27

3
2
1

1 2 3 0.5 1.0 1.5 2.0 2.5

m = 100000

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients

SLIDE 28

TAYLOR SERIES ROOT FINDING

SLIDE 29

Consider a Taylor series

f(z) =

∞

k=0

ˆ fkzk

We want to find the all roots of this function in the unit circle, i.e., all z such that

f(z) = 0

The first step is to truncate and scale the Taylor series, so that we are dealing with monic polynomials:

f(z) ˆ fm ≈ p(z) =

m−1

k=0

ˆ fk ˆ fm zk + zm =

m−1

k=0

ˆ pkzk + zm

SLIDE 30

We will rephrase the root finding problem as a eigenvalue problem
Consider the companion matrix
We will show that the eigenvalues of the companion matrix are the roots of p

     1 ... 1 −ˆ p0 −ˆ p1 · · · −ˆ pm−1     

SLIDE 31

Given z, we have

     1 ... 1 −ˆ p0 −ˆ p1 · · · −ˆ pm−1           1 . . . zm−2 zm−1      =      z . . . zm−1 −ˆ p0 − ˆ p1z − · · · − ˆ pm−1zm−1     

If z is a root of p, then this simplifies:

     z . . . zm−1 −ˆ p0 − ˆ p1z − · · · − ˆ pm−1zm−1      =      z . . . zm−1 zm      = z      1 . . . zm−2 zm−1     

Eigenvalue! Eigenvector!

SLIDE 32

Thus if z1,…,zm denote the m roots of p, we have

     1 ... 1 −ˆ p0 −ˆ p1 · · · −ˆ pm−1      V = V    z1 ... zm   

for

V =    1 · · · 1 . . . ... . . . zm−1

1

· · · zm−1

m

  

SLIDE 33

Why is this useful? Because there are well-developed, reliable algorithms for computing eigenvalues of matrices
For example, eigs in MATLAB
Thus we have a very simple algorithm for calculating roots of analytic functions:
Step 1: approximate the Taylor series of f on the unit circle using the discrete Taylor transform
Step 2: construct the companion matrix A
Step 3: call eigs(A)

SLIDE 34

Example: sin 4x =

4
2

2 4

4
2

2 4

m = 4

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 35

4
2

2 4

4
2

2 4

m = 8

Example: sin 4x =

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 36

Example: sin 4x =

4
2

2 4

4
2

2 4

m = 16

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 37

Example: sin 4x =

4
2

2 4

4
2

2 4

m = 32

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 38

Example: sin 4x =

4
2

2 4

4
2

2 4

m = 40

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 39

Example: sin 4x =

4
2

2 4

4
2

2 4

m = 50

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 40

4
2

2 4

4
2

2 4

m = 100

Example: sin 4x =

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 41

LAURENT SERIES ROOT FINDING

SLIDE 42

In the last example, we saw that the roots were accurate in the domain that the Taylor

series was accurate

Now, we consider root finding of Laurent series
Where can we expect the roots to be accurate?
On the unit circle itself
We simply convert an approximate Laurent series to a polynomial and use the previous

approach:

z−α ˆ fβ f(z) ≈ z−α ˆ fβ

ˆ

fαzα + · · · + ˆ fβzβ = ˆ fα ˆ fβ + · · · + ˆ fβ−1 ˆ fβ zm−2 + zm−1 = p(z)

SLIDE 43

We will adapt this to computing roots of functions on the periodic interval
Recall for g defined on the unit circle, f(θ) = g(θ) was defined on the periodic

interval

Let f ∈ ∞[T]. We want to convert it to a g defined on the unit circle
This is straightforward:

g(z) = f(− z) so that g(θ) = f(θ)

The periodicity of f implies that the choice of branch cut in z doesn't matter
We can thus find the roots z1, . . . , zr of g that lie on the unit circle
Then − z1, . . . , − zr are roots of f on the periodic interval

SLIDE 44

3
2
1

1 2 3

1.0
0.5

0.5 1.0

f(θ) = cos(50 cos2 θ) cos θ

1.5
1.0
0.5

0.5 1.0 1.5

1.5
1.0
0.5

0.5 1.0 1.5

f(−i log z)

Numerical roots of

SLIDE 45

3
2
1

1 2 3

1.0
0.5

0.5 1.0

f(θ) = cos(50 cos2 θ) cos θ f(−i log z)

Pruned numerical roots of

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

SLIDE 46

f(θ) = cos(50 cos2 θ) cos θ

f(−i log z)

Pruned numerical roots of

1.0
0.5

0.5 1.0

1.0
0.5

0.5 1.0

3
2
1

1 2 3

1.0
0.5

0.5 1.0

Roots mapped to periodic interval

SLIDE 47

CHEBYSHEV SERIES ROOT FINDING

SLIDE 48

48

We want to find the roots of

0 = f(x) =

∞

k=0

ˇ fkTk(x)

The easy way:

Recall the Joukowsky map J(z) = 1

2

z + 1

z

and let

g(z) = f(J(z)) = 1 2

−1

k=−∞

ˇ f−kzk + ˇ f0 + 1 2

∞

k=1

ˇ fkzk Truncate g and find the roots z1, . . . , zr that lie on the lower half of the unit circle The roots of f are then J(z1), . . . , J(zr)

The problem: this requires a matrix of size 2n × 2n instead of n × n

SLIDE 49

THREE-TERM RECURRENCE RELATIONSHIP

SLIDE 50

Suppose a family of polynomials p0(x), p1(x), . . . (where pk is of degree k) are
rthogonal with respect to an inner product

f, g = 1

−1

¯ f(x)g(x)w(x) x where w(x) is a weighting function In our case, pk(x) = Tk(x) and w(x) =

1 √ 1−x2

We will show that there exists constants

, and such that

First note that

is a polynomial of degree and is degree , so we have where is of degree We want to show that

SLIDE 51

Suppose a family of polynomials p0(x), p1(x), . . . (where pk is of degree k) are
rthogonal with respect to an inner product

f, g = 1

−1

¯ f(x)g(x)w(x) x where w(x) is a weighting function In our case, pk(x) = Tk(x) and w(x) =

1 √ 1−x2

We will show that there exists constants Ak, Bk and Ck such that

pk+1(x) = (Akx + Bk)pk Ckpk−1

First note that

is a polynomial of degree and is degree , so we have where is of degree We want to show that

SLIDE 52

Suppose a family of polynomials p0(x), p1(x), . . . (where pk is of degree k) are
rthogonal with respect to an inner product

f, g = 1

−1

¯ f(x)g(x)w(x) x where w(x) is a weighting function In our case, pk(x) = Tk(x) and w(x) =

1 √ 1−x2

We will show that there exists constants Ak, Bk and Ck such that

pk+1(x) = (Akx + Bk)pk Ckpk−1

First note that pk+1 is a polynomial of degree k+1 and pk is degree k, so we have

pk+1(x) = Akxpk(x) + Bkpk(x) Ckpk−1(x) + qk−2(x) where qk−2 is of degree qk−2 We want to show that qk−2 = 0

SLIDE 53

We can write

qk−2(x) =

k−2

j=0

cjpj(x) where cj = pj, qk−2 pj2 .

Rearrange terms we also have
For

we have

SLIDE 54

We can write

qk−2(x) =

k−2

j=0

cjpj(x) where cj = pj, qk−2 pj2 .

Rearrange terms we also have

qk−2(x) = pk+1(x) + Ckpk−1(x) Akxpk(x) Bkpk(x)

For

we have

SLIDE 55

55

We can write

qk−2(x) =

k−2

j=0

cjpj(x) where cj = pj, qk−2 pj2 .

Rearrange terms we also have

qk−2(x) = pk+1(x) + Ckpk−1(x) Akxpk(x) Bkpk(x)

For j = 0, . . . , k 2 we have

pj, qk−2 = pj, pk+1 + Ck pj, pk−1 Ak pj, xpk Bk pj, pk

Zero due to orthogonality!

SLIDE 56

For Chebyshev polynomials, the three-term recurrence is very simple:

T1(x) = xT0(x) Tk+1(x) = 2xTk(x) − Tk−1(x)

We prove this by letting

and recalling :

SLIDE 57

For Chebyshev polynomials, the three-term recurrence is very simple:

T1(x) = xT0(x) Tk+1(x) = 2xTk(x) − Tk−1(x)

We prove this by letting x = J(z) and recalling Tk(J(z)) = 1

2

zk + z−k

: 2J(z)Tk(J(z)) − Tk−1(J(z)) = 1 2

(z + z−1)(zk + z−k) − (zk−1 + z1−k)

SLIDE 58

For Chebyshev polynomials, the three-term recurrence is very simple:

T1(x) = xT0(x) Tk+1(x) = 2xTk(x) − Tk−1(x)

We prove this by letting x = J(z) and recalling Tk(J(z)) = 1

2

zk + z−k

: 2J(z)Tk(J(z)) − Tk−1(J(z)) = 1 2

(z + z−1)(zk + z−k) − (zk−1 + z1−k)
= 1

2

zk+1 + z−k−1 + zk−1 + z1−k − (zk−1 + z1−k)
1

SLIDE 59

For Chebyshev polynomials, the three-term recurrence is very simple:

T1(x) = xT0(x) Tk+1(x) = 2xTk(x) − Tk−1(x)

We prove this by letting x = J(z) and recalling Tk(J(z)) = 1

2

zk + z−k

: 2J(z)Tk(J(z)) − Tk−1(J(z)) = 1 2

(z + z−1)(zk + z−k) − (zk−1 + z1−k)
= 1

2

zk+1 + z−k−1 + zk−1 + z1−k − (zk−1 + z1−k)
= 1

2

zk+1 + z−k−1

= Tk+1(J(z))

SLIDE 60

COMRADE/COLLEAGUE MATRIX

SLIDE 61

We first normalize:
We will show that the eigenvalues of the companion matrix are the roots of p
The comrade matrix is defined as

       1

1 2 1 2

... ... ...

1 2 1 2

−ˇ p0 · · · −ˇ pn−3

1 2 − ˇ

pn−2 −ˇ pn−1        f(x) 2 ˇ fn ≈ p(x) =

n−1

k=0

ˇ pkTk(x) + 1 2Tn(x)

SLIDE 62

Given x which is a root of p, we have

       1

1 2 1 2

... ... ...

1 2 1 2

−ˇ p0 · · · −ˇ pn−3

1 2 − ˇ

pn−2 −ˇ pn−1               T0(x) T1(x) . . . Tn−2(x) Tn−1(x)        =        T1(x)

1 2(T0(x) + T2(x))

. . .

1 2(Tn−3(x) + Tn−1(x))

−ˇ p0T0(x) − · · · − ˇ pn−1Tn−1(x) + 1

2Tn−2(x)

        

SLIDE 63

Given x which is a root of p, we have

       1

1 2 1 2

... ... ...

1 2 1 2

−ˇ p0 · · · −ˇ pn−3

1 2 − ˇ

pn−2 −ˇ pn−1               T0(x) T1(x) . . . Tn−2(x) Tn−1(x)        =        T1(x)

1 2(T0(x) + T2(x))

. . .

1 2(Tn−3(x) + Tn−1(x))

−ˇ p0T0(x) − · · · − ˇ pn−1Tn−1(x) + 1

2Tn−2(x)

       =        xT0(x) xT1(x) . . . xTn−2(x)

1 2Tn(x) + 1 2Tn−2(x)

        

SLIDE 64

Given x which is a root of p, we have

       1

1 2 1 2

... ... ...

1 2 1 2

−ˇ p0 · · · −ˇ pn−3

1 2 − ˇ

pn−2 −ˇ pn−1               T0(x) T1(x) . . . Tn−2(x) Tn−1(x)        =        T1(x)

1 2(T0(x) + T2(x))

. . .

1 2(Tn−3(x) + Tn−1(x))

−ˇ p0T0(x) − · · · − ˇ pn−1Tn−1(x) + 1

2Tn−2(x)

       =        xT0(x) xT1(x) . . . xTn−2(x)

1 2Tn(x) + 1 2Tn−2(x)

       = x        T0(x) T1(x) . . . Tn−2(x) Tn−1(x)       

SLIDE 65

1.0
0.5

0.5 1.0

0.4
0.2

0.2 0.4

f(x) = ex cos 100x2

Numerical roots of p

1.0
0.5

0.5 1.0

2
1

1 2

1.0
0.5

0.5 1.0

2
1

1 2