在Python中使用DFT查找一阶导数 [英] Finding first derivative using DFT in Python

问题描述

我想使用离散傅里叶变换在区间[0, 2/pi]上找到exp(sin(x))的一阶导数.基本思想是，首先在给定的时间间隔内评估exp(sin(x))的DFT，给出v_k，然后计算ikv_k的反DFT，为您提供所需的答案.实际上，由于以编程语言实现了Fourier变换，因此您可能需要在某处对输出进行重新排序和/或在此处与此处乘以不同的因子.

I want to find the first derivative of exp(sin(x)) on the interval [0, 2/pi] using a discrete Fourier transform. The basic idea is to first evaluate the DFT of exp(sin(x)) on the given interval, giving you say v_k, followed by computing the inverse DFT of ikv_k 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.

我首先在Mathematica中做到了，那里有一个选项FourierParameters，它使您可以为变换指定一个约定.首先，我获得了一个高斯的傅里叶级数，以查看我必须乘以什么然后再求导数的归一化因子.不幸的是，此后将我的Mathematica代码转换成Python(由此我再次进行了高斯的傅里叶级数-这是成功的)，但我没有得到相同的结果.这是我的代码:

I first did it in Mathematica, where there is an option FourierParameters, 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:

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()

如果人们希望我可以发布Mathematica代码.

I can post the Mathematica code, if people want me to.

推荐答案

原来的问题是np.zeros为您提供了一个实零数组，而不是复杂的零数组，因此此后的赋值不会更改任何内容，例如他们是虚构的.

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

因此解决方案非常简单

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

从而将np.zeros替换为1j*np.zeros

这篇关于在Python中使用DFT查找一阶导数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！