Create Presentation
Download Presentation

Download Presentation
## Environmental Data Analysis with MatLab

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -

**Environmental Data Analysis with MatLab**Lecture 9: • Fourier Series**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**temporal periodicitiesand their periods**astronomical rotation daily revolution yearly other natural ocean waves a few seconds anthropogenic electric power 60 Hz**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**Air temperature**Black Rock Forest 365 days 1 year**Air temperature**Black Rock Forest 1 day time, days**Stream FlowNeuse River**365 days 1 year discharge, cfs time, days**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)**sinusoidal oscillation**f(t) = C cos{ 2π (t-t0) / T } amplitude, C d(t) period, T time, t delay, t0**Fourier Serieslinear model containing nothing butsines and**cosines**A’s and B’s aremodel parameters**ω’s are auxiliary variables**Data series: spacing Dt, length T Mean of series in**cos(0*t) termLowest frequency is 1 cycle/TFreq. spacing Dw: 1,2,3,... cycles/T Highest frequency is (1 cycle)/2Dt(cosine part only – a zigzag) Nyquist frequency**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)**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.**d1 (t) = cos(w1t), with w1=2Dw**d1(t) time, t d2(t) = cos{w2t}, with w2=(2+N)Dw, d2(t) time, t**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**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**power spectral density**it is big at a frequency ω when when sine or cosine at that frequency has a large coefficient**alternatively, plot**amplitude spectral density**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**purpose of the lecture**switch from Fourier Series containing sines and cosines to Fourier Series containing complex exponentials**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)**complex numbera = ar + iai**imaginary part real part**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**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**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**MatLabhandles complex numbers completely transparentlya = 2**+ 3*i;b = 4 + 6*i;c = a+b;works just fine**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**or use the alternate notationa = complex(2,3);b =**complex(4,6);c = a+b;which is safer**Euler’s Formulaexp(iz) = cos(z) + i sin(z)**… where does that come from ???**any complex number can be writtenz = r exp(iθ) where r = |z**| and θ=tan-1(zi/zr)**Euler’s Formulascomplex exponentials can be written as**sines and cosines**or reverse themsine and cosines can be written as complex**exponentials**so a Fourier Seriesalternatively can be writtenas a sum of**sines and cosinesor a sum of complex exponentials**old-style Fourier Series**non-negative frequencies only, from 0 to ωny sin=0 for 0 and ωny**new-style Fourier Seriesor “Inverse Discrete Fourier**Transform” added for compatibility with MatLab non-negative frequencies negative frequencies**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**least-squares solution for the Fourier coefficients, Cnor**“Discrete Fourier Transform” … derivation requires complex version of least-squares. See text for details …**MatLabFourier Coefficients Cj from time series dnc =**fft(d); vector of N data vector of N complex Fourier coefficients**MatLabtime series dnfrom Fourier Coefficients Cjd = ifft(c);**vector of N complex Fourier coefficients vector of N data**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;**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;**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;**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