3
$\begingroup$

I'm trying to learn how to compute accurate derivatives using the FFT. In the code at the end of this question I'm trying to compute derivatives of

$$ f(x) = \exp(-10(x-1)^2) ,\, \, x \in [0,2] $$

This works out because the function is infinitely differentiable, decays smoothly towards the end of the interval and can be extended periodically. However, when I try to compute derivatives of

$$ f(x) = 3x ,\, \, x \in [0,2] $$

I see oscillations, perhaps because the function exhibits a discontinuity when extended periodically. So, my question is, what can I do to be able to compute the derivative of $f(x)=3x$ accurately using the FFT? Perhaps zero-padding? Perhaps multiplication with a function which decays smoothly to zero the end of the interval and is 1.0 in the center? If you can point me to a treatment which explains with equations, even better.

Derivative of $f(x) = \exp(-10(x-1)^2)$

enter image description here

Derivative of $f(x) = 3x$

enter image description here

import numpy as np
import matplotlib.pyplot as plt

from scipy.fft import fft, ifft


nn = 128
LL = 1.0
dx = 2*LL/(nn-1)
xx = np.linspace(0.0,2*LL,nn,dtype='float64')


ff = np.exp(-10*(xx-LL)**2)
df = ff*(-20*(xx-LL))

#ff = 3.0*xx
#df = 3.0*np.ones_like(xx)

dfFD=np.zeros(len(df), dtype='complex_') 

# Compute central finite difference derivative
for kappa in range(1,len(df)-1):
    dfFD[kappa]=(ff[kappa+1]-ff[kappa-1])/(2.0*dx) 
dfFD[-1] = dfFD[-2]
dfFD[0]  = dfFD[1]


# Compute FFT derivative
fhat   = fft(ff)
kappa1 = (2*np.pi/(2*LL))*np.arange(-nn/2,nn/2)
kappa  = np.fft.fftshift(kappa1) 
dfhat  = kappa*fhat*(1j)


dfFFT  = np.real(ifft(dfhat))

plt.plot(xx, ff, color='r', linewidth=2, label='Function')
plt.plot(xx, df.real, color='k', linewidth=2, label='True Derivative')
plt.plot(xx, dfFD.real, 'x', color='b', linewidth=1.5, label='Finite Difference')
plt.plot(xx, dfFFT.real, 'o-', color='c', markerfacecolor='none',linewidth=1.5, label='Spectral Derivative')

plt.xlabel('X values')
plt.ylabel('Y values')
plt.legend()
plt.show()

plt.savefig('exponential.png')

Edit

Basically the idea suggested by @lightxbulb, I extended the function outside the interval and made it periodic.

enter image description here

$\endgroup$
3
  • 3
    $\begingroup$ Sample $3x$ in a larger window than $[0,1]$. Since $3x$ is not periodic, in the continuous case you would compute the Fourier transform of it and not its Fourier series, which means you would technically need to sample it in $(-\infty, \infty)$ in the discrete case. Practically you cannot do the latter, but you could extend it sufficiently and then ignore results near the boundary. $\endgroup$
    – lightxbulb
    Feb 20 at 12:13
  • $\begingroup$ Another option is to use a different spectral basis on $[0,1]$. Chebyshev polynomials are often used for nonperiodic functions, but of course you can no longer compute everything using an FFT $\endgroup$
    – whpowell96
    Feb 20 at 15:34
  • $\begingroup$ Seems arxiv.org/pdf/2211.15960.pdf has a method for extending the function periodically. I will investigate it. $\endgroup$
    – NNN
    Feb 21 at 7:21

2 Answers 2

10
$\begingroup$

There is little you can do about it beyond just using more terms in the Fourier expansion. What you observe is called Gibbs phenomenon. It, in essence, says that if you have a function that is discontinuous (like the periodic extension of $f(x)=x$ is), then the Fourier series (i) has over- and undershoots near the places of discontinuity, and (ii) these disappear only very slowly as you increase the number of terms in the Fourier series. Because you approximate the derivative of $f(x)$ by the derivative of the Fourier series with its over- and undershoots, you should expect the approximation to be poor -- which is exactly what you observe.

$\endgroup$
3
  • $\begingroup$ Thanks. I was hoping that there was something I could do. I have a situation where I need to compute spatial derivatives of the output of a physics based neural network using the FFT. The output is not guaranteed to be continuously periodically extendable. I was hoping that through zero padding, or reflecting the data or using some sort of mollifying function (something which decays to zero at the boundary, while being infinitely differentiable) I could compute the derivatives. $\endgroup$
    – NNN
    Feb 21 at 3:25
  • $\begingroup$ Just don't use the FFT for derivatives. Use a finite difference approximation. That's likely cheaper and more stable anyway. $\endgroup$ Feb 21 at 13:19
  • $\begingroup$ Yep, I'm going to go this route. $\endgroup$
    – NNN
    Feb 22 at 8:57
1
$\begingroup$

Savitzky-Golay filters work fairly well to approximate derivatives and a there's some cool work here too http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.