12.006J Lectures 12–14: Spectral Analysis for Dynamical Systems

Lecture notes for 12.006J/18.353J/2.050J, Nonlinear Dynamics: Chaos Contents 1 Spectral analysis for dynamical systems 1.1 Fourier transforms . . . . . . . . . . . . . . . . . . . 1.1.1 Continuous Fourier transform . . . . . . . . . 1.1.2 Discrete-time signals . . . . . . . . . . . . . . 1.1.3 Discrete Fourier transform . . . . . . . . . . . 1.1.4 Inverse discrete Fourier tranform . . . . . . . 1.2 The autocorrelation function and the power spectrum 1.3 Power spectrum of a periodic signal . . . . . . . . . . 1.3.1 Sinusoidal signal . . . . . . . . . . . . . . . . 1.3.2 Non-sinusoidal signal . . . . . . . . . . . . . . 1.3.3 tmax /T 6= integer . . . . . . . . . . . . . . . . 1.3.4 Conclusion . . . . . . . . . . . . . . . . . . . . 1.4 Quasiperiodic signals . . . . . . . . . . . . . . . . . . 1.5 Aperiodic signals . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 2 2 3 5 7 10 10 11 12 15 15 19 Spectral analysis for dynamical systems Power spectra provide one of the most important tools for analyzing the behavior of dynamical systems, both theoretically and experimentally. 1.1 Fourier transforms The precise oscillatory nature of an observed time series x(t) is usually not identifiable from x(t) alone. We may ask 1 • How well-defined is the the dominant frequency of oscillation? • How many frequencies of oscillation are present? • What are the relative contributions of all frequencies? The analytic tool for answering these and myriad related questions is the Fourier transform. 1.1.1 Continuous Fourier transform We first state the Fourier transform for functions that are continuous with time. The Fourier transform of a function f (t) is Z ∞ 1 F (ω) = √ f (t)e−iωt dt 2π −∞ Similarly, the inverse Fourier transform is Z ∞ 1 f (t) = √ F (ω)eiωt dω. 2π −∞ That the second relation is the inverse of the first may be proven, but we save that calculation for the discrete transform, below. 1.1.2 Discrete-time signals We are interested in the analysis of observational or experimental data, which is almost always discrete. Thus we specialize to discrete Fourier transforms. In modern data, one almost always observes a discretized signal xj , j = {0, 1, 2, . . . , n − 1} 2 We take the sampling interval—the time between samples—to be ∆t. Then xj = x(j∆t). The discretization process is pictured as x x(t) ∆t j−1 j t j+1 A practical question concerns the choice of ∆t. To choose it, we must know the highest frequency, fmax , contained in x(t). The shortest period of oscillation is Tmin = 1/fmax Pictorially, Tmin x t We require at least two samples per period. Therefore Tmin 1 ∆t ≤ = . 2 2fmax 1.1.3 Discrete Fourier transform The discrete Fourier transform (DFT) of a time series xj , j = 0, 1, . . . , n − 1 is   n−1 X 2πjk x̂k = xj exp −i k = 0, 1, . . . , n − 1 n j=0 3 To gain some intuitive understanding, consider the range of the exponential multiplier. • k = 0 ⇒ exp(−i2πjk/n) = 1. Then x̂0 = X xj j Thus x̂0 is n times the mean of the xj ’s. This is the “DC” component of the transform. Question: Suppose a seismometer measures ground motion. What would x̂0 6= 0 mean? • k = n/2 ⇒ exp(−i2πjk/n) = exp(−iπj). Then X x̂n/2 = xj (−1)j j = x0 − x1 + x2 − x3 . . . Frequency index n/2 is clearly the highest accessible frequency. • The frequency indices k = 0, 1, . . . , n/2 correspond to frequencies fk = k/tmax , i.e., k oscillations per tmax , the period of observation. Index k = n/2 then corresponds to n  1  1 fmax = = 2 n∆t 2∆t But if n/2 is the highest frequency that the signal can carry, what is the significance of x̂k for k > n/2? For real xj , frequency indices k > n/2 are redundant, being related by x̂k = x̂∗n−k where z ∗ is the complex conjugate of z (i.e., if z = a + ib, z ∗ = a − ib). 4 We derive this relation as follows. From the definition of the DFT, we have   n−1 X 2πj(n − k) ∗ x̂n−k = xj exp +i n j=0 n−1 X  −i2πjk = xj exp (i2πj) exp | {z } n j=0  1 n−1 X  −i2πjk = xj exp n j=0  = x̂k where the + in the first equation derives from the complex conjugation, and the last line again employs the definition of the DFT. Note that we also have the relation x̂∗−k = x̂∗n−k = x̂k . The frequency indices k > n/2 are therefore sometimes referred to as negative frequencies 1.1.4 Inverse discrete Fourier tranform The inverse DFT is given by   n−1 1X 2πjk xj = x̂k exp +i n n j = 0, 1, . . . , n − 1 k=0 We proceed to demonstrate this inverse relation. 5 We begin by substituting the DFT for x̂k , using dummy variable j 0 :       n−1 X n−1 0 X 1 2πj k 2πkj   exp +i xj = xj 0 exp −i n n n 0 k=0 j =0   n−1 n−1 X 2πk(j 0 − j) 1X exp −i = xj 0 n 0 n j =0 k=0 n−1  1X xj 0 × = n 0 j =0 = n, j 0 = j 0, j 0 = 6 j 1 (nxj ) n = xj P The third relation derives from the fact that the previous k amounts to a vanishing sum over the unit circle in the complex plane, except when j 0 = j. To see why the sum over the circle vanishes, consider the example of j 0 − j = 1 and n = 4. The elements of the sum are then just the four points on the unit circle that intersect the real and imaginary axes, i.e.,   3 X 2πk(j 0 − j) exp −i = e0 + e−iπ/2 + e−iπ + e−i3π/2 4 k=0 = 1−i−1+i = 0. Finally, note that the DFT relations imply that xj is periodic in n, so that xj+n = xj . Consequently a finite time series is treated as if it were recurring: 6 x(t) −tmax 1.2 tmax 0 2 tmax t The autocorrelation function and the power spectrum Assume that the time series xj has zero mean and that it is periodic, i.e., xj+n = xj . Define the autocorrelation function ψ: ψm = n−1 X x∗j xj+m j=0 where ψm = ψ(m∆t) The autocorrelation function measures the degree to which a signal resembles itself over time. Thus it measures the predictability of the future from the past. To gain some intuition: • Consider, for example, m = 0 and real xj . Then ψ0 = n−1 X x2j , j=0 which is n times the mean squared value of xj . • Alternatively, if m∆t is much less than the dominant period of the data, ψm should not be too much less than ψ0 . • Last, if m∆t is much greater than the dominant period of the data, |ψm | is relatively small. 7 A typical ψm looks like Ψm m The power spectrum of a time series is the magnitude squared of its Fourier transform:  2 n−1 X 2πjk xj exp −i |x̂k |2 = . n j=0 The Wiener-Khintchin theorem states that power spectrum = Fourier transform of the autocorrelation. In symbols, n−1 X   2πkm |x̂k | = ψm exp −i n m=0 2 We also have the inverse relation   n−1 2πkm 1X |x̂k |2 exp +i ψm = n n k=0 To prove the latter relation, we first substitute the inverse DFT for xj and xj+m in the definition of ψm : ψm = n−1 X x∗j xj+m j=0 " n−1  # " X  # n−1 n−1 0 X X 1 2πkj 1 2πk (j + m) = x̂∗k exp −i x̂k0 exp i n n n 0 n j=0 k=0 k =0 8 We then change the order of the summations and simplify as follows:   n−1   n−1 n−1 2πmk 0 X 2πj(k 0 − k) 1 XX ∗ x̂k x̂k0 exp i exp i ψm = 2 n n n 0 j=0 k=0 k =0 | {z } 0 = n, k = k = 0, k 0 6= k   n−1 2πmk 1X ∗ x̂k x̂k exp i = n n k=0 which is the Wiener-Khintchin relation. By Fourier transforming ψm we also prove the inverse relation: the power spectrum is the Fourier transform of the autocorrelation. For a real time series {xj }, we can use the previously derived relation x̂∗k = x̂n−k = x̂−k to show that |x̂k |2 = x̂k x̂∗k = x̂k x̂n−k = x̂∗n−k x̂n−k = |x̂n−k |2 . This redundancy results from the fact that neither the autocorrelation nor the power spectrum contain information on any “phase lags” in either xj or its individual frequency components. Thus while the DFT of an n-point time series results in n independent quantities (2 ×n/2 complex numbers), the power spectrum yields only n/2 independent quantities. One may therefore show that there are an infinite number of time series that have the same power spectrum, but that each time series uniquely defines its Fourier transform, and vice-versa. Consequently a time series cannot be reconstructed from its power spectrum or autocorrelation function. 9 1.3 Power spectrum of a periodic signal Consider a periodic signal  2π x(t) = x(t + T ) = x t + ω  Consider the extreme case where the period T is equal to the duration of the signal: T = tmax = n4t The Fourier components are separated by ∆f = 1 tmax i.e. at frequencies 0, 1/T, 2/T, . . . , (n − 1)/T. 1.3.1 Sinusoidal signal In the simplest case, x(t) is a sine or cosine, i.e.,   2πt x(t) = sin . tmax What is the Fourier tranform? Pictorially, we expect xk 2 x(t) tmax t 1/T 10 k ∆ f = k/T We calculate the power spectrum analytically, beginning with the DFT:   X −i2πjk x̂k = xj exp n j     X −i2πjk 2πj∆t exp = sin t n max j        i2πj∆t 1 X −i2πj∆t −i2πjk exp = − exp exp 2i j tmax tmax n        ∆t k k 1 X ∆t exp i2πj − − exp −i2πj + = 2i j tmax n tmax n = ± Thus n 2i when k = ±n∆t . tmax n2 |x̂k | = 4 2 1.3.2 for k = ±1. Non-sinusoidal signal Consider now a non-sinusoidal yet periodic signal, similar to the relaxation oscillations seen in the van der Pol limit cycle. The non-sinusoidal character of such oscillations implies that it contains higher-order harmonics, i.e., integer multiples of the fundamental frequency 1/T . Thus, pictorially, we expect xk 2 x(t) tmax t 1/T 2/T 3/T 11 fundamental harmonics k ∆ f = k/T Now suppose tmax = pT , where p is an integer. The non-zero components of the power spectrum must still be at frequencies 1/T, 2/T, . . . . But since 1 1 tmax pT the frequency resolution is p times greater. Contributions to the power spectrum would remain at integer multiples of the frequency 1/T , but spaced p samples apart on the frequency axis. ∆f = 1.3.3 = tmax /T 6= integer If tmax /T is not an integer, the (effectively periodic) signal looks like x(t) tmax t We calculate the power spectrum of such a signal, assuming the sinusoidal function   2πt x(t) = exp i T which has the discrete form   2πj∆t xj = exp i . T The DFT is n−1 X     2πj∆t 2πjk x̂k = exp i exp −i . T n j=0 Set φk = ∆t k − . T n 12 Then x̂k = n−1 X exp (i2πφk j) . j=0 Recall the identity n−1 X xn − 1 x = . x − 1 j=0 Then x̂k = j exp(i2πφk n) − 1 . exp(i2πφk ) − 1 The power spectrum is |x̂k |2 = x̂k x̂∗k = 1 − cos(2πφk n) 1 − cos(2πφk ) sin2 (πφk n) = . sin2 (πφk ) Note that tmax n∆t −k = −k T T is the difference between a DFT index k and the “real” non-integral frequency index tmax /T . nφk = Assume that n is large and k is close to that “real” frequency index such that n∆t − k  n. T Consequently φk  1, so we may also assume nφk = πφk  1. Then 2 |x̂k | sin2 (πφk n) ' (πφk )2 sin2 (πφk n) = n (πφk n)2 2 sin2 z ∝ z2 13 where  z = nπφk = π n∆t −k T   =π  tmax −k . T Thus |x̂k |2 is no longer a simple spike. Instead, as a function of z = nπφk it appears as sin 2 z / z2 1 −3π −2π −π π 0 2π z=n πφ k 3π The plot gives the kth component of the power spectrum of ei2πt/T as a function of π(tmax /T − k). To interpret the plot, let k0 be the integer closest to tmax /T . There are then two extreme cases: 1. tmax is an integral multiple of T: tmax − k0 = 0. T The spectrum is perfectly sharp: xk 2 sin 2z/z k0 0 k 2. tmax /T falls midway between two frequencies. Then tmax 1 − k0 = . T 2 The spectrum is smeared: 14 z xk 2 sin 2z/z k0 0 k z The smear decays like 1 (k − tmax 1.3.4 /T )2 ∼ 1 k2 Conclusion The power spectrum of a periodic signal of period T is composed of: 1. a peak at the frequency 1/T 2. a smear (sidelobes) near 1/T 3. possibly harmonics (integer multiples) of 1/T 4. smears near the harmonics. 1.4 Quasiperiodic signals Let y be a function of r independent variables: y = y(t1 , t2 , . . . , tr ). y is periodic, of period 2π in each argument, if y(t1 , t2 , . . . , tj + 2π, . . . , tr ) = y(t1 , t2 , . . . , tj , . . . , tr ), j = 1, . . . , r y is called quasiperiodic if each tj varies with time at a different rate (i.e., different “clocks”). We have then tj = ωj t, j = 1, . . . , r. 15 The quasiperiodic function y has r fundamental frequencies: fj = ωj 2π and r periods Tj = 1 2π = . fj ωj Example: The astronomical position of a point on Earth’s surface changes due to • rotation of Earth about axis (T1 = 24 hours). • revolution of Earth around sun (T2 ' 365 days). At long time scales, we also have changes in precession (26 Kyr), obliquity (41 Kyr), and eccentricity (∼100 Kyr). Considering just two oscillations (e.g, rotation and revolution), we can conceive of such a function on a 2-D torus T 2 , existing in a 3-D space. T2 T1 Here we think of a disk spinning with period T1 while it revolves along the circular path with period T2 . Such behavior can be conceived as a trajectory on the surface of a doughnut or inner tube, or a torus T2 in R3 . 16 f2 f1 What is the power spectrum of a quasiperiodic signal x(t)? There are two possibilities: 1. The quasiperiodic signal is a linear combination of independent periodic functions. For example: x(t) = r X xi (ωi t). i=1 Because the Fourier transform is a linear transformation, the power spectrum of x(t) is a set of peaks at frequencies f1 = ω1 /2π, f2 = ω2 /2π, . . . and their harmonics m1 f1 , m2 f2 , . . . (m1 , m2 , . . . positive integers). 2. The quasiperiodic signal x(t) depends nonlinearly on periodic functions. For example, x(t) = sin(2πf1 t) sin(2πf2 t) = 1 1 cos(|f1 − f2 |2πt) − cos(|f1 + f2 |2πt). 2 2 The fundamental frequencies are |f1 − f2 | and |f1 + f2 |. The harmonics are m1 |f1 − f2 | and m2 |f1 + f2 |, 17 m1 , m2 positive integers. The nonlinear case requires more attention. In general, if x(t) depends nonlinearly on r periodic functions, then the harmonics are |m1 f1 + m2 f2 + . . . + mr fr |, mi arbitrary integers. In what follows, we specialize to r = 2 frequencies, and forget about finite ∆f . Each nonzero component of the spectrum of x(ω1 t, ω2 t) is a peak at f = |m1 f1 + m2 f2 |, m1 , m2 integers . There are two cases: 1. f1 /f2 rational ⇒ sparse spectrum. 2. f1 /f2 irrational ⇒ dense spectrum. To understand this, rewrite f as f = f2 m1 f1 + m2 . f2 In the rational case, f1 integer = . f2 integer Then m1 integer 1 f1 + m2 = + integer = integer multiple of . f2 f2 f2 Thus the peaks of the spectrum must be separated (i.e., sparse). Alternatively, if f1 /f2 is irrational, then m1 and m2 may always be chosen so that f1 m1 + m2 is not similarly restricted. f2 These distinctions have further implications. 18 In the rational case, n1 f1 = , f2 n2 n1 , n2 integers. Since n1 n2 = f1 f2 the quasiperiodic function is periodic with period T = n1 T1 = n2 T2 . All spectral peaks must then be harmonics of the fundamental frequency f0 = 1 f1 f2 = = . T n1 n2 Thus the rational quasiperiodic case is in fact periodic, and some writers restrict quasiperiodicity to the irrational case. Note further that, in the irrational case, the signal never exactly repeats itself. One may consider, as an example, the case of a child walking on a sidewalk, attempting with uniform steps to never step on a crack. Then if x(t) were the distance from the closest crack at each step, it would only be possible to avoid stepping on a crack if the ratio step size crack width were rational. 1.5 Aperiodic signals Aperiodic signals are neither periodic nor quasiperiodic. Aperiodic signals appear random, though they may have a deterministic foundation. 19 An example is white noise, which is a signal that is “new” and unpredictable at each instant, e.g., x(t) t Statistically, each sample of a white-noise signal is independent of the others, and therefore uncorrelated to them. The power spectrum of white noise is, on average, flat: xk 2 k The flat spectrum of white noise is a consequence of its lack of harmonic structure (i.e., one cannot recognize any particular tone, or dominant frequency). We proceed to derive the spectrum of a white noise signal x(t). Rather than considering only one white-noise signal, we consider an ensemble of such signals, i.e., x(1) (t), x(2) (t), . . . where the superscipt denotes the particular realization within the ensemble. Each realization is independent of the others. Now discretize each signal so that j = 0, . . . , n − 1 xj = x(j∆t), We take the signal to have finite length n but consider the ensemble to contain an infinite number of realizations. 20 We use angle brackets to denote ensemble averages. The ensemble-averaged mean of the jth sample is then p 1 X (i) hxj i = lim xj p→∞ p i=1 Similarly, the mean-square value of the jth sample is p 1 X  (i) 2 2 xj = lim xj p→∞ p i=1 Now assume stationarity: hxj i and x2j are independent of j. We take these mean values to be hxi and x2 , respectively, and assume hxi = 0. Recall the autocorrelation ψm : ψm = n−1 X xj xj+m . j=0 By definition, each sample of white noise is uncorrelated with its past and future. Therefore * + X hψm i = xj xj+m j = n x2 δm where  1 m=0 0 else We obtain the power spectrum from the autocorrelation function by the Wiener-Khintchine theorem:   n−1 X 2πmk |x̂k |2 = hψm i exp −i n m=0 δm = = n−1 X n x2 m=0 = n x2 = constant. 21  2πmk δm exp −i n  Thus for white noise, the spectrum is indeed flat, as previously indicated: xk 2 k A more common case is “colored” noise: a continuous spectrum, but not constant: xk 2 k In such (red) colored spectra, there is a relative lack of high frequencies. The signal is still apparently random, but only beyond some interval ∆t. The autocorrelation of colored noise is broader, e.g., Ψm m 22