3
$\begingroup$

I wrote the following code to compute the approximate derivative of a function using FFT:

from scipy.fftpack import fft, ifft, dct, idct, dst, idst, fftshift, fftfreq
from numpy import linspace, zeros, array, pi, sin, cos, exp
import matplotlib.pyplot as plt

N = 100
x = linspace(0,2*pi,N)

dx = x[1]-x[0]
y = sin(2*x)+cos(5*x)
dydx = 2*cos(2*x)-5*sin(5*x)

k = fftfreq(N,dx)
k = fftshift(k)

dydx1 = ifft(-k*1j*fft(y)).real

plt.plot(x,dydx,'b',label='Exact value')
plt.plot(x,dydx1,'r',label='Derivative by FFT')
plt.legend()
plt.show()

However, it is giving unexpected results, which I believe is related to the incorrect input of the wavenumbers given by the array k:

Comparsion exact and numeric derivative

I know that different implementations of the FFT handle the wavenumbers order differently, so what am I missing here? Any ideas will be very appreciated.

$\endgroup$
4
  • $\begingroup$ Welcome to Scientific Computing SE. Can you please edit your question to elaborate why you want to use the FFT to compute the derivative? This seems like a very roundabout way of doing this. $\endgroup$ – Wrzlprmft May 30 '20 at 14:44
  • $\begingroup$ I must implement a FFT solver for the Poisson equation, however I must be able to solve a simpler problem like this one first. $\endgroup$ – Leonardo Araujo May 30 '20 at 14:49
  • $\begingroup$ If you run through the error analysis of the FFT, you'll see that this is an inaccurate way to compute the numerical derivative. B-splines have better spectral properties for numerical differentiation. $\endgroup$ – user14717 May 30 '20 at 18:16
  • 2
    $\begingroup$ You can't expect accurate results from the FFT if your function is non-periodic. Hoeever, ig it is periodic and, moreover, smooth, it will yield exponential accuracy for the derivative. Your test function is suitable for spectral differentiation, so it seems to be an implementation issue. This is properly adressed in the answer by @MaximUmansky. $\endgroup$ – davidhigh Jun 3 '20 at 22:44
8
$\begingroup$

FFT returns a complex array that has the same dimensions as the input array. The output array is ordered as follows:

  1. Element 0 contains the zero frequency component, F0.
  2. The array element F1 contains the smallest, nonzero positive frequency, which is equal to 1/(Ni Ti), where Ni is the number of elements and Ti is the sampling interval.
  3. F2 corresponds to a frequency of 2/(Ni Ti).
  4. Negative frequencies are stored in the reverse order of positive frequencies, ranging from the highest to lowest negative frequencies.
  5. For an even number of points, the frequencies corresponding to the returned complex values are: 0, 1/(NiTi), 2/(NiTi), ..., (Ni/2–1)/(NiTi), 1/(2Ti), –(Ni/2–1)/(NiTi), ..., –1/(NiTi) where 1/(2Ti) is the Nyquist critical frequency.

  6. For an odd number of points, the frequencies corresponding to the returned complex values are: 0, 1/(NiTi), 2/(NiTi), ..., (Ni–1)/2)/(NiTi), –(Ni–1)/2)/(NiTi), ..., –1/(NiTi)

Using this information we can construct the proper vector of frequencies that should be used for calculating the derivative. Below is a piece of Python code that does it all correctly. Note that the factor 2$\pi$N cancels out due to normalization of FFT.

from scipy.fftpack import fft, ifft, dct, idct, dst, idst, fftshift, fftfreq
from numpy import linspace, zeros, array, pi, sin, cos, exp, arange
import matplotlib.pyplot as plt


N = 100
x = 2*pi*arange(0,N,1)/N #-open-periodic domain                                                   

dx = x[1]-x[0]
y = sin(2*x)+cos(5*x)
dydx = 2*cos(2*x)-5*sin(5*x)


k2=zeros(N)

if ((N%2)==0):
    #-even number                                                                                   
    for i in range(1,N//2):
        k2[i]=i
        k2[N-i]=-i
else:
    #-odd number                                                                                    
    for i in range(1,(N-1)//2):
        k2[i]=i
        k2[N-i]=-i

dydx1 = ifft(1j*k2*fft(y))

plt.plot(x,dydx,'b',label='Exact value')
plt.plot(x,dydx1, color='r', linestyle='--', label='Derivative by FFT')
plt.legend()
plt.show()

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Your approach works, but I think that it would be good if you explain why it does. $\endgroup$ – nicoguaro May 31 '20 at 4:55
  • $\begingroup$ @nicoguaro I tried to give some more details in my answer. $\endgroup$ – Socob Oct 28 '20 at 14:30
4
$\begingroup$

Maxim Umansky’s answer describes the storage convention of the FFT frequency components in detail, but doesn’t necessarily explain why the original code didn’t work. There are three main problems in the code:

  1. x = linspace(0,2*pi,N): By constructing your spatial domain like this, your x values will range from $0$ to $2\pi$, inclusive! This is a problem because your function y = sin(2*x)+cos(5*x) is not exactly periodic on this domain ($0$ and $2\pi$ correspond to the same point, but they’re included twice). This causes spectral leakage and thus a small deviation in the result. You can avoid this by using x = linspace(0,2*pi,N, endpoint=False) (or x = 2*pi*arange(0,N,1)/N, as Maxim Umansky did; this is what he is referring to with “open-periodic domain”).
  2. k = fftshift(k): As Maxim Umansky explained, your k values need to be in a specific order to match the FFT convention. fftshift sorts the values (from small/negative to large/positive), which can be useful e. g. for plotting, but is incorrect for computations.
  3. dydx1 = ifft(-k*1j*fft(y)).real: scipy defines the FFT as y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum(), i. e. with a factor of $2\pi$ in the exponential, so you need to include this factor when deriving the formula for the derivative. Also, for scipy’s FFT convention, the k values shouldn’t get a minus sign.

So, with these three changes, the original code can be corrected as follows:

from scipy.fftpack import fft, ifft, dct, idct, dst, idst, fftshift, fftfreq
from numpy import linspace, zeros, array, pi, sin, cos, exp
import matplotlib.pyplot as plt

N = 100
x = linspace(0,2*pi,N, endpoint=False) # (1.)

dx = x[1]-x[0]
y = sin(2*x)+cos(5*x)
dydx = 2*cos(2*x)-5*sin(5*x)

k = fftfreq(N,dx)
# (2.)

dydx1 = ifft(2*pi*k*1j*fft(y)).real # (3.)

plt.plot(x,dydx,'b',label='Exact value')
plt.plot(x,dydx1,'r',label='Derivative by FFT')
plt.legend()
plt.show()
$\endgroup$
1
  • $\begingroup$ Good explanations! $\endgroup$ – Maxim Umansky Oct 28 '20 at 14:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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