Is there a way to numerically integrate in Fourier space with scipy.fft?

Question

I am interested in integrating in Fourier space after using scipy to take an fft of some data. I have been following along with this stack exchange post numerical integration in Fourier space with numpy.fft but it does not properly integrate a few test cases I have been working with. I have added a few lines to address this issue but still am not recovering the correct integrals. Below is the code I have been using to integrate my test cases. At the top of the code are the 3 test cases I have been using.

import numpy as np
import scipy.special as sp
from scipy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt

#set number of points in array
Ns = 2**16 
#create array in space
x = np.linspace(-np.pi, np.pi, Ns)

#test case 1 from stack exchange post
# y = np.exp(-x**2)                          # function f(x)
# ys = np.exp(-x**2) * (-2 *x)               # derivative f'(x)

#test case 2 
# y = np.exp(-x**2) * x - 1/2 *np.sqrt(np.pi)*sp.erf(x)
# ys =  np.exp(-x**2) * -2*x**2

#test case 3
y = np.sin(x**2) + (1/4)*np.exp(x)
ys = 1/4*(np.exp(x) + 8*x*np.cos(x**2))


#find spacing in space array
ss = x[1]-x[0]

#definte fft integration function
def fft_int(N,s,dydt):
    #create frequency array
    f = fftfreq(N,s)

    # integration step ignoring divide by 0 errors
    Fys = fft(dydt)
    with np.errstate(divide="ignore", invalid="ignore"):
        modFys = Fys / (2*np.pi*1j*f) 

    #set DC term to 0, was a nan since we divided by 0
    modFys[0] = 0

    #take inverse fft and subtract by integration constant
    fourier = ifft(modFys)
    fourier = fourier-fourier[0]
    
    #tilt correction if function doesn't approach 0 at its ends
    tilt = np.sum(dydt)*s*(np.arange(0,N)/(N-1) - 1/2)
    fourier = fourier + tilt
    
    return fourier

Test case 1 was from the stack exchange post from above. If you copy paste the code from the top answer and plot you'll get something like this:

with the solid blue line being the fft integration method and the dashed orange as the analytic solution. I account for this offset with the following line of code:

fourier = fourier-fourier[0]

since I don't believe the code was setting the constant of integration.

Next for test case 2 I get a plot like this:

again with the solid blue line being the fft integration method and the dashed orange as the analytic solution. I account for this tilt in the solution using the following lines of code

tilt = np.sum(dydt)*s*(np.arange(0,N)/(N-1) - 1/2)
fourier = fourier + tilt

Finally we arrive at test case 3. Which results in the following plot:

again with the solid blue line being the fft integration method and the dashed orange as the analytic solution. This is where I'm stuck, this offset has appeared again and I'm not sure why.

TLDR: How do I correctly integrate a function in fourier space using scipy.fft?

Answer 1

The tilt component makes no sense. It fixes one function, but it's not a generic solution of the problem.

The problem is that the FFT induces periodicity in the signal, meaning you compute the integral of a different function. Multiplying the FFT of the signal by 1/(2*np.pi*1j*f) is equivalent to a circular convolution of the signal with ifft(1/(2*np.pi*1j*f)). "Circular" is the key here. This is just a boundary problem.

Padding the function with zeros is one way to attempt to fix this:

import numpy as np
import scipy.special as sp
from scipy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt

def fft_int(s, dydt, N=0):
    dydt_padded = np.pad(dydt, (0, N))
    f = fftfreq(dydt_padded.shape[0], s)
    F = fft(dydt_padded)
    with np.errstate(divide="ignore", invalid="ignore"):
        F = F / (2*np.pi*1j*f) 
    F[0] = 0
    y_padded = np.real(ifft(F))
    y = y_padded[0:dydt.shape[0]]
    return y - np.mean(y)

N = 2**16
x = np.linspace(-np.pi, np.pi, N)
s = x[1] - x[0]

# Test case 3
y = np.sin(x**2) + (1/4)*np.exp(x)
dy = 1/4*(np.exp(x) + 8*x*np.cos(x**2))
plt.plot(y - np.mean(y))
plt.plot(fft_int(s, dy))
plt.plot(fft_int(s, dy, N))
plt.plot(fft_int(s, dy, 10*N))
plt.show()

(Blue is expected output, computed solution without padding is orange, and with increasing amount of padding, green and red.)

Here I've solved the "offset" problem by plotting all functions with their mean removed. Setting the DC component to 0 is equal to subtracting the mean. But after cropping off the padding the mean changes, so fft_int subtracts the mean again after cropping.

Anyway, note how we get an increasingly better approximation as the padding increases. To get the exact result, one would need an infinite amount of padding, which of course is unrealistic.

Test case #1 doesn't need padding, the function reaches zero at the edges of the sampled domain. We can impose such a behavior on the other cases too. In Discrete Fourier analysis this is called windowing. This would look something like this:

def fft_int(s, dydt):
    dydt_windowed = dydt * np.hanning(dydt.shape[0])
    f = fftfreq(dydt.shape[0], s)
    F = fft(dydt_windowed)
    with np.errstate(divide="ignore", invalid="ignore"):
        F = F / (2*np.pi*1j*f) 
    F[0] = 0
    y = np.real(ifft(F))
    return y

However, here we get correct integration results only in the middle of the domain, with increasingly suppressed values towards to ends. So this is not a practical solution either.

My conclusion is that no, this is not possible to do. It is much easier to compute the integral with np.cumsum:

yp = np.cumsum(dy) * s

plt.plot(y - np.mean(y))
plt.plot(yp - np.mean(yp))
plt.show()

(not showing output: the two plots overlap perfectly.)