python - 如何修复 FFT 驱动的微分程序中的移位和缩放错误？

Question

因此，我正在尝试，像我之前的许多其他人一样，筛选通过傅立叶变换工作进行微分所需的所有相移和归一化系数的天体对齐。

我试图使用尽可能少的代码，在使用数组操作时完全使用 numpy 功能来保持代码干净，有人告诉我，这比显式 python 循环更快。

从函数关系表的第 106 行，一维，在关于傅里叶变换的维基百科文章的重要傅里叶变换部分中，我得到傅里叶变换具有这个属性，其中实空间中的微分等同于频率空间中的乘法这样，的傅里叶变换d^nf/dx^n = ifft[F*(i*w)^n]在哪里。Ff

为了寻找更可靠的参考来定义将是什么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 进行数值微分的帖子中（是的，多余的，但我会到达那里），@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

为了评估二阶导数，我们看到由于移位引起的所有误差都消失了，只剩下缩放。我不知道这是否是一般 FFT 驱动的二阶导数的属性，或者仅仅是因为我使用正弦函数作为我的 y(x)，特别是。

重点，终于

我很高兴端点处奇怪、丑陋的摆动消失了，这意味着我必须以某种方式告诉我的傅立叶微分器我正在处理周期性信号。但是我的微分估计的这些奇怪的变化和缩放仍然很麻烦，而且我一辈子都无法指出哪里出了问题。正如您可能已经猜到的那样，部分原因是我几乎没有玩过 DFT 中所有移动部件的经验。

以这种身份，如果您愿意帮助我解决这个问题，我有两个问题要问社区。

1）如果我确实以某种方式告诉我的傅立叶微分器我正在使用周期函数，那么我到底是怎么做到的？有效代码和无效代码背后的原因是什么？

和

2) 分析导数和估计导数之间这些奇怪的尺度/位移不协调的根源是什么？我该如何修复它们？

标签： pythonnumpyfftdiscrete-mathematicsdifferentiation

Answer 1

我们先回答问题的症结所在：

分析导数和估计导数之间这些奇怪的尺度/偏移不协调的根源是什么？我该如何修复它们？

问题从以下行开始：

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在每个后续值处递增。重复这些值以获得周期性的扩展 period 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 帖子中的解决方案，看起来它只实现了时移操作。对于可能看起来不太糟糕但与微分无关的正弦信号。

最后，回答部分

如果我确实以某种方式告诉我的傅立叶微分器我正在使用周期函数，那么我到底是怎么做到的？

是的，这是真的，您在开始使用离散傅里叶变换（及其 FFT 实现）的那一刻就做到了。