Fourier Analysis of Time Series by Dr. R. L. Herman, UNC Wilmington Friday, September 20, 2002 This is a work in progress.

Introduction Often one is interested in determining the frequency content of signals. Signals are typically represented as time dependent functions. Real signals are continuous, or analog signals. However, through sampling the signal by gathering data, the signal does not contain high frequencies and is finite in length. The data is then discrete and the corresponding frequencies are discrete and bounded. Thus, in the process of gathering data, one seriously affects the frequency content of the signal. This is true for simple a superposition of signals with fixed frequencies. The situation becomes more complicated if the data has an overall non-constant trend or even exists in the presence of noise.

From Transforms to Series To be completed.

Fourier Series As described in the last section (hopefully), we have seen that by restricting our data to a time interval [0, T] for period T, and extending the data to ( −∞ , ∞ ) , one generates a periodic function of infinite duration at the cost of losing data outside the fundamental range. This is not unphysical, as the data is typically taken over a finite period of time. Thus, any physical results in the analysis can be obtained be restricting the outcome to the given period. In typical problems one seeks a representation of the signal, valid for t ∈ [0, T ], as ∞ 1 f (t ) = A0 + ∑ [ A p cos(ω p t ) + B p sin( ω p t )], 2 p =1 where the angular frequency is given by 2π ω p = 2πf p = p. T Note that f(t) is periodic with period P=T: f (t + P) = f (t) . Most signals have an infinite period, so we already have made a first approximation. Thus, this has restricted the physical data to 0 < t < T . [That is, by limiting t, we are forced to using discrete frequencies, or the above sum as opposed to an integral in the Fourier Transform in the last section.]

One can extract the Fourier Coefficients ( A p , B p ) using the orthogonality of the trigonometric basis. Namely, such orthogonality for a set {φ n ( x), n = 1K ∞} over the interval [ a, b] with weight ρ ( x) is given by the condition b

∫φ

n

( x)φ m ( x) ρ ( x ) dx = 0, m ≠ n.

a

For the trigonometric functions, this is given by the relations  0, m ≠ n T  ∫0 cos(ω mt )cos(ω nt ) dt = T , m = n.  2  0, m ≠ n T  ∫0 sin(ω mt )sin(ω nt) dt = T , m = n.  2 T

∫ cos(ω

m

t ) sin( ω n t ) dt = 0.

0

These relations are proved in Appendix A. The coefficients are found to in Appendix B as Ap =

2T y (t ) cos(ω p t ) dt , p = 0,1, 2, K T ∫0

Bp =

2T y (t ) sin( ω p t ) dt , p = 1,2,K T ∫0

Discrete Series In the previous analysis we had restricted time to the interval [0,T], leading to a Fourier series with discrete frequencies and a periodic function of time. In reality, taking data can only be done at certain frequencies, thus eliminating high frequencies. Such a restriction on the frequency will lead to a discretization of the data in time. Another way to view this is that when recording data we sample at a finite number of time steps, limiting our ability to collect data with large oscillations. Thus, we not only have discrete frequencies but we also have discrete times. This can be captured in the following representation: M 1 A0 + ∑ [ Ap cos(ω p t ) + B p sin(ω pt )], 2 p =1 where the time, time step and trigonometric arguments are given by T 2πpn t n = n∆t , ∆t = , and ω p t = . N T

f (t n ) =

We need to determine M and the unknown coefficients. As for the Fourier series, we rely on an orthogonality principle, but this time replacing the integral by a sum. Again we will determine the unknowns in terms of the given samples of the function f(t). The orthogonality is provided in Appendix C. If we take N samples, we can the determine N N unknown coefficients A0 , A1 ,K, AN / 2 and B1 ,K , BN / 2 −1 . Thus, we can fix M = . Often 2 the coefficients B0 and BN / 2 are included for symmetry. Note that the corresponding sine function factors evaluate to zero, leaving these two coefficients arbitrary. Thus, we can take them to be zero. The full set of coefficients are found to be (in Appendix D) Ap =

2 N

Bp =

N

∑ y(t )cos( n

n =1

N

∑ y (t

2 N

n =1

n

2

), p = 1,K N

2

−1

) sin( 2πNpn ), p = 1, 2, K N − 1 2

A0 = AN =

2 π pn N

1 N

1 N

N

∑ y( t

n

),

n=1

N

∑ y( t

n

) cos(nπ ),

n=1

B0 = B N = 0 2

