Environmental Data Analysis with MatLab

# Environmental Data Analysis with MatLab

## Environmental Data Analysis with MatLab

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -
##### Presentation Transcript

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. A’s and B’s aremodel parameters ω’s are auxiliary variables

12. 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

13. 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)

14. 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.

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

16. 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

17. 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

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

19. alternatively, plot amplitude spectral density

20. 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

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

22. 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)

23. review of complex numbers

24. imaginary unitisuch thati2 = -1

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

26. 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

27. 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

28. 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

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

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

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

32. 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

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

34. end of review

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

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

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

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

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

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

41. 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

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

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

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

45. 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;

46. 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;

47. 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;

48. 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