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

Lesson 9 S IGNAL SMOOTHING AND ROOT FINDING

• 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

W HITE NOISE

• Return to the case of sampling a function at evenly spaced points on the periodic interval f ( θ 1 ) f ( θ m ) … - 3 - 2 - 1 0 1 2 3 • 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: r 1 … r m - 3 - 2 - 1 0 1 2 3

Instances of white noise m = 5 1.0 • Function values 0.5 - 3 - 2 - 1 1 2 3 - 0.5 - 1.0

Instances of white noise m = 10 1.0 • Function values 0.5 - 3 - 2 - 1 1 2 3 - 0.5 - 1.0

Instances of white noise m = 20 1.0 • Function values 0.5 - 3 - 2 - 1 1 2 3 - 0.5 - 1.0

Instances of white noise m = 40 1.0 • Function values 0.5 - 3 - 2 - 1 1 2 3 - 0.5 - 1.0

Instances of white noise m = 100 1.0 • Function values 0.5 - 3 - 2 - 1 1 2 3 - 0.5 - 1.0

White noise DFT coefficients m = 5 1.00 0.50 0.20 0.10 0.05 - 1 0 1 2 Computed Fourier coefficients

White noise DFT coefficients m = 10 1.00 0.70 0.50 0.30 0.20 0.15 0.10 - 4 - 2 0 2 4 Computed Fourier coefficients

White noise DFT coefficients m = 20 1.00 0.50 0.20 0.10 0.05 0.02 - 5 0 5 Computed Fourier coefficients

White noise DFT coefficients m = 40 1.00 0.50 0.20 0.10 0.05 0.02 0.01 - 10 0 10 Computed Fourier coefficients

White noise DFT coefficients m = 100 1.00 0.50 0.20 0.10 0.05 0.02 0.01 - 40 - 20 0 20 40 Computed Fourier coefficients

White noise DFT coefficients m = 1000 1.000 0.500 0.100 0.050 0.010 0.005 - 400 - 200 0 200 400 Computed Fourier coefficients

White noise DFT coefficients m = 10000 1.000 0.500 0.100 0.050 0.010 0.005 0.001 - 4000 - 2000 0 2000 4000 Computed Fourier coefficients

S IGNAL DENOISING

• 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

Original Signal Signal with noise m = 50 m = 50 3 3 2 2 1 1 - 3 - 2 - 1 1 2 3 - 3 - 2 - 1 1 2 3 - 1 - 1

DFT coefficients m = 50 Noisy signal 0.01 10 - 5 10 - 8 10 - 11 10 - 14 Original signal - 20 - 10 0 10 20

DFT coefficients m = 100 1 Noisy signal 0.001 10 - 6 10 - 9 10 - 12 10 - 15 Original signal - 20 - 10 0 10 20

DFT coefficients m = 1000 1 Noisy signal 0.001 10 - 6 10 - 9 10 - 12 10 - 15 Original signal - 20 - 10 0 10 20

DFT coefficients m = 10000 1 Noisy signal 10 - 4 10 - 8 10 - 12 10 - 16 Original signal - 20 - 10 0 10 20

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients m = 100 2.5 2.0 1.5 1.0 0.5 - 3 - 2 - 1 1 2 3

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients m = 1000 2.5 2.0 1.5 1.0 0.5 - 3 - 2 - 1 1 2 3

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients m = 10000 2.5 2.0 1.5 1.0 0.5 - 3 - 2 - 1 1 2 3

We denoise the signal by taking only the –10,…,10 computed coefficients of the noise signal, and throwing out all other coefficients m = 100000 2.5 2.0 1.5 1.0 0.5 - 3 - 2 - 1 1 2 3

T AYLOR SERIES ROOT FINDING

• Consider a Taylor series ∞ ˆ � f k z k f ( z ) = k =0 • 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: m − 1 m − 1 ˆ f ( z ) f k z k + z m = p k z k + z m � � ≈ p ( z ) = ˆ ˆ ˆ f m f m k =0 k =0

• We will rephrase the root finding problem as a eigenvalue problem • Consider the companion matrix   1 ...       1   − ˆ − ˆ − ˆ p 0 p 1 p m − 1 · · · • We will show that the eigenvalues of the companion matrix are the roots of p

• Given z , we have       1 1 z . . ...       . . . .  =             z m − 2 z m − 1 1      z m − 1 p m − 1 z m − 1 − ˆ − ˆ − ˆ − ˆ p 0 − ˆ p 1 z − · · · − ˆ p 0 p 1 p m − 1 · · · • If z is a root of p , then this simplifies:       1 z z . . .       . . . . . .  =  = z             z m − 1 z m − 1 z m − 2     p m − 1 z m − 1 z m − 1 − ˆ p 0 − ˆ p 1 z − · · · − ˆ z m Eigenvalue! Eigenvector!

• Thus if z 1 ,…, z m denote the m roots of p , we have   1   z 1 ...   ...  V = V         1  z m − ˆ p 0 − ˆ p 1 − ˆ p m − 1 · · · for   1 1 · · · . . ... . . V =   . .   z m − 1 z m − 1 · · · 1 m

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

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 4 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 8 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 16 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 32 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 40 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 50 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

1.0 0.5 Example: sin 4 x = - 1.0 - 0.5 0.5 1.0 - 0.5 m = 100 - 1.0 4 2 - 4 - 2 2 4 - 2 - 4

L AURENT SERIES ROOT FINDING

• 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 ( z ) ≈ z − α f α f β − 1 f α z α + · · · + ˆ z m − 2 + z m − 1 = p ( z ) � f β z β � ˆ = + · · · + ˆ ˆ ˆ ˆ f β f β f β f β

• 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 z 1 , . . . , z r of g that lie on the unit circle • Then − � ��� z 1 , . . . , − � ��� z r are roots of f on the periodic interval

f ( θ ) = cos(50 cos 2 θ ) cos θ f ( − i log z ) Numerical roots of 1.0 1.5 1.0 0.5 0.5 - 3 - 2 - 1 1 2 3 - 1.5 - 1.0 - 0.5 0.5 1.0 1.5 - 0.5 - 0.5 - 1.0 - 1.5 - 1.0

f ( θ ) = cos(50 cos 2 θ ) cos θ f ( − i log z ) Pruned numerical roots of 1.0 1.0 0.5 0.5 - 1.0 - 0.5 0.5 1.0 - 3 - 2 - 1 1 2 3 - 0.5 - 0.5 - 1.0 - 1.0

f ( θ ) = cos(50 cos 2 θ ) cos θ f ( − i log z ) Pruned numerical roots of 1.0 1.0 0.5 0.5 - 1.0 - 0.5 0.5 1.0 - 3 - 2 - 1 1 2 3 - 0.5 - 0.5 - 1.0 - 1.0 • Roots mapped to periodic interval

C HEBYSHEV SERIES ROOT FINDING

• We want to find the roots of ∞ ˇ � 0 = f ( x ) = f k T k ( x ) k =0 • The easy way: � Recall the Joukowsky map J ( z ) = 1 and let � z + 1 � z 2 − 1 ∞ g ( z ) = f ( J ( z )) = 1 f 0 + 1 f − k z k + ˇ ˇ ˇ � � f k z k 2 2 k = −∞ k =1 48 � Truncate g and find the roots z 1 , . . . , z r that lie on the lower half of the unit circle � The roots of f are then J ( z 1 ) , . . . , J ( z r ) • The problem: this requires a matrix of size 2 n × 2 n instead of n × n

T HREE - TERM RECURRENCE RELATIONSHIP