因此,我想像我之前的许多人一样,仔细筛选通过傅立叶变换工作进行微分所必需的所有相移和归一化系数的天体对准。
我正在尝试使用尽可能少的代码,在使用其数组操作时充分利用numpy功能保持代码干净,据我了解,这比显式python循环要快。
从一维函数关系表的第106行中,在有关傅立叶变换的Wikipedia文章中有关Important Fourier Transforms的部分中,我得到了傅立叶变换具有的这种属性,其中实空间的微分等于频率空间的乘法这样的d^nf/dx^n = ifft[F*(i*w)^n]
,其中F
是f
的傅立叶变换。
在寻找更坚实的参考来定义w
在d^nf/dx^n = ifft[F*(i*w)^n]
中的含义时,我找到了一篇有关使用DFT进行区分的数学的PDF论文。在本文中,作者讨论了给定函数的无限多个所谓的三角插值,并展示了一种获得“较少摆动”解的方法,这是唯一的。
因此,我使用了他的定义并编写了一个小的python代码进行区分。我出现的代码就是这一行(请注意最后两行,因为它们会更改)
N = 51
a = 0
b = 2*np.pi
X = np.linspace(a,b, N)
L = abs(b-a)
TwoPiOverL = 2.0*np.pi/L
freqs = np.fft.fftfreq(N,1./N)
# THESE COME FROM [1] #
D1 = TwoPiOverL*freqs*1j
D2 = -(TwoPiOverL*freqs)**2
#=========================#
S = np.sin(X)
fftS=np.fft.fft(S)
C = np.cos(X)
# COMPUTING THE DERIVATIVES #
Y1 = np.fft.ifft(D1*fftS).real
Y2 = np.fft.ifft(D2*fftS).real
在生成的图形的图中,符号
~d/dx sin
表示根据傅立叶变换乘法评估的正弦函数的导数,而不是解析表达式。测试1
首次测试,仅使用PDF中的定义:
因此,第一个测试是一团糟。导数在域的中间进行很好的插值,但是端点很疯狂。在PDF本身中,第6节开始时,作者说:
以上所有内容都适用于周期函数y(x)的谱微分。如果函数是非周期性的,
由于端点处的隐式不连续性,伪像将出现在导数中。
听起来不错,这意味着我正在搞砸,因为如果整篇文章都写成讨论周期性信号,并且如果我的信号碰巧属于非周期性类,那么如果我可以使用周期性信号。
但是我正在使用DFT,这意味着我正在对数组应用数字过程。离散数据,不包含有关周期性的信息。所以我想:
好了,我怎么能告诉我的傅立叶微分器信号是周期性的?这在代码中意味着什么?我必须在代码中的哪里更改才能将该信息传达给我的程序?
似乎一个答案来自stackoverflow。在关于使用FFT进行数值微分的this post中(是的,冗余,但我会到达那里),@ dkato提供了我已解释为将代码中的最后两行从
Y1 = np.fft.ifft(D1*fftS).real
Y2 = np.fft.ifft(D2*fftS).real
至
Y1 = np.fft.ifft(np.exp(D1)*fftS).real
Y2 = np.fft.ifft(np.exp(D2)*fftS).real
所以我做到了,然后测试了我会得到什么。
测试2
第二次测试,使用@dkato更正:
从上图可以看出,通过使用@dkato提供的校正,可以消除端点摆动,但是对导数的估计仍然相距甚远。我正在使用51个数据点,所以我猜测这不是由于扩展不准确所致。看来我计算导数的方法不正确。
解析的一阶导数与其傅里叶变换估计之间准确性损失最大的原因似乎是向右偏移,而二阶导数不仅偏移而且缩放比例也很差。有趣的是,当我们使用
Y2 = np.fft.ifft(-np.exp(D2)*fftS).real
要评估二阶导数,我们看到all of the error due to shifting goes away, leaving only the scaling。我不知道这是否是一般FFT驱动的二阶导数的属性,或者仅是因为我正使用Sine函数作为y(x)。
最后一点
我很高兴端点上的怪异,丑陋的摆动消失了,这意味着我必须以某种方式告诉我的傅立叶微分器,我正在使用周期性信号。但是我的微分估计的这些奇怪的变化和标度仍然令人不安,而且我一生都无法专心于错误所在。您可能已经猜到了,部分原因是我对DFT中所有活动部件的使用经验很少。
以这个身份,如果您愿意帮助我解决这个问题,我对社区有两个问题。
1)如果我确实以某种方式告诉傅立叶微分器我正在使用周期函数,那我是怎么做到的呢?起作用的代码和不起作用的代码背后的原因是什么?
和
2)分析和估计导数之间这些怪异的比例/位移不一致性的来源是什么?我该如何解决?
最佳答案
让我们首先回答问题的症结:
分析和估计的导数之间这些怪异的比例/位移不一致性的来源是什么?我该如何解决?
问题从以下几行开始:
X = np.linspace(a,b, N)
默认情况下,
numpy.linspace
包括端点,由于端点被重复,因此端点在计算的周期函数扩展中引入了不连续性。为了说明这可能是一个问题,您可以查看
sin(X)
的计算函数N=4
:sin(0)
sin(2*np.pi/3)
sin(4*np.pi/3)
sin(2*np.pi)
其相位在每个后续值处增加
2*np.pi/3
。重复这些值以获得周期为N=4
的周期性扩展,将得出sin(0)
sin(2*np.pi/3)
sin(4*np.pi/3)
sin(2*np.pi)
sin(0) # == sin(2*np.pi)
sin(2*np.pi/3) # == sin(2*np.pi + 2*np.pi/3)
sin(4*np.pi/3) # == sin(2*np.pi + 4*np.pi/3)
sin(2*np.pi) # == sin(2*np.pi + 2*np.pi)
您可以在每个步骤中看到相位以
2*np.pi/3
递增,但在块边界处保持不变。这种小的不连续性使得用三角插值基础很难很好地逼近。然后,所得的内插包含明显的振荡,该振荡会通过推导操作进一步放大。要解决此问题,您应该将终结点与排除在一起(否则保留您的初始解决方案不变):
X = np.linspace(a,b, N, endpoint=False)
就@dkato's post中的解决方案而言,它似乎仅实现了时移操作。对于看起来不太差的正弦信号,但与微分无关。
最后,请回答部分
如果确实以某种方式告诉我的傅立叶微分器我正在使用周期函数,那我该怎么做呢?
是的,这是真的,您在开始使用离散傅立叶变换(及其FFT实现)的那一刻就做到了。
关于python - 如何解决由FFT驱动的微分程序中的移位和缩放错误?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/52224285/