2
$\begingroup$

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)

enter image description here

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')

enter image description here

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...

$\endgroup$
2
  • 1
    $\begingroup$ I'm not really sure cross validated is the best place for this question. There is a dedicated digital signal processing stack exchange: dsp.stackexchange. Or you could try the mathematics stack exchange: math.stackexchange. $\endgroup$
    – user27119
    Sep 21, 2021 at 10:06
  • 1
    $\begingroup$ ok i'll do that, cheers $\endgroup$
    – invictus
    Sep 21, 2021 at 17:54

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ OK. This is cool. But is it just the same as doing an IFFT on the FFT? I was hoping to end up with the smoothed approximation function, not the original time series. $\endgroup$
    – invictus
    Sep 21, 2021 at 17:44
  • $\begingroup$ @invictus: I didn't look at the detail of Adam's solution ( thanks Adam. I plan to try to understand what you did at some point ) but a periodic series always series has an exact fourier series representation so nothing is being approximated. Therefore one would expect to get back the original series if things are done correctly like Adam did them. See Silov's book on fourier series for a nice discussion and examples. $\endgroup$
    – mlofton
    Sep 21, 2021 at 18:12
  • $\begingroup$ The equation you used in your question is actually just the inverse Fourier transform. By approximate all you mean is that you remove some of the noisy frequencies. There are two reasons your outcome looks weird. 1) You didn’t get the original time series back because you removed half the coefficients you need to recover it. 2) The smoothing/approximating works by applying the inverse transform but with the smallest frequencies/coefficients set equal to zero. But you’ll notice on your coefficients, they are basically all zero since your time series is so simple to begin with. $\endgroup$
    – Adam Kells
    Sep 21, 2021 at 20:19
  • $\begingroup$ Try adding some noise to the original time series. Then go through your steps and you’ll be able to approximate it by setting the small non-zero coefficients to zero. $\endgroup$
    – Adam Kells
    Sep 21, 2021 at 20:20
  • $\begingroup$ Adam: Is that not correct that the reason is because a periodic series can be represented exactly by its fourier series coefficients. I'm no DSP whiz so definitely opened to corrections.. Thanks. $\endgroup$
    – mlofton
    Sep 21, 2021 at 21:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.