Terminology
Nothing in this question is specific to the fast Fourier transform (FFT).
The FFT is a particular algorithm for computing the discrete Fourier transform (DFT), so I'm going to say "DFT" instead of "FFT".
In this answer, m indicates a discrete time index and k indicates a discrete frequency index.
What is the DFT?
There are several problems here, all of which are of mathematical nature coming from misunderstanding how the DFT works.
Taken from the numpy.fft module docstring, numpy defines the discrete Fourier transform as
A_k = \sum_{m=0}^{n-1} a_m \exp[-2 \pi i (m k / n)]
That's LaTeX notation saying that the discrete Fourier transform is a linear combination of complex exponentials exp[2 pi i m k / n] where n is the total number of points and m is the discrete time index.
In your notation this would be exp[2 pi i m k / N] because you're using N to mean the total number of points.
Use exp, not sine
Note that the DFT uses exponentials; these are not sine functions.
If you want to build up the time domain signal from the discrete Fourier transform coefficients, you have to use the same functions that the DFT itself used!
Therefore, you could try this:
FFF = zeros(N)
for i in range(300):
FFF += fftcoeffs[i]*np.exp(2*pi*i*freqs[i]*t)
plt.plot(t, FFF)
plt.show()
However, that, too, will fail in a way which might puzzle you.
Aliasing
The final puzzle piece has to do with an effect called aliasing.
Suppose you take the DFT of a signal exp[2 pi i (N + p) m / N].
If you compute it, you find that all of the A_k are zero except for A_p.
In fact, you get the same thing that you would get if you took the DFT of exp[2 pi i p m / N].
You can see that any complex exponential with frequency larger than N shows up as if it were a lower frequency.
In particular, any complex exponential with frequency q + b N where b is any integer looks like it has frequency q.
Now suppose we have a time domain signal cos(2 pi p m / N).
That's equal to
(1/2)[ (exp(2 pi i p m / N) + exp(-2 pi i p m / N) ].
That negative frequency has interesting consequences for the DFT.
The frequency -p can be written as (N-p) + N.
This has the form q + b N with q = N - p and b=1.
So, that negative frequency -p winds up looking like N - p!
The numpy function fftfreq knows this.
Take at look at the output of fftfreq and you'll see that it starts at zero, runs up to half the sampling frequency (called the Nyquist frequency), and then goes negative!
This is to help you deal with the aliasing effect we just described.
The upshot of all this, is that if you want to approximate a function by summing the lowest frequency Fourier components, you do not want to take the lowest few elements from fftfreq.
Instead, you want to sort fftfreq by absolute value and then sum up complex exponentials with those frequencies.
Also take a look at np.fft.hfft.
This function is designed to handle real valued functions and the aliasing issues associated with them.
Code
Since this is a fairly difficult issue to discuss verbally, here's a script which does exactly what you want.
Note that I put comments after the block of code those comments describe.
Make sure you have matplotlib installed (in your virtual environment... you use virtual environments, right?).
If you have questions leave a comment.
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
pi = np.pi
def square_function(N, square_width):
"""Generate a square signal.
Args:
N (int): Total number of points in the signal.
square_width (int): Number of "high" points.
Returns (ndarray):
A square signal which looks like this:
|____________________
|<-- square_width -->
| ______________
|
|^ ^ ^
index |0 square_width N-1
In other words, the output has [0:N]=1 and [N:]=0.
"""
signal = np.zeros(N)
signal[0:square_width] = 1
return signal
def check_num_coefficients_ok(N, num_coefficients):
"""Make sure we're not trying to add more coefficients than we have."""
limit = None
if N % 2 == 0 and num_coefficients > N // 2:
limit = N/2
elif N % 2 == 1 and num_coefficients > (N - 1)/2:
limit = (N - 1)/2
if limit is not None:
raise ValueError(
"num_coefficients is {} but should not be larger than {}".format(
num_coefficients, limit))
def test(N, square_width, num_coefficients):
"""Test partial (i.e. filtered) Fourier reconstruction of a square signal.
Args:
N (int): Number of time (and frequency) points. We support both even
and odd N.
square_width (int): Number of "high" points in the time domain signal.
This number must be less than or equal to N.
num_coefficients (int): Number of frequencies, in addition to the dc
term, to use in Fourier reconstruction. This is the number of
positive frequencies _and_ the number of negative frequencies.
Therefore, if N is odd, this number cannot be larger than
(N - 1)/2, and if N is even this number cannot be larger than
N/2.
"""
if square_width > N:
raise ValueError("square_width cannot be larger than N")
check_num_coefficients_ok(N, num_coefficients)
time_axis = np.linspace(0, N-1, N)
signal = square_function(N, square_width)
ft = np.fft.fft(signal)
reconstructed_signal = np.zeros(N)
reconstructed_signal += ft[0] * np.ones(N)
# Adding the dc term explicitly makes the looping easier in the next step.
for k in range(num_coefficients):
k += 1 # Bump by one since we already took care of the dc term.
if k == N-k:
reconstructed_signal += ft[k] * np.exp(
1.0j*2 * pi * (k) * time_axis / N)
# This catches the case where N is even and ensures we don't double-
# count the frequency k=N/2.
else:
reconstructed_signal += ft[k] * np.exp(
1.0j*2 * pi * (k) * time_axis / N)
reconstructed_signal += ft[N-k] * np.exp(
1.0j*2 * pi * (N-k) * time_axis / N)
# In this case we're just adding a frequency component and it's
# "partner" at minus the frequency
reconstructed_signal = reconstructed_signal / N
# Normalize by the number of points in the signal. numpy's discete Fourier
# transform convention puts the (1/N) normalization factor in the inverse
# transform, so we have to do it here.
plt.plot(time_axis, signal,
'b.', markersize=20,
label='original')
plt.plot(time_axis, reconstructed_signal.real,
'r-', linewidth=3,
label='reconstructed')
# The imaginary part is zero anyway. We take the real part to
# avoid matplotlib warnings.
plt.grid()
plt.legend(loc='upper right')
