Manually inverting FFT using Numpy

Question

I have a little script for calculating the Fourier Transform of a square wave which works well and returns the square wave correctly when I invert the fft using numpy.fft.ifft(). However, I am unable to invert the transform by manually adding up harmonics after multiplying them by their respective coefficients that I obtain from numpy.fft.fft() Below is my script and I believe you'll see my intention.

from numpy import zeros, concatenate, sin, pi, linspace
from numpy.fft import fft, fftfreq, ifft
import numpy as np
import matplotlib.pyplot as plt

N = 1024 # samples
T = 1 # period
dt = T/N # sampling period
fs = 1/dt # sampling frequency
t = linspace(0, T, N) # time points
functime = .... # square wave

funcfft = fft(functime) # fft
fftcoeffs = np.abs(funcfft)/N # coefficients, divide by N to get actual coeff.s(I believe?)
freqs = fftfreq(N, dt) # frequencies

plt.plot(freqs, fftcoeffs) # gives me reasonable output
plt.show()

FF = ifft(funcfft)
plt.plot(t, FF) # plots exactly the same function as functime defined above
plt.show()

All is well until this far. Now my question is, shouldn't I converge to the original function if I run the below script after the above script?:

FFF = zeros(N)
for i in range(300):
    FFF += fftcoeffs[i]*sin(2*pi*freqs[i]*t)
plt.plot(t, FFF)
plt.show()

Assume that range(300) is good enough for convergence. Now when I do this, FFF is different than my original function. I thought that if I multiplied harmonics of respective frequencies by their corresponding coefficients, which I think are stored in fftcoeffs, I would then converge to the original function. What am I doing wrong?

Update: According to DanielSank's suggestions, I have updated my for loop as below, which unfortunately did not give me the desired results:

freqs2 = np.abs(freqs)
freqs2 = np.sort(freqs2)
for k in range(300):
    FFF += fftcoeffs[k]*exp(2j*pi*freqs2[k]*t/N)

I am not sure if I am doing the "sort fftfreq by absolute value" part right here.

Answer 1

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