Matlab Implementation In this section we provide implementations of the discrete trigonometric transform in Matlab. The first implementation is a straightforward one which can be done in most programming languages. The second implementation makes use of matrix computations that can be performed in Matlab. Sums can be done with matrix multiplication, as describes in Appendix I. This eliminates the loops in the program. Direct Implementation without special use of matrix algebra % % DFT in a direct implementation % % Enter Data in y y=[7.6 7.4 8.2 9.2 10.2 11.5 12.4 13.4 13.7 11.8 10.1 ... 9.0 8.9 9.5 10.6 11.4 12.9 12.7 13.9 14.2 13.5 11.4 10.9 8.1]; % Get length of data vector or number of samples N=length(y); % Compute Fourier Coefficients for p=1:N/2+1 A(p)=0; B(p)=0;

for n=1:N A(p)=A(p)+2/N*y(n)*cos(2*pi*(p-1)*n/N)'; B(p)=B(p)+2/N*y(n)*sin(2*pi*(p-1)*n/N)'; end end A(N/2+1)=A(N/2+1)/2; % Reconstruct Signal - pmax is number of frequencies used in increasing order pmax=13; for n=1:N ynew(n)=A(1)/2; for p=2:pmax ynew(n)=ynew(n)+A(p)*cos(2*pi*(p-1)*n/N)+B(p)*sin(2*pi*(p-1)*n/N); end end % Plot Data plot(y,'o') % Plot reconstruction over data hold on plot(ynew,'r') hold off

Compact Implementation This implementation uses matrix products and is described in Appendix H. % % DFT in a compact implementation % % Enter Data in y y=[7.6 7.4 8.2 9.2 10.2 11.5 12.4 13.4 13.7 11.8 10.1 ... 9.0 8.9 9.5 10.6 11.4 12.9 12.7 13.9 14.2 13.5 11.4 10.9 8.1]; N=length(y); % Compute the matrices of trigonometric functions p=1:N/2+1; n=1:N; C=cos(2*pi*n'*(p-1)/N); S=sin(2*pi*n'*(p-1)/N); % Compute Fourier Coefficients A=2/N*y*C; B=2/N*y*S A(N/2+1)=A(N/2+1)/2; % Reconstruct Signal - pmax is number of frequencies used in increasing order pmax=13; ynew=A(1)/2+C(:,2:pmax)*A(2:pmax)'+S(:,2:pmax)*B(2:pmax)'; % Plot Data plot(y,'o') % Plot reconstruction over data

hold on plot(ynew,'r') hold off

Appendices A. Orthogonality of Trigonometric Basis 2π m , m = 1,2, K. ) T  0, m ≠ n T  ∫0 cos(ω mt )cos(ω nt ) dt = T , m = n.  2  0, m ≠ n T  ∫0 sin(ω mt )sin(ω nt) dt = T , m = n.  2

