3.5 The Fourier Transform

Although presentation of Fourier coefficients via sines and cosines has intuitive appeal, we can present the same ideas in a more compact manner using complex exponentials. This approach will allow for a simple presentation of the fast Fourier transform (FFT) algorithm in the following section.

The Fourier transform of a time series \(y_t\) for frequency \(p\) cycles per \(n\) observations can be written as \[ z_p = \sum_{t=0}^{n-1} y_t\exp(-2\pi i \,p\,t / n). \] for \(p = 0, \dots, n-1\). Note that the above expression differs slightly from what we have presented in the previous sections but is consistent with how R computes the Fourier transform.

From the expression above, it’s clear that \(z_p\) is a complex number. In this case, the real part corresponds to the cosine variation and the imaginary part corresponds to the sine variation. We can see this because of Euler’s relation \[ \exp(i\theta) = \cos\theta + i\sin\theta \] Expanding the definition of \(z_p\) above gives us \[\begin{eqnarray*} z_p & = & \sum_{t=0}^{n-1} y_t\exp(-2\pi i p t/ n)\\ & = & \sum_{t=0}^{n-1} y_t(\cos(-2\pi p t / n) + i\sin(-2\pi p t/n))\\ & = & \sum_{t=0}^{n-1} y_t\cos(-2\pi p t/ n) + i\sum_{t=0}^{n-1}\,y_t\sin(-2\pi p t/n)\\ & = & \alpha_p + i\beta_p \end{eqnarray*}\]

So \(z_p\) for \(p=0,\dots,n-1\) represents the (complex) Fourier coefficient associated with frequency \(p\).

The inverse Fourier transform is then \[ y_t = \frac{1}{n}\sum_{p=0}^{n-1} z_p\exp(2\pi i p t/n). \]

Note the change of sign in the exponential when inverting the Fourier transform.

For real-valued time series, the full set of Fourier coefficients is redundent because \(z_p = \bar{z}_{n-p}\). But the FFT algorithm in R (for example) computes all of the coefficients for any series.

3.5.1 Example

We can take a look at how R computes the Fourier transform via the fft() function (the details of which will be discussed in the next section).

Here is a simple time series.

y <- rnorm(10)
y
 [1] -0.22407225 -1.20389294 -0.07886933  0.80119409  1.47853611 -0.19530257
 [7]  0.32241714 -0.95226477  0.14767474 -0.55742245

And here is its Fourier transform.

z <- fft(y)
z
 [1] -0.4620022+0.0000000i -2.8427601-1.7517461i -0.3405720+2.8781810i
 [4]  0.8941484+0.4127935i -0.4768639-0.8235592i  3.7533751-0.0000000i
 [7] -0.4768639+0.8235592i  0.8941484-0.4127935i -0.3405720-2.8781810i
[10] -2.8427601+1.7517461i

The first coefficient is always the sum of the data.

sum(y)
[1] -0.4620022

The second coefficient has a real and imaginary component. For the real component, we can see that it is equal to

sum(y * cos(2 * pi * 1 * 0:9 / 10))
[1] -2.84276

The imaginary component is then equal to

sum(y * sin(-2 * pi * 1 * 0:9 / 10))
[1] -1.751746

We can even use complex exponentials directly by computing

z1 <- sum(y * exp(-2 * pi * 1i * 1 * 0:9 / 10))
z1
[1] -2.84276-1.751746i

Note that the imaginary number \(i\) is denoted by 1i in R and the real and imaginary components of a complex number can be obtained via the Re() and Im() functions.

Re(z1)
[1] -2.84276
Im(z1)
[1] -1.751746

If we look at the sequence of complex Fourier coefficients carefully we notice that there is a symmetry to it.

 [1] -0.4620022+0.0000000i -2.8427601-1.7517461i -0.3405720+2.8781810i
 [4]  0.8941484+0.4127935i -0.4768639-0.8235592i  3.7533751-0.0000000i
 [7] -0.4768639+0.8235592i  0.8941484-0.4127935i -0.3405720-2.8781810i
[10] -2.8427601+1.7517461i

For example look at

z[2 + 1]
[1] -0.340572+2.878181i
z[10 - 2 + 1]
[1] -0.340572-2.878181i

The latter number is the complex conjugate of the first number. (Note that the indices are written with + 1 because the mathematics wants to index starting at \(0\) but R wants to start at \(1\).)

When doing spectral analysis, we care about the size of the complex Fourier coefficient,

\[ R=\sqrt{\alpha_p^2 + \beta_p^2} \]

which can be computed with the Mod() function.

Mod(z)
 [1] 0.4620022 3.3391465 2.8982607 0.9848349 0.9516559 3.7533751 0.9516559
 [8] 0.9848349 2.8982607 3.3391465

Note that because we square the coefficients here, their signs do not matter and so the values beyond index \(n/2 + 1\) are a mirror image of the values before them.

Finally, what if we have our Fourier coefficients and want to retrieve our original series? We can compute the inverse Fourier transform via fft(inverse = TRUE), so

Re(fft(z, inverse = TRUE) / 10)
 [1] -0.22407225 -1.20389294 -0.07886933  0.80119409  1.47853611 -0.19530257
 [7]  0.32241714 -0.95226477  0.14767474 -0.55742245

Note that we had to divide by \(n\) and take the real component explicitly because the fft() function always returns complex numbers.

We can verify directly, that we can retrieve \(y_0\) if we compute the sum

mean(z * exp(2 * pi * 1i * 0:19 * 0 / 10))
[1] -0.2240722-0i

Similarly, we can retrieve \(y_1\) by computing

mean(z * exp(2 * pi * 1i * 0:19 * 1 / 10))
[1] -1.203893+0i

Direct computation like that shown above is illustrative but not efficient. In the next section we will describe the fast Fourier transform which provides a dramatically faster algorithm for computing the Fourier coefficients.