Lecture 15 Orthogonal Functions Fourier Series LGA mean

Lecture 15 Orthogonal Functions Fourier Series

LGA mean daily temperature time series is there a global warming signal?

Model that includes annual variability T(t) = a + bt + A 1 cos(2 pf 1 t) + B 1 sin(2 pf 1 t) + A 2 cos(2 pf 2 t) + B 2 sin(2 pf 2 t) + A 3 cos(2 pf 3 t) + B 3 sin(2 pf 3 t) + … with f 1 = 1 cycle per year f 2 = 2 cycles per year etc

Why both sines and cosines? cos{2 pf 1(t-t 0)}

Why both sines and cosines? cos{2 pf 1(t-t 0)} Cosine does not ‘start’ at t=0 But remember cos(a+b)=cos(a)cos(b)-sin(a)sin(b)

cos(a+b)=cos(a) cos(b) - sin(a) sin(b) cos{2 pf 1(t-t 0)} = cos(2 pf 1 t 0) cos(2 pf 1 t) – sin(2 pf 1 t 0) sin(2 pf 1 t) = A cos(2 pf 1 t) + B sin(2 pf 1 t)

cos(a+b)=cos(a) cos(b) - sin(a) sin(b) cos{2 pf 1(t-t 0)} = cos(2 pf 1 t 0) cos(2 pf 1 t) – sin(2 pf 1 t 0) sin(2 pf 1 t) = A cos(2 pf 1 t) + B sin(2 pf 1 t) So using both sines and cosines moves the delay, t 0, out of the cosine, and into the coefficients of the sines and cosines. This trick ‘linearizes’ the unknown, t 0.

Why more than one frequency? f 1 = 1 cycle per year f 2 = 2 cycles per year etc Allows us to represent non-sinusoidal shape of annual cycle.

cos(ft) 0. 3 cos(2 ft) sum: cos(ft)+0. 3 cos(2 ft) exactly periodic, but shape not exactly sinusoidal

Least-squares fit to LGA data (up to f 8) data fit constant term, a error of fit, e linear term, bt

Statistics of linear term, bt b = 0. 31 degrees F per decade sd = [ e. Te / N ]1/2 = 7 deg F Cm = sd 2 [GTG]-1 sb = [ sd 2 Cmb, b ]1/2 = 0. 05 degrees F per decade 95% confidence b = 0. 31± 0. 1 degrees F per decade So LGA is warming

sines and cosines are “orthogonal” functions T(t) = A 0 + A 1 cos(2 pf 1 t) + B 1 sin(2 pf 1 t) + A 2 cos(2 pf 2 t) + B 2 sin(2 pf 2 t) + A 3 cos(2 pf 3 t) + B 3 sin(2 pf 3 t) + … Called a “Fourier Series” with f 2=2 f 1, f 3=3 f 1, etc

Standard least-squares G matrix 1 cos(2 pf 1 t 1) sin(2 pf 1 t 1) cos(2 pf 2 t 1) sin(2 pf 2 t 1) … 1 cos(2 pf 1 t 2) sin(2 pf 1 t 2) cos(2 pf 2 t 2) sin(2 pf 2 t 2) … G= 1 cos(2 pf 1 t 3) sin(2 pf 1 t 3) cos(2 pf 2 t 3) sin(2 pf 2 t 3) … 1 cos(2 pf 1 t 4) sin(2 pf 1 t 4) cos(2 pf 2 t 4) sin(2 pf 2 t 4) … 1 cos(2 pf 1 t 5) sin(2 pf 1 t 5) cos(2 pf 2 t 5) sin(2 pf 2 t 5) … 1 cos(2 pf 1 t 6) sin(2 pf 1 t 6) cos(2 pf 2 t 6) sin(2 pf 2 t 6) … …

With the proper choice of f 1 the matrix GTG is diagonal dot product of any pair of columns of G is zero columns of G are orthogonal

The proper choice of f 1 Suppose the time-series is N data points long, with spacing Dt. Then the lowest frequency must be f 1 = 1 / (NDt) one oscillation over the length of the time-series And the highest frequency must be f. N/2 = 1 / (2 Dt) one-half oscillation per sampling interval

f 1 = 1 / (2 NDt) f 1 = 1 / (2 Dt) note sine is zero

