FFT with actual frequencies

real-fft.py

# This gist shows how to use numpy's FFT functions to use the discrete FFT
# of a uniformly-spaced data set to get what corresponds to the actual Fourier
# transform as you'd write it down mathematically.
# In particular, the functio uniform_dft returns the corresponding frequencies
# in terms of the coordinate units.

def uniform_dft(interval_length, data):
    """Discrete Fourier Transform of real-valued data, interpreted for a uniform series
    over an interval of a given length, including starting and end point.
    The function returns the frequencies in units of the coordinate.
    """
    X = numpy.fft.rfft(data)
    n = len(data)
    # The input data is assumed to cover the entire time interval, i.e., including start
    # and end point. The data produced from RFFT however assumes that the end point is
    # excluded. Hence, stretch the interval_length such that cutting off the end point
    # results in the interval of length interval_length.
    interval_length *= n / float(n - 1)
    # Note that the following definition of the frequencies slightly differs from the
    # output of np.fft.freqs which is
    #
    #     freqs = [i / interval_length / n for i in range(n//2 + 1)].
    freqs = numpy.arange(n // 2 + 1) / interval_length
    # Also note that the angular frequency is  omega = 2*pi*freqs.
    #
    # With RFFT, the amplitudes need to be scaled by a factor of 2.
    X /= n
    X[1:-1] *= 2
    if n % 2 != 0:
        X[-1] *= 2
    assert len(freqs) == len(X)
    return freqs, X


def uniform_idft(X, n):
    # Equivalent:
    # n = 2 * len(freqs) - 1
    # data = numpy.zeros(n)
    # t = numpy.linspace(0.0, len_interval, n)
    # for x, freq in zip(X, freqs):
    #     alpha = x * numpy.exp(1j * 2 * numpy.pi * freq * t)
    #     data += alpha.real

    # transform back
    X[1:-1] /= 2
    if n % 2 != 0:
        X[-1] /= 2
    X *= n
    data = numpy.fft.irfft(X, n)
    return data