6
$\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, 2020 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$ May 30, 2020 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, 2020 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, 2020 at 22:44

3 Answers 3

10
$\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 self-explanatory 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$
4
  • 1
    $\begingroup$ Your approach works, but I think that it would be good if you explain why it does. $\endgroup$
    – nicoguaro
    May 31, 2020 at 4:55
  • $\begingroup$ @nicoguaro I tried to give some more details in my answer. $\endgroup$
    – Socob
    Oct 28, 2020 at 14:30
  • $\begingroup$ How would your approach be used to obtain the gradient of a 2D function? $\endgroup$
    – Graham G
    Jul 31, 2022 at 10:36
  • $\begingroup$ @Graham G In 2D (or beyond) it is similar; basically you need to come up with a 2D Fourier interpolation of your data, and take the derivative of the Fourier series analytically. I suggest checking out this article first, en.wikipedia.org/wiki/… $\endgroup$ Jul 31, 2022 at 20:29
9
$\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$
6
  • $\begingroup$ I saw you used 2*pi to make sure your signal is periodic. If I want use the method to a signal (for instance a ramp signal plus a sine signal plus white noise) whose period is unknown, how can I avoid spectral leakage? $\endgroup$
    – John
    Feb 21 at 9:04
  • $\begingroup$ @John I’m not sure what exactly you mean – are you referring to point (1.)? The original code in the question already used 2*pi, so that’s not something I changed. It turns out that the function does have the period 2π, so that’s why it works well here. If you want to ask how to deal with spectral leakage in general, when the signal period is unknown (or perhaps the signal is not periodic at all), you should probably ask a separate question. $\endgroup$
    – Socob
    Feb 21 at 21:12
  • $\begingroup$ Using your code, I got inaccurate 2nd derivative result after changing 2*pi to 6. I wonder if there has to be an assumption that signal needs to have integer number of period so that this spectral quadratic weighting method (to get 2nd derivative) will work. $\endgroup$
    – John
    Feb 22 at 1:23
  • $\begingroup$ @John Yes, this is related to the spectral leakage I mentioned (even for the first derivative, there will be some deviation). Again, there are methods to deal with this, but comments are not the place to discuss completely separate questions like this. If you are interested, you should open a new question. $\endgroup$
    – Socob
    Feb 22 at 15:56
  • 1
    $\begingroup$ @John Yes it does, as I state in the text. If it didn’t, there wouldn’t have been any need to change x from a closed to an open interval in order to get the “correct” period. $\endgroup$
    – Socob
    Feb 23 at 22:08
1
$\begingroup$

Here we have $L=NT=2\pi$ (the total duration for which the signal was sampled), with the fundamental frequency $ω_o=\frac{2\pi}{NT}=\frac{2\pi}{L}=1$, slight modification of the code yields the correct derivative values computed with DFT (using the property $F(df/dx) = iωF(f)$).

L = 2*pi
dx = L/N 
x = linspace(-L/2, L/2-dx, N)

y = sin(2*x)+cos(5*x)
dydx = 2*cos(2*x)-5*sin(5*x)

k = (2*np.pi/L)*np.arange(-N/2, N/2)  # fundatmental frequency = 2π/L
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()

enter image description here

$\endgroup$

Your Answer

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

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