Count of unknowns The constant term, one unknown plus 2 coefficients per frequency, N/2 frequencies so N unknowns minus One unknown since the f. N/2 term, which has no sine term equals N unknowns, same as number of data

Mat. Lab Code N = 100; dt = 0. 5; tmin = 0. 0; t = tmin + dt*[0: N-1]'; tmax = tmin + dt*(N-1); % times vector df = 1/(N*dt); M = N; % frequency spacing % number of unknowns same as data G = zeros(N, M); G(: , 1)=ones(N, 1); for p = 2*[1: M/2 -1] G(: , p) = cos(pi*p*df*t); G(: , p+1) = sin(pi*p*df*t); end p=M/2; G(: , M) = cos(2*pi*p*df*t); % set up least-squares G matrix

GTG for N=100 [GTG]11= [GTG]NN=N Other diagonal elements [GTG]ii=N/2 Off diagonal elements are zero

So least-squares solution is m = [GTG]-1 GTd = = diag( N-1, 2/N, … 2/N, N-1 ) GT d NO matrix inversion required!

Example: Neuse River Hydrograph (100 days) GTG for N=100 data, d d=Gm with m=[GTG]-1 GTd d=Gm with m=DGTd where D=diag( N-1, 2/N, … 2/N, N-1 )

“spectrum” amount of power at different frequencies 2 si = 2 Ai 2 si fp + 2 Bi time-series has a lot of energy at frequency fp fi

Spectrum of Neuse data set for N=4380

2 mo 3 mo 4 mo 6 mo 12 mo Close up of low frequencies Big annual cycle in Neuse hydrograph

Error Estimates for Fourier Series Assume uncorrelated, normally-distributed data, d, with variance sd 2 The problem Gm=d is linear, so the unknowns, m, (the coefficients of the cosines and sines, Ai and Bi) are also normally-distributed. Since sines and cosines are orthogonal, GTG is diagonal and Cm= sd 2 [GTG]-1 is diagonal, too So that m’s have uncorrelated errors. All but the first and last have variance sm 2= 2 sd 2/N. The spectrum si 2=Ai 2+Bi 2 is the sum of two uncorrelated, normally distributed random variables and is thus c 22 distributed. The c 22 -distribution has a variance of 4, so that ss 2= 8 sd 2/N

Switching to complex numbers nothing different in principle but calculations become easier

But first Lets switch to angular frequency measured in radians per second wi = 2 p fi Beats writing all those 2 p’s !

Remember Euler’s formula exp( iwt ) = cos( wt ) + i sin( wt ) ?

exp( iwt ) = cos( wt ) + i sin( wt ) exp( -iwt ) = cos( wt ) - i sin( wt ) cos( wt ) = (1/2) [exp( iwt ) + exp( -iwt )] sin( wt ) = (1/2 i) [exp( iwt ) - exp( -iwt )]

Let’s compare =1 =0 T(t) = A 0 cos(w 0 t) + B 0 sin(w 0 t) + with wp=p. Dw A 1 cos(w 1 t) + B 1 sin(w 1 t) + w-p= -wp A 2 cos(w 2 t) + B 2 sin(w 2 t) + A 3 cos(w 3 t) + B 3 sin(w 3 t) + … with First, if T is real, then we must have C-p = Cp* T(t) =. . . + C-2 exp(-iw 2 t) + C-1 exp(-w 1 t) + C 0 exp(iw 0 t) + C 1 exp(iw 1 t) + C 2 exp(iw 2 t) + C 3 exp(iw 3 t) + … (Cpr-i. Cpi) [cos(wpt) - i sin(wpt)] + Then C-p exp(-wpt) + Cp exp(wpt) = (Cpr+i. Cpi) [cos(wpt) + i sin(wpt)] = 2 Cpr cos(wpt) - 2 Cpi sin(-w-pt)] So Ap= 2 Cpr and Bp= -2 Cpi So these two representations are equivalent

