Environmental Data Analysis with MatLab

1. Environmental Data Analysis with MatLab Lecture 9: • Fourier Series

2. SYLLABUS Lecture 01 Using MatLabLecture 02 Looking At DataLecture 03Probability and Measurement ErrorLecture 04 Multivariate DistributionsLecture 05Linear ModelsLecture 06 The Principle of Least SquaresLecture 07 Prior InformationLecture 08 Solving Generalized Least Squares ProblemsLecture 09 Fourier SeriesLecture 10 Complex Fourier SeriesLecture 11 Lessons Learned from the Fourier Transform Lecture 12 Power SpectraLecture 13 Filter Theory Lecture 14 Applications of Filters Lecture 15 Factor Analysis Lecture 16 Orthogonal functions Lecture 17 Covariance and AutocorrelationLecture 18 Cross-correlationLecture 19 Smoothing, Correlation and SpectraLecture 20 Coherence; Tapering and Spectral Analysis Lecture 21 InterpolationLecture 22 Hypothesis testing Lecture 23 Hypothesis Testing continued; F-TestsLecture 24 Confidence Limits of Spectra, Bootstraps

3. temporal periodicitiesand their periods astronomical rotation daily revolution yearly other natural ocean waves a few seconds anthropogenic electric power 60 Hz

4. spatial periodicitiesand their wavelengths natural sand dunes hundreds of meters • tree rings • a few millimeters anthropogenic furrows plowed in a field few tens of cm

5. Air temperature Black Rock Forest 365 days 1 year

6. Air temperature Black Rock Forest 1 day time, days

7. Stream FlowNeuse River 365 days 1 year discharge, cfs time, days

8. lingo temporal f(t) = C cos{ 2π t / T } spatial f(x) = C cos{ 2π x / λ } amplitude, C amplitude, C period, T wavelength, λ frequency, f=1/T "cycles per (year,day,s)" - angular frequency, ω=2 π/T "radians per ..." wavenumber, k=2 π/ λ • f(t) = C cos(ωt) f(x) = C cos(kx)

9. sinusoidal oscillation f(t) = C cos{ 2π (t-t0) / T } amplitude, C d(t) period, T time, t delay, t0

10. pairing sines and cosinesto avoid using time delays

11. Fourier Serieslinear model containing nothing butsines and cosines

12. A’s and B’s aremodel parameters ω’s are auxiliary variables

13. Data series: spacing Dt, length T Mean of series in cos(0*t) termLowest frequency is 1 cycle/TFreq. spacing Dw: 1,2,3,... cycles/T Highest frequency is (1 cycle)/2Dt(cosine part only – a zigzag) Nyquist frequency

14. frequency 0 = mean offset 32 numbers in physical space  32 numbers in spectral space cos and sin of (2pt/32 x 1,2,3,...15) (1,2,3...15 cycles/record-length) cos(2pt/32 *16) cos(0t) cos(Δωt) sin(Δωt) cos(2Δω t) sin(2Δω t)

15. problem of aliasingfrequencies higher than Nyquist(that is, periods shorter than 2Dt)masquerade as low frequenciesexamplesolution:use averages over Dt, not instantaneous samples spaced Dt apart.

16. d1 (t) = cos(w1t), with w1=2Dw d1(t) time, t d2(t) = cos{w2t}, with w2=(2+N)Dw, d2(t) time, t

17. cos(ω t) has same shape as cos(-ω t) • and • sin(ω t) has same shape as sin(-ω t) • so really only the • 0 toωny • part of the ω-axis is unique symmetry of sines and cosines

18. cos(ω t) = cos(ω t +2p) = cos(ω t -2p) = ... sin(ω t) = sin(ω t +2p) = sin(ω t -2p) = ... equivalent points on the ω-axis • so really only the • 0 toωny • part of the ω-axis is unique w -wny 0 wny 2wny 3wny

19. power spectral density it is big at a frequency ω when when sine or cosine at that frequency has a large coefficient

20. alternatively, plot amplitude spectral density

21. Stream FlowNeuse River all interesting frequencies near origin amplitude spectral density frequency, cycles per day 365.2 d = 12 cycles per 4500d 60.0 days amplitude spectral density 182.6 days 5 cycles per 4500d 6 per 4500d • line is misleading! • psd is discrete • (should be bar plot) period, days 7 • can plot period, T=1/f instead

22. purpose of the lecture switch from Fourier Series containing sines and cosines to Fourier Series containing complex exponentials

