Finding first derivative using DFT in Python
<p>I want to find the first derivative of <code>exp(sin(x))</code> on the interval <code>[0, 2/pi]</code> using a discrete Fourier transform. The basic idea is to first evaluate the DFT of <code>exp(sin(x))</code> on the given interval, giving you say <code>v_k</code>, followed by computing the inverse DFT of <code>ikv_k</code> giving you the desired answer. In reality, due to the implementations of Fourier transforms in programming languages, you might need to reorder the output somewhere and/or multiply by different factors here and there.</p>

<p>I first did it in Mathematica, where there is an option <code>FourierParameters</code>, which enables you to specify a convention for the transform. Firstly, I obtained the Fourier series of a Gaussian, in order to see what the normalisation factors are that I have to multiply by and then went on finding the derivative. Unfortunately, translating my Mathematica code into Python thereafter (whereby again I first did the Fourier series of a Gaussian - this was successful), I didn't get the same results. Here is my code:</p>

<pre><code>N=1000
xmin=0
xmax=2.0*np.pi
step = (xmax-xmin)/(N)
xdata = np.linspace(xmin, xmax-step, N)
v = np.exp(np.sin(xdata))
derv = np.cos(xdata)*v
vhat = np.fft.fft(v)
kvals1 = np.arange(0, N/2.0, 1)
kvals2 = np.arange(-N/2.0, 0, 1)
what1 = np.zeros(kvals1.size+1)
what2 = np.empty(kvals2.size)
it = np.nditer(kvals1, flags=['f_index'])
while not it.finished:
    np.put(what1, it.index, 1j*(2.0*np.pi)/((xmax-xmin))*it[0]*vhat[[int(it[0])]])
    it.iternext()
it = np.nditer(kvals2, flags=['f_index'])
while not it.finished:
    np.put(what2, it.index, 1j*(2.0*np.pi)/((xmax-xmin))*it[0]*vhat[[int(it[0])]])
    it.iternext()
xdatafull = np.concatenate((xdata, [2.0*np.pi]))
what = np.concatenate((what1, what2))
w = np.real(np.fft.ifft(what))

fig = plt.figure()
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))

plt.plot(xdata, derv, color='blue')
plt.plot(xdatafull, w, color='red')
plt.show()
</code></pre>

<p>I can post the Mathematica code, if people want me to.</p>

<p>Turns out the problem is that <code>np.zeros</code> gives you an array of real zeroes and not complex ones, hence the assignments after that don't change anything, as they are imaginary.</p>

<p>Thus the solution is quite simply</p>

<pre><code>import numpy as np

N=100
xmin=0
xmax=2.0*np.pi
step = (xmax-xmin)/(N)
xdata = np.linspace(step, xmax, N)
v = np.exp(np.sin(xdata))
derv = np.cos(xdata)*v
vhat = np.fft.fft(v)

what = 1j*np.zeros(N)
what[0:N/2.0] = 1j*np.arange(0, N/2.0, 1)
what[N/2+1:] = 1j*np.arange(-N/2.0 + 1, 0, 1)
what = what*vhat
w = np.real(np.fft.ifft(what))

# Then plotting
</code></pre>

<p>whereby the <code>np.zeros</code> is replaced by <code>1j*np.zeros</code></p>

python fft derivative