How does numpy.fft.fft work?

Question

I am currently trying to understand the fft-function from numpy. For that I tested the following assumption:
I have two functions, f(x) = x^2 and g(x) = f'(x) = 2*x. According to the fourier transformation laws and wolfram alpha it should be that G(w) = 2pi*i*F(w) (prefactors can vary, but there should only be a constant factor). When implementing that in python, I write

import numpy as np
def x2(x):
    return x*x
def nx(x):
    return 2*x

a = np.linspace(-3, 3, 16)
a1 = x2(a)
a2 = nx(a)

b1 = np.fft.fft(a1)
b2 = np.fft.fft(a2)

c = b1/b2

Now I am expecting a nearly constant value for c, but I get

array([  1.02081592e+16+0.j        ,   1.32769987e-16-1.0054679j ,
         4.90653893e-17-0.48284271j,  -1.28214041e-16-0.29932115j,
        -1.21430643e-16-0.2j       ,   5.63664751e-16-0.13363573j,
        -5.92271642e-17-0.08284271j,  -4.21346622e-16-0.03978247j,
        -5.55111512e-16-0.j        ,  -5.04781597e-16+0.03978247j,
        -6.29288619e-17+0.08284271j,   8.39500693e-16+0.13363573j,
        -1.21430643e-16+0.2j       ,  -0.00000000e+00+0.29932115j,
        -0.00000000e+00+0.48284271j,   1.32769987e-16+1.0054679j ])

Where is my mistake, and what can I do to use the fft as intended?

Answer 1

The properties you give apply to the Continuous Fourier transform (CFT). What is computed by the FFT is the Discrete Fourier transform (DFT), which is related to the CFT but is not exactly equivalent.

It's true that the DFT is proportional to the CFT under certain conditions: namely with sufficient sampling of a function that is zero outside the sample limits (see e.g. Appendix E of this book).

Neither condition holds for the functions you propose above, so the DFT is not proportional to the CFT and your numerical results reflect that.

---

Here's some code that confirms via the FFT the relationship you're interested in, using an appropriately sampled band-limited function:

import numpy as np

def f(x):
    return np.exp(-x ** 2)
def fprime(x):
    return -2 * x * f(x)

a = np.linspace(-10, 10, 100)
a1 = f(a)
a2 = fprime(a)

b1 = np.fft.fft(a1)
b2 = np.fft.fft(a2)
omega = 2 * np.pi * np.fft.fftfreq(len(a), a[1] - a[0])

np.allclose(b1 * 1j * omega, b2)
# True