T(t) =. . . + C-2 exp(-iw 2 t) + C-1 exp(-w 1 t) + C 0 exp(iw 0 t) + C 1 exp(iw 1 t) + C 2 exp(iw 2 t) + C 3 exp(iw 3 t) + … T 0 T 1 T 2 T 3 T 4 … = … exp(-iw 2 t 0) … exp(-iw 2 t 1) … exp(-iw 2 t 2) … exp(-iw 2 t 3) Implies a simple form of the equation d=Gm exp(-iw 1 t 0) exp(-iw 1 t 1) exp(-iw 1 t 2) exp(-iw 1 t 3) exp(iw 0 t 0) exp(iw 0 t 1) exp(iw 0 t 2) exp(iw 0 t 3) exp( iw 1 t 0) exp( iw 1 t 1) exp( iw 1 t 2) exp( iw 1 t 3) exp(iw 2 t 0) … exp( iw 2 t 1) … exp( iw 2 t 2) … exp( iw 2 t 3) … … exp(-iw 2 t 4) exp(-iw 1 t 4) exp(iw 0 t 4) exp( iw 1 t 4) exp( iw 2 t 4) … … C-2 C-1 C 0 C 1 C 2 …

Least-squares with complex numbers real numbers: complex nos: given Gm =d minimize E=e. Te implies m=[GTG]-1 GT d The Hermitian transpose, that is, the transpose of the complex conjugate. given Gm =d minimize E=e. He where e. H = e*T implies m=[GHG]-1 GH d The formula m=[GHG]-1 GHd is not hard to work out using the standard minimization procedure, but we don’t have time to work it out in class.

d=Gm T 0 T 1 T 2 T 3 T 4 … = … exp(-iw 2 t 0) … exp(-iw 2 t 1) … exp(-iw 2 t 2) … exp(-iw 2 t 3) exp(-iw 1 t 0) exp(-iw 1 t 1) exp(-iw 1 t 2) exp(-iw 1 t 3) exp(iw 0 t 0) exp(iw 1 t 0) exp(iw 0 t 1) exp(iw 1 t 1) exp(iw 0 t 2) exp(iw 1 t 2) exp(iw 0 t 3) exp(iw 1 t 3) exp(iw 2 t 0) exp(iw 2 t 1) exp(iw 2 t 2) exp(iw 2 t 3) … … … exp(-iw 2 t 4) exp(-iw 1 t 4) exp(iw 0 t 4) exp(iw 1 t 4) exp(iw 2 t 4) … … C-2 C-1 C 0 C 1 C 2 … Note T 2 Si Ci exp(+iwit 2) … C-2 C-1 C 0 C 1 C 2 … m=N-1 GHm … exp(iw 2 t 0) exp(iw 2 t 1) exp(iw 2 t 2) … exp(iw 1 t 0) exp(iw 1 t 1) exp(iw 1 t 2) =N-1 … exp(iw 0 t 0) exp(iw 0 t 1) exp(iw 0 t 2) … exp(-iw 1 t 0) exp(-iw 1 t 1) exp(-iw 1 t 2) exp(iw 2 t 3) exp(iw 2 t 4) … exp(iw 1 t 3) exp(iw 1 t 4) … exp(iw 0 t 3) exp(iw 0 t 4) … exp(-iw 1 t 3) exp(-iw 1 t 4) … … exp(-iw 2 t 0) exp(-iw 2 t 1) exp(-iw 2 t 2) exp(-iw 2 t 3) exp(-iw 2 t 4) … Note C 2 Si Ti exp(-iw 2 ti) T 0 T 1 T 2 T 3 T 4 …

d=Gm T 0 T 1 T 2 T 3 T 4 … = … exp(-iw 2 t 0) … exp(-iw 2 t 1) … exp(-iw 2 t 2) … exp(-iw 2 t 3) exp(-iw 1 t 0) exp(-iw 1 t 1) exp(-iw 1 t 2) exp(-iw 1 t 3) exp(iw 0 t 0) exp(iw 1 t 0) exp(iw 0 t 1) exp(iw 1 t 1) exp(iw 0 t 2) exp(iw 1 t 2) exp(iw 0 t 3) exp(iw 1 t 3) exp(iw 2 t 0) exp(iw 2 t 1) exp(iw 2 t 2) exp(iw 2 t 3) … exp(-iw 2 t 4) exp(-iw 1 t 4) exp(iw 0 t 4) exp(iw 1 t 4) exp(iw 2 t 4) … Note T 2 Si Ci exp(+iwit 2) … C-2 C-1 C 0 C 1 C 2 … … … M=N-1 GHm … exp(iw 2 t 0) exp(iw 2 t 1) exp(iw 2 t 2) … exp(iw 1 t 0) exp(iw 1 t 1) exp(iw 1 t 2) =N-1 … exp(iw 0 t 0) exp(iw 0 t 1) exp(iw 0 t 2) … exp(-iw 1 t 0) exp(-iw 1 t 1) exp(-iw 1 t 2) Opposite signs exp(iw 2 t 3) exp(iw 2 t 4) … exp(iw 1 t 3) exp(iw 1 t 4) … exp(iw 0 t 3) exp(iw 0 t 4) … exp(-iw 1 t 3) exp(-iw 1 t 4) … … exp(-iw 2 t 0) exp(-iw 2 t 1) exp(-iw 2 t 2) exp(-iw 2 t 3) exp(-iw 2 t 4) … Note C 2 Si Ti exp(-iw 2 ti) … C-2 C-1 C 0 C 1 C 2 … T 0 T 1 T 2 T 3 T 4 …