23. purpose of the lecture switch from Fourier Series containing sin(ωt) and cos(ωt) to Fourier series containing exp(-iωt) and exp(+iωt)

24. review of complex numbers

25. imaginary unitisuch thati2 = -1

26. complex numbera = ar + iai imaginary part real part

27. adding complex numbersa = ar + iai b = br + i bic = a+b = (ar + iai)+ (br + i bi ) = (ar + br )+ i(ai+bi ) ci cr … just add real and imaginary parts, separately

28. subtracting complex numbersa = ar + iai b = br + i bic = a-b = (ar + iai)-(br + i bi ) = (ar - br )+ i(ai- bi ) ci cr … just subtract real and imaginary parts, separately

29. multiplying complex numbersa = ar + iai b = br + i bic = ab = (ar + iai)(br + i bi ) == arbr + iar bi + iaibr + i2aibi =(arbr - aibi )+ i(ar bi + aibr ) cr ci … like multiplying polynomials

30. complex conjugate, a*a = ar + iai a* = ar - iai

31. absolute value, |a ||a|= [ ar2 + ai2 ]½note|a|2 = a* a

32. MatLabhandles complex numbers completely transparentlya = 2 + 3*i;b = 4 + 6*i;c = a+b;works just fine

33. Warning!accidentally resetting i to somethingother than iis so easyi=100; (and then you get nonsense) so execute aclear i;at the top of your script if you plan to use i

34. or use the alternate notationa = complex(2,3);b = complex(4,6);c = a+b;which is safer

35. end of review

36. Euler’s Formulaexp(iz) = cos(z) + i sin(z) … where does that come from ???

37. any complex number can be writtenz = r exp(iθ) where r = |z | and θ=tan-1(zi/zr)

38. Euler’s Formulascomplex exponentials can be written as sines and cosines

39. or reverse themsine and cosines can be written as complex exponentials

40. so a Fourier Seriesalternatively can be writtenas a sum of sines and cosinesor a sum of complex exponentials

41. old-style Fourier Series non-negative frequencies only, from 0 to ωny sin=0 for 0 and ωny

42. new-style Fourier Seriesor “Inverse Discrete Fourier Transform” added for compatibility with MatLab non-negative frequencies negative frequencies

43. why the weird ordering of frequencies? • ωn = ( 0, Δω, 2Δω, …,½N Δω, -(½N-1) Δω, …, -2 Δω, -Δω ) • same as • ωn = ( 0, Δω, 2Δω, …,½N Δω, (½N+1) Δω, …, (N-1)Δω, NΔω ) non-negative frequencies negative frequencies non-negative frequencies up to Nyquist non-negative frequencies above the Nyquist

44. least-squares solution for the Fourier coefficients, Cnor “Discrete Fourier Transform” … derivation requires complex version of least-squares. See text for details …

45. MatLabFourier Coefficients Cj from time series dnc = fft(d); vector of N data vector of N complex Fourier coefficients

46. MatLabtime series dnfrom Fourier Coefficients Cjd = ifft(c); vector of N complex Fourier coefficients vector of N data

47. standard setup M=N; tmax=Dt*(N-1); t=Dt*[0:N-1]'; fmax=1/(2.0*Dt); df=fmax/(N/2); f=df*[0:N/2,-N/2+1:-1]'; Nf=N/2+1; dw=2*pi*df; w=dw*[0:N/2,-N/2+1:-1]'; Nw=Nf;

48. standard setup same number M of Fourier Coefficients as number N of data M=N; tmax=Dt*(N-1); t=Dt*[0:N-1]'; fmax=1/(2.0*Dt); df=fmax/(N/2); f=df*[0:N/2,-N/2+1:-1]'; Nf=N/2+1; dw=2*pi*df; w=dw*[0:N/2,-N/2+1:-1]'; Nw=Nf;

49. standard setup maximum time, for N data sampled at Dt M=N; tmax=Dt*(N-1); t=Dt*[0:N-1]'; fmax=1/(2.0*Dt); df=fmax/(N/2); f=df*[0:N/2,-N/2+1:-1]'; Nf=N/2+1; dw=2*pi*df; w=dw*[0:N/2,-N/2+1:-1]'; Nw=Nf;

50. standard setup M=N; tmax=Dt*(N-1); t=Dt*[0:N-1]'; fmax=1/(2.0*Dt); df=fmax/(N/2); f=df*[0:N/2,-N/2+1:-1]'; Nf=N/2+1; dw=2*pi*df; w=dw*[0:N/2,-N/2+1:-1]'; Nw=Nf; time column-vector