Using Fourier Coefficients to approximate a time series

Question

I'm trying to understand approximating periodic functions with Fourier Transforms. I'm using R, but I suppose the question is language agnostic.

Say I have this discrete periodic series. It is some process that is measured each second and cycles every 5 seconds:

y = rep(c(1,2,1,-1,0), 5)

The FFT will give me the Fourier Coefficients to approximate the observed series as a sum of complex exponentials:

$$ X(t) \approx \frac{1}{N} \sum_{k} Re(a_k)cos(2\pi\frac{k}{N}t) + Im(a_k)sin(2\pi\frac{k}{N}t) $$

where $N$ is the sample size, $t$ it sime, $a_k$ is the k-th Fourier coefficient, $Re(a_k)$ is its real component, and $Im(a_k)$ is the imaginary component.

So I do the FFT:

# do fft
z = fft(y)
# it's symmetric just keep first half (ignoring first value, which is sum of the series)
z1 = z[2:(length(z)/2 + 1)]

and checking the coefficients I see action at the fifth and tenth coefficients:

> round(z1)
[1]  0+ 0i  0+ 0i  0+ 0i  0+ 0i  8-15i  0+ 0i  0+ 0i  0+ 0i  0+ 0i -3+ 4i  0+ 0i  0+ 0i

so does that mean I can approximate the series with:

$$ 8cos(2\pi\frac{5}{N}t) -15sin(2\pi\frac{5}{N}t) - 3cos(2\pi\frac{10}{N}t) + 4sin(2\pi\frac{10}{N}t) $$

because when I do that, I get something that sort of looks like the series, but not really:

t = 1:length(y)
f1 = 1/N * (Re(z[5]) * cos(2*pi*(5/N)*t) - Im(z[5]) * sin(2*pi*(5/N)*t))
f2 = 1/N * (Re(z[10]) * cos(2*pi*(10/N)*t) - Im(z[10]) * sin(2*pi*(10/N)*t))
f = f1 + f2
plot(f, type = 'l')

I'm sure I'm doing something wrong but I've read a bunch of tutorials and still can't figure it out. I'm not sure if I'm correctly accounting for period, frequency, phase...

Answer 1

The equation you want to use is that the k-th point in your original series $x_k$ is reproduced as:

\begin{equation} x_k = \frac{1}{N}\sum_{n=0}^{N-1}X_n e^{i\frac{2\pi}{N}kn} \end{equation}

Where $X_n$ is the n-th Fourier coefficient.

This can be implemented in Python as:

from scipy.fft import fft
import numpy as np

x = [1,2,1,-1,0]*5

# compute fourier coefficients
coef = fft(x)

n = len(x)
# implement the inverse transform equation
z = [np.sum(coef * np.exp(2j * np.pi * k * np.arange(n)/n)) / n for k in range(n)]

# disregard small imaginary component and round to integer form
z = [np.int(np.real(z_i)) for z_i in z]

Doing this we find that $z$ is equal to the original series $x$.

I think your error may be down to throwing out values from $X_n$, even though there are repeated values you need all of them to reconstruct your original time series.