Discrete Fourier Transform Find the coefficients C given the data, T Equivalent to m = GHd Note normalization factor of N-1 has been omitted Ck = Sn=-N/2 N/2 Tn exp(± 2 pikn/N ) with k=-½N, …, ½N Discrete Inverse Fourier Transform Find the data T given the coefficients, C Equivalent to d = N-1 Gm Note normalization factor of N-1 has been added Tn = N-1 Sk=-N/2 N/2 Ck exp( 2 pikn/N ) with n=-½N, …, ½N ± Warnings: 1) no one can agree on signs 2) no one can agree on normalizations

Counting unknowns frequencies from –(N/2)Dw to (N/2)Dw in steps of Dw So N+1 complex numbers, Cp So 2 N+2 real and imaginary parts, Cpr and Cpi But C-p = Cp*, so really only N/2+1 unknown complex numbers So N+2 real and imaginary parts , Cpr and Cpi (p 0) But C 0 i=0 and CN/2 i=0 (always) So N unknowns, matching N data

% standard fft setup. The standard implementation of the digital fourier % transform is VERY INFLEXIBLE. Learn these rules: N=256; % you can choose the length N of the time series % in some implementations N can be any positive % integer, but in others it MUST be a % power-or-two. I set it here to 256, which % is two-to-the-eigth-power. dt=1. 0; % and you can choose the sampling interval dt % but then the following variables are set tmax=dt*(N-1); % we presume the time series starts at t=0, so % the maximum time is tmax t=dt*[0: N-1]; % time then goes from 0 to (N-1)*dt fmax=1/(2. 0*dt); % the maximum frequency in the fft calculation is % called the Nyquist frequency. It is % determined by the two-points-per wavelength % rule df=fmax/(N/2); % the frequency spacing, df, assumes that a N-point % time series is reperesented by an N-point fourier % transform f=df*[0: N/2, -N/2+1: -1]'; % The fourier transform has N values, from a negative % frequency of -(fmax-df) through zero freqency, to % positive frequency of fmax. But note the weird order. The % zero and positive frequencies are put in the first % half of the array and the negative frequencies are % put in the second half.

% p is the timeseries whose transform is being computed w 0 = 2*pi*fmax/10; % sample p, a simple sinusoid of frequency w 0 p = sin(w 0*t); % fourier transform using Mat. Lab's fft function. The help function % says that it uses the NEGATIVE sign in the exponential. pt=fft(p); % these are the coefficients, C, of the complex exponential % presumably one would do something with the fourier transform % at this point - apply a filter, for example. But I do nothing. % Inverse fourier transform using the Matlab function ifft. % Help says it uses the POSITIVE sign in the exponential, and that % it has the right normalization that ifft(x))=x. But % BE WARNED, that doesn't mean that the normalization on fft % is 1 and that the normalization on ifft is (1/N), like % I had it in class. You can put any constant, b, in front % of the fft integral, as long as you put 1/b in front of % the IFFT integral. But judging by the Help, I think that % Matlab used b=1. pr=ifft(pt); % this reconstructs the function from the coefficients

The fast Fourier transform algorithm The Fourier Transform equation m = [GTG]-1 GTd = diag(N-1, 2/N, … 2/N, N-1) GT d has N multiplications for each of N unknowns, So N 2 in total. For example, for N=1024, N 2=1, 048, 576 But in the special case of N being a power of two, there is an algorithm – the fast fourier transform algorithm - that can compute m in only Nlog 2 N multiplications For example, N=1024, Nlog 2 N=10, 240 A substantial savings! Mat. Lab implements it. Use it!