Computational Fourier Analysis Mathematics, Computing and Nonlinear - PowerPoint PPT Presentation

Computational Fourier Analysis Mathematics, Computing and Nonlinear Oscillations Rubin H Landau Sally Haerer, Producer-Director Based on A Survey of Computational Physics by Landau, Páez, & Bordeianu with Support from the National Science Foundation Course: Computational Physics II 1 / 1

Outline 2 / 1

Applied Math: Approximate Fourier Integral Numerical Integration � + ∞ dt y ( t ) e − i ω t Transform Y ( ω ) = √ (1) 2 π −∞ N h y ( t i ) e − i ω t i � ≃ √ (2) 2 π i = 0 Approximate Fourier integral → finite Fourier series Consequences to follow 3 / 1

Experimental Constraints Too! � + ∞ dt y ( t ) e − i ω t Transform, Spectral: Y ( ω ) = √ (3) 2 π −∞ � + ∞ d ω Y ( ω ) e + i ω t Inverse, Synthesis: y ( t ) = √ (4) 2 π −∞ Real World: Data Restrict Us Measured y ( t ) only @ N times ( t i ’s) Discrete not continuous & not −∞ ≤ t ≤ + ∞ Can’t measure enough data to determine Y ( ω ) The inverse problem with incomplete data DFT: one possible solution 4 / 1

–0.5 0 t฀=฀0 t฀=฀T –1.0 0.0 –6 –4 –2 2 1.0 4 6 t 12 8฀ 10 14฀ 0.5 Algorithm with Discrete & Finite Times Measure: N signal values y y y y y 0 1 2 3 h N Uniform time steps ∆ t = h , t k = kh y k = y ( t k ) , k = 0 , 1 , . . . , N Finite T ⇒ ambiguity Integrate over all t; y ( t < 0 ) , y ( t > T ) = ? Assume periodicity y ( t + T ) = y ( t ) (removes ambiguity) ⇒ Y ( ω ) at N discrete ω i ’s ⇒ y 0 ≡ y N repeats! ⇒ N + 1 values, N independent 5 / 1