We want to prove the relations (for ω m =

T

∫ cos(ω

m

t) sin(ω n t ) dt = 0.

0

These are based upon the trigonometric identities in Appendix F. For ω m ≠ ±ω n T

∫ cos(ω mt )sin(ω nt ) dt = 0

1T sin((ω n + ω m )t ) + sin((ω n − ω m )t ) dt 2 ∫0 T

1  − cos((ω n + ω m )t ) − cos((ω n − ω m )t )  =  +  2 ωn +ωm ωn − ωm 0

=0 since (ωn ± ωm )T = 2π ( n ± m) and cos(2π ( n ± m)) = 1 for n and m integers. Similarly, the other integrals can be computed. T

∫ cos(ω mt )cos(ω nt ) dt = 0

1T cos((ω n + ω m )t ) + cos((ω n − ω m )t ) dt 2 ∫0 T

1  sin((ωn + ωm ) t ) sin((ωn − ωm ) t )  =  +  2  ωn +ωm ωn −ωm 0 =0 since sin(2π (n ± m )) = 0. Finally,

T

1T ∫0 sin(ωmt )sin(ωnt) dt = 2 ∫0 cos((ωn − ωm )t) − cos((ωn + ωm )t ) dt T

1  sin((ωn − ωm )t ) sin((ωn + ωm )t )  =  +  2  ωn −ωm ωn + ωm 0 =0 For ω m = ω n , we have instead T

T

0

0

2 ∫ cos(ω nt )cos(ω nt ) dt = ∫ cos (ω nt ) dt

1T [1 + cos(2ω n t ) ]dt 2 ∫0

=

T

1  sin(2ωn t)  = t + 2 2ω n  0 =

T . 2

T

T

0

0

2 ∫ sin(ωnt )sin(ωnt) dt = ∫ sin (ωnt ) dt

=

1T [1 − cos(2ω n t ) ]dt 2 ∫0 T

1  sin(2ω n t )  = t − 2 2ω n  0 =

T . 2

Here we have used the trigonometric identities cos 2 θ = sin 2 θ =

1 (1 + cos2 2θ ) and 2

1 (1 − cos 2 2θ ) , which are obtainable from the identities in Appendix F. 2

B. Derivation of Fourier Coefficients We now look at the Fourier series representation ∞ 1 f (t ) = A0 + ∑ [ Ap cos(ω p t ) + B p sin(ω pt )]. 2 p =1 We can use the orthogonality relations in the last Appendix in order to arrive at expressions for the Fourier coefficients. First, we integrate the series over one period.

T

∫ 0

T ∞ 1  f ( t ) dt = ∫  A0 + ∑ [ Ap cos(ω p t ) + B p sin(ω pt )]  dt p =1 0 2  T T ∞   1 T = A0 ∫ dt + ∑  Ap ∫ cos(ω pt ) dt + Bp ∫ sin(ω p t ) dt  2 0 p =1  0 0  T T ∞  sin(ω p t ) − cos(ω p t )  1  = A0 T + ∑  Ap + Bp 2 ωp 0 ωp  p =1  0  1 = A0T 2

This will give the expression A0 =

2T f ( t ) dt . T ∫0

Now multiply the series by cos(ω q t ) for some q = 1,2, K and integrate. T

∫ 0

T ∞ 1  f ()cos( t ω q t ) dt = ∫  A 0 + ∑ [ Ap cos(ω p t ) + B p sin(ω p t )] cos(ω q t ) dt p =1 0 2  T T T ∞   1 = A0 ∫ cos(ω q t ) dt + ∑  Ap ∫ cos(ω pt )cos(ω qt ) dt + Bp ∫ sin(ω pt )cos(ω q t ) dt  2 0 p =1  0 0  T

∞ 1 sin(ω q t )  T  = A0 + ∑  Ap δ p ,q + Bp ⋅ 0  2 ω q 0 p =1  2 

=

T A. 2 q

Here use was made of the known orthogonality relations from the last section. Also, the Kronecker delta was used, defined as 0, m ≠ n δ n, m =  . 1, m = n This renders the sum above a sum with all zero terms except for one, that for p = q. For example, if q = 3, we have T

t ω t ) dt = A ⋅ 0 + A ⋅ 0 + A ∫ f ()cos( 3

1

2

3

0

Thus, one can solve for the Aq 's to obtain Aq = Similarly, one has

2 T

T T + A4 ⋅ 0 + L = A3 . 2 2

T

t ω t) dt , ∫ f ()cos( q

0

q = 1,2,K.

T

∫ 0

T ∞ 1  f ()sin( t ω q t ) dt = ∫  A0 + ∑ [ Ap cos(ω p t ) + B p sin(ω p t )] sin(ω q t ) dt p =1 02  T T ∞   1 T = A0 ∫ sin(ω q t )dt + ∑  Ap ∫ cos(ω p t )sin(ω q t ) dt + Bp ∫ sin(ωp t )sin(ωq t ) dt  2 0 p =1  0 0  T

∞ 1 − cos(ω q t ) T   = A0 + ∑  Ap ⋅ 0 + Bp δ p ,q  2 ωq 2  p =1  0

=

T B. 2 q

2T So, Bq = ∫ f ()sin( t ωq t ) dt . T0

C. Discrete Orthogonality The derivation of the discrete Fourier coefficients can be done using the discrete orthogonality of the discrete trigonometric basis similar to the derivation of the above Fourier coefficients for the Fourier series. We first prove the following  2π nk   0, k = 1,K , N −1 = k = 0, N N  N , n =1 N 2π nk  sin  ∑  = 0, k = 0,K , N  N  n =1 N

∑ cos 

This can be done more easily using the exponential form,  2π nk  + i N sin  2π nk  = N e 2 π ink / N , cos ∑  N  ∑  N  ∑   n =1   n =1 n =1 iθ by using Euler's formula, e = cos θ + i sin θ for each term in the sum. N

The exponential sum is the sum of a geometric progression, which can be summed as done in Appendix G. Thus, we have N

∑e n =1

2π ink / N

= ∑ (e N

n =1

2π ik / N

)

n

= e2 π ik / N + ( e2 π ik / N ) + L + ( e2 π ik / N ) 2

N e 2π ik / N 1 − ( e 2π ik / N )    = 2π ik / N 1− e 2π ik / N 1 − e 2π ik  e = . 1 − e 2π ik / N

N

As long as k ≠ 0, N the numerator is 0. In the special cases that k = 0, N , we have e 2π ink / N = 1. So, N

N

n =1

n =1

∑ e 2π ink / N = ∑ 1 = N . Therefore, N

 2π nk  N  2π nk   0, k = 1,K , N −1 and the result follows.  + i ∑ sin  = k = 0, N N  n =1  N   N ,

∑ cos  n =1

We can use this to establish orthogonality relations of the following type for N p , q = 0,1,K , : 2  2π pn   2π qn  1  N  2π ( p − q ) n   2π ( p + q )n   cos cos ∑     =  ∑ cos   + cos   . N N  N   N  2  n =1     n =1 N

Splitting the above sum into two sums and then evaluating the separate sums from earlier in this section, N  2π ( p − q ) n   0, p ≠ q cos  ∑  = N, p = q, N    n =1 N

 2π ( p + q )n   0,  = N , N  

∑ cos  n= 1

p+ q≠ N p+q=N

we obtain  N /2,  2π pn   2π qn   cos  ∑  cos  N  =  N ,  N     n =1  0, N

p = q≠ N/2 p = q = N /2. otherwise

Similarly, we find N

∑ sin( n =1

2π pn 2π qn 1 N 2π ( p − q ) n 2π ( p + q) n  )cos( ) =  ∑ sin( ) + sin( )  = 0. N N 2  n =1 N N 

and N

∑ sin( n =1

2π pn 2π qn 1N 2π ( p − q ) n 2π ( p + q ) n  )sin( ) = ∑ cos( ) − cos( ) N N 2  n =1 N N   N /2, =  0,

p = q≠ N/2 . otherwise

D. Discrete Transform The derivation of the coefficients for the DFT is now easily obtained using the discrete orthogonality in the last section. We start with the expansion f (t n ) =

N /2 1  2π pn   2π pn  A0 + ∑ [ Ap cos  + B p sin    ]. 2  N   N  p =1

We first sum over n: N

∑ n =1

N /2 1 2π pn  2π pn   A + [ Ap cos  + Bp sin   ∑ ∑ 0  ]  N   N  n =1  2 p =1 N N /2 N 1  N  2π pn   2π pn   = A0 ∑1 + ∑  Ap ∑ cos  + Bp ∑ sin    2 n =1  N   N  p =1  n =1 n =1

f ( tn ) =

N

=

N /2 1 A0 N + ∑ { Ap ⋅ 0 + Bp ⋅ 0} 2 p =1

=

1 A N. 2 0

2π qn  Now, we can multiply both sides of the expansion by cos   and sum over n.  N  N N /2  2π qn  N  1  2π pn   2π pn    2π qn  f (t n )cos  = ∑  A0 + ∑ [ Ap cos  + Bp sin  ]  cos  ∑      N  n =1  2  N   N    N  n =1 p =1 =

N 1  2π qn  A0 ∑ cos   2 n =1  N 

N /2 N  N  2π pn   2π qn   2π pn   2π qn   + ∑  Ap ∑ cos  cos + B    p ∑ sin   cos    N   N   N   N  p =1  n =1 n =1  N /2  N   ∑  Ap 2 δ p , q + B p ⋅ 0  q ≠ N / 2   p =1  = N / 2  { A Nδ q= N/2 p p , N / 2 + Bp ⋅ 0}  ∑ p =1 1  A N, q ≠ N / 2 = 2 q .  AN / 2 N , q = N / 2 So, we have found that 2 N 2π qn  N Aq = ∑ f ( tn )cos  , q≠ ,  N n =1 2  N  1 N 2π n( N / 2 )  1 N AN / 2 = ∑ f ( t n )cos   = N ∑ f ( tn )cos (π n ) . N n =1 N   n =1

Similarly, N

∑ f (t n =1

n

N /2  2π qn  N  1  2π pn   2π pn    2π qn  )sin  = ∑  A0 + ∑ [ Ap cos  + B p sin     ] sin    N  n =1  2 p =1  N   N    N 

=

N 1 N  2π qn  N / 2  N  2π pn   2π qn   2π pn   2π qn   A0 ∑ sin  + ∑  Ap ∑ cos  sin  + B p ∑ sin      sin   2 n =1  N  p =1  n =1  N   N   N   N  n =1

N /2 N = ∑  A p ⋅ 0 + B p δ p, q  2  p=1  1 = Bq N . 2

Finally, we have 2 N 2π qn  N Bq = ∑ f (t n )sin  , q = 1, K, − 1.  N n =1 2  N 

E. The Discrete Exponential Transform The derivation of the coefficients for the DFT was obtained using the discrete orthogonality in the last section. However, this is not the form used in Matlab for spectral analysis. Matlab allows for the computation of the Fast Fourier Transform (FFT) and its description in the help section does not involve sines and cosines. Namely, Matlab defines the transform and inverse transform as "For length N input vector x, the DFT is a length N vector X, with elements N X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1