本文介绍了用Numpy.fft实现傅里叶空间的数值积分的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我想用傅里叶空间中的数值积分对函数进行积分。
以下代码显示了一个工作示例：
import numpy as np
from pylab import *
from numpy.fft import fft, ifft, fftshift, ifftshift
N = 2**16
x = np.linspace(- np.pi , np.pi,N)
y = np.exp(-x**2) # function f(x)
ys = np.exp(-x**2) * (-2*x) # derivative f'(x)
T = x[-1] - x[0] # the whole range
w = (np.arange(N) - N /2.) / T + 0.00000001 # slightly shifted
# integration
fourier = ifft(ifftshift(1./ ( 2 * np.pi * 1j * w) * fftshift(fft(ys)) ) )
# differentiation
fourier2 = ifft(ifftshift(( 2 * np.pi * 1j * w) * fftshift(fft(fourier)) ) )
+ 0.00000001w。我需要它，因为否则我将生成ZeroDivisionError或麻木警告。这是一种变通方法，对于上面的例子似乎还可以，但对于我遇到的一个更复杂的问题，它却失败了。一位同事告诉我，如果我得到移频值的FFT(np.arange(N) - N /2. + 1./2) / T，我就可以简单地避免它。如何在麻木中做到这一点？有没有办法指定NumPy FFT的输出网格？
谢谢！
推荐答案
问题是w包含0(这是应该的)，并且除以w。w中的0是"DC"频率；它对应于傅里叶级数的常量项。
w时，将得到(0+0j)/0，即(nan+nanj)。如果输入的DC分量不为零，则会得到(inf+nanj)。无论采用哪种方法，解决方案都是忽略您得到的任何内容，并在使用ifft进行反转之前将DC傅立叶系数设置为0。
有几种方法可以实现这一点。一种方法是更改以下行：
w = (np.arange(N) - N /2.) / T + 0.00000001 # slightly shifted
# integration
fourier = ifft(ifftshift(1./ ( 2 * np.pi * 1j * w) * fftshift(fft(ys)) ) )
到这里(我添加了几个中间变量)：
w = (np.arange(N) - N /2.) / T
# integration
Fys = fft(ys)
with np.errstate(divide="ignore", invalid="ignore"):
modFys = ifftshift(1./ (2 * np.pi * 1j * w) * fftshift(Fys))
# modFys[0] will hold the result of dividing the DC component of y by 0, so it
# will be nan or inf. Setting modFys[0] to 0 amounts to choosing a specific
# constant of integration.
modFys[0] = 0
fourier = ifft(modFys).real
我也接受了ifft结果的真实部分。理论上，虚部应该都是0；实际上，由于正常的不精确浮点运算，虚部会很小，但不是零。
顺便说一句，如果您不想实现此技术的您自己的版本，可以使用scipy.fftpack.diff。
这篇关于用Numpy.fft实现傅里叶空间的数值积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持编程学习网！
本站部分内容来源互联网,如果有图片或者内容侵犯您的权益请联系我们删除！