“A” Solution to Indeterminant Problem Discrete y ( t i ) s ⇒ Discrete ω i N independent y ( t i ) measured 1 ) ⇒ N independent Y ( ω i ) ( Y 0 ω i = ? Total time T = Nh ⇒ min ω 1 : -1 0 20 ω ω 1 = 2 π T = 2 π (5) Nh Only N ω i s (+1) ω k = k ω 1 , k = 0 , 1 , . . . , N (6) ω 0 = 0 × ω 1 = 0 (DC) (7) 6 / 1

Algorithm: Discrete Fourier Transform (DFT) f ( t ) dt ≃ � h f ( t i ) � Trapezoid Rule: � T � + ∞ dt e − i ω n t dt e − i ω n t Y ( ω n ) = √ y ( t ) ≃ √ y ( t ) (8) 2 π 2 π 0 −∞ N h y ( t k ) e − i ω n t k � √ (9) ≃ 2 π k = 1 Symmetrize notation, substitute ω n : N e − 2 π ikn / N Y n = 1 � hY ( ω n ) = y k (10) √ 2 π k = 1 7 / 1

14฀ –2 –0.5 0.0 0.5 1.0 –6 –4 0 ω 2 4 6 t 12 8฀ 10 –1.0 Discrete Inverse Transform (Function Synthesis) Trapezoid Rule for Inverse Frequency Step ∆ ω = ω 1 = 2 π T = 2 π Nh � + ∞ N d ω e i ω t e i ω n t 2 π � y ( t )= √ Y ( ω ) ≃ √ Y ( ω n ) (11) Nh 2 π 2 π −∞ n = 1 Trig functions ⇒ periodicity for y and Y : y ( t k + N ) = y ( t k ) , Y ( ω n + N ) = Y ( ω n ) (12) 1 Y( ) 0 0 20 40 8 / 1

0.0 0 –0.5 0.5 1.0 –6 –4 –2 2 ω 4 6 t 12 8฀ 10 14฀ –1.0 Consequences of Discreteness N N h y ( t k ) e − i ω n t k e i ω n t 2 π � � Y ( ω n ) ≃ √ y ( t ) ≃ √ Y ( ω n ) Nh 2 π 2 π k = 1 n = 1 1 Y( ) 0 0 20 40 ∆ ω = ω 1 = 2 π T = 2 π Ethical question Nh Finer ω ⇔ larger T = Nh Synthetic y ( t ) = bad t → T Finer ω → smoother Y ( ω ) Periodicity ↓ as T → ∞ “Pad” y ( t ) ⇒ smoother Aliases and ghosts: Y ( ω ) special lecture 9 / 1

Concise & Efficient DFT Computation Compute 1 Complex number Z = e − 2 π i / N = cos 2 π N − i sin 2 π (13) N N 1 � Z nk y k Y n = √ (Transform) (14) 2 π k = 1 √ N 2 π � Z − nk Y n (Synthesis, TF − 1 ) y k = (15) N n = 1 Z rotates y into complex Y and visa versa Compute only powers of Z (basis of FFT) Z nk ≡ ( Z n ) k = cos 2 π nk − i sin 2 π nk (16) N wN 10 / 1

DFT Fourier Series: y ( t ) = � ∞ n = 0 a n cos ( n ω t ) ≃ ? Discrete: Sample y ( t ) @ N times y ( t = t k ) ≡ y k , k = 0 , 1 , . . . , N (17) Repeat period T = Nh ⇒ y 0 = y N N independent a n s ⇒ Use trapezoid rule, ω 1 = 2 π/ Nh : � T N a n = 2 dt cos ( k ω t ) y ( t ) ≃ 2 � cos ( nk ω 1 ) y k (18) T N 0 k = 1 Truncate sum: N � y ( t ) ≃ a n cos ( n ω 1 t ) (19) n = 1 11 / 1

Code Implementation: DFTcomplex.py Example N = 1000; twopi = 2.*math.pi; sq2pi = 1./math.sqrt(twopi); def fourier(dftz): for n in range(0, N): zsum = complex(0.0, 0.0) for k in range(0, N): zexpo = complex(0, twopi*k*n/N) zsum += signal[k]*exp(-zexpo) dftz[n] = zsum * sq2pi 12 / 1

Assessment: Test Where know Answers Simple Analytic Cases y ( t ) = 3 cos ( ω t ) + 2 cos ( 3 ω t ) + cos ( 5 ω t ) e.g. (20) Known: Y 1 : Y 2 : Y 3 = 3 :2 :1 (9 :4 :1 power spectrum) 1 Resum to input signal? (> graph, idea of error) 2 Effect of: time step h , period T = Nh 3 Sample mixed signal: 4 y ( t ) = 4 + 5 sin ( ω t + 7 ) + 2 cos ( 3 ω t ) + sin ( 5 ω t ) (21) Effects of varying 4 & 7? 5 13 / 1

Physics Assessment Harmonic Anharmonic V(x) Harmonic p = 2 p = 6 V V Anharmonic x Linear Unbound x x Nonlinear Linear Nonlinear Determine Spectrum and Check Inversion Nonlinearly perturbed oscillator: 1 V ( x ) = 1 2 kx 2 � 1 − 2 � 3 α x (1) Determine when > 10 % higher harmonics ( b n > 1 ≥ 10 % ) 2 Highly nonlinear oscillator: 3 V ( x ) = kx 12 (2) Compare to sawtooth. 4 14 / 1

Summary Represent periodic or nonperiodic functions with DFT. Finiteness of measurements → ambiguities ( T ) Infinite series or integral not practical algorithm or in experiment. Approximate integration → simplicity & approximations Better high frequency components: smaller h , same T . Smoother transform: larger T , same h (padding). Less periodicity: more measurements. DFT is simple, elegant and powerful. Rotation between signal and transform space. ( e i φ ) n → Fast Fourier Transform. 15 / 1