numerical integration in Fourier space with numpy.fft

Question

I want to integrate a function with the numerical integration in the Fourier space.

The following code shows a working example:

import numpy as np
from pylab import *
from numpy.fft import fft, ifft, fftshift, ifftshift

N = 2**16
x = np.linspace(- np.pi , np.pi,N)
y = np.exp(-x**2)                          # function f(x)
ys = np.exp(-x**2) * (-2*x)                # derivative f'(x)
T = x[-1] - x[0] # the whole range
w = (np.arange(N) - N /2.) / T + 0.00000001 # slightly shifted
# integration
fourier =  ifft(ifftshift(1./ ( 2 * np.pi * 1j * w) * fftshift(fft(ys)) )  ) 
# differentiation 
fourier2 =  ifft(ifftshift(( 2 * np.pi * 1j * w) * fftshift(fft(fourier)) )  ) 

You might notice the + 0.00000001 in the definition of the frequencies w. I need it, because otherwise I would generate a ZeroDivisionError or a numpy warning. This is a workaround, which seems to be okay for the example above, but it fails for a more complex problem I am having. A co-worker told me, that I could simply avoid it, if I got the fft for shifted frequency values (np.arange(N) - N /2. + 1./2) / T. How to do that in numpy? Is there a way to specify the output grid of the numpy fft?

Thanks!

Answer 1

The problem is that w contains 0 (as it should), and you divide by w. The 0 in w is the "DC" frequency; it corresponds to the constant term of the Fourier series.

If you integrate a function with a nonzero DC component with coefficient A0, the resulting function includes a term of the form A0*t, which is not in the space of periodic functions to which this Fourier technique applies. So you must assume that the DC component of the input is 0. In this case, when you divide by w, you'll get (0+0j)/0, which is (nan+nanj). If the DC component of the input is not zero, you'll get (inf+nanj). Either way, the solution is to simply ignore whatever you get, and set the DC Fourier coefficient to 0 before inverting with ifft.

There are several ways to implement this. One way is to change these lines:

w = (np.arange(N) - N /2.) / T + 0.00000001 # slightly shifted
# integration
fourier =  ifft(ifftshift(1./ ( 2 * np.pi * 1j * w) * fftshift(fft(ys)) )  ) 

to this (I added a couple intermediate variables):

w = (np.arange(N) - N /2.) / T

# integration
Fys = fft(ys)
with np.errstate(divide="ignore", invalid="ignore"):
    modFys = ifftshift(1./ (2 * np.pi * 1j * w) * fftshift(Fys))

# modFys[0] will hold the result of dividing the DC component of y by 0, so it
# will be nan or inf.  Setting modFys[0] to 0 amounts to choosing a specific
# constant of integration.
modFys[0] = 0

fourier = ifft(modFys).real

I also took the real part of the result of ifft. Theoretically, the imaginary parts should all be 0; in practice they will be very small but nonzero because of normal inexact floating point arithmetic.

By the way, if you'd rather not implement your own version of this technique, you could use scipy.fftpack.diff.