傅里叶变换浮点问题答案

【问题标题】：Fourier transform floating point issues傅里叶变换浮点问题
【发布时间】：2013-12-12 18:46:54
【问题描述】：

我正在为图像实现一个传统的（这意味着不快）分离傅里叶变换。我知道在浮点中，等间隔样本中的一个 sin 或 cos 周期的总和不是完全为零，这对于传统变换来说比快速变换更成问题。

该算法适用于二维双数组并且是正确的。反转是在内部完成的（在使用不对称公式时通过双符号标志和条件检查），而不是在外部进行共轭。结果几乎 100% 符合预期，所以这是一个关于细节的问题：

当我执行正向变换、将对数幅值和角度保存到图像、重新加载它们并进行逆变换时，我会遇到不同类型的舍入误差以及不同类型的实现公式：

F(u,v) = Sum(x=0->M-1) Sum(y=0->N-1) f(x,y) * e^(-i*2*pi*u *x/M) * e^(-i*2*pi*v*y/N)

f(x,y) = 1/M*N *（如上）
F(u,v) = 1/sqrt(M*N) *（如上）

f(x,y) = 1/sqrt(M*N) * （如上）

所以第一个是不对称的变换对，第二个是对称的。对于非对称对，舍入误差更多地出现在图像的亮点中（某些像素在值范围之外略微舍入（例如 256））。使用对称对，误差更多地出现在图像的恒定中间范围区域（没有超出值范围！）。总的来说，对称对似乎产生了更多的舍入误差。

然后，它还取决于输入：当图像存储在 [0,255] 中时，舍入误差不同于在 [0,1] 中时的舍入误差。

所以我的问题是：应该如何实现最优、最准确的算法（理论上，没有代码）：非对称/对称对？ [0,255] 或 [0,1] 中的输入值范围？在将对数保存到文件之前如何线性放大结果？

编辑：

我的算法只是计算分离的非对称或对称 DFT 公式。使用欧拉恒等式将因子分解为实部和虚部，然后分别展开总结为实部和虚部：

sum_re += f_re * cos(-mode*pi*((2.0*v*y)/N)) - // mode = 1 for forward, -1
          f_im * sin(-mode*pi*((2.0*v*y)/N));  // for inverse transform
// sum_im permutated in the known way and + instead of -

这个值在 cos 和 sin 内分组在我看来应该是最小的舍入误差（与例如 cos(-mode*2*pi*v*y/N) 相比），因为不会多次乘/除显着错误的舍入超越 pi，但只有一次。不是吗？

比例因子1/M*N 或1/sqrt(M*N) 分别在最内层和的每个分隔外部之后应用。里面好点？还是在两次分离结束时完全合并？

为了进行更深入的分析，我已经退出了input->transform->save-to-file->read-from-file->transform^-1->output 工作流程并选择直接以双精度进行比较：input->transform->transform^-1->output。

这里是真实生活中 704x528 8 位图像的结果（增量 = 输入和输出的实部之间的最大绝对差）：

在 [0,1] 内输入和不对称公式：delta = 2.6609e-13（对应于 [0,255] 范围的 6.785295e-11）。
输入 insde [0,1] 和对称公式：delta = 2.65232e-13（对应于 [0,255] 范围的 6.763416e-11）。
在 [0,255] 内输入和不对称公式：delta = 6.74731e-11。
在 [0,255] 内输入和对称公式：delta = 6.7871e-11。

这些并不是真正的显着差异，但是，具有非对称变换的全范围输入表现最佳。我认为 16 位输入的值可能会变差。

但总的来说，我发现，我遇到的问题更多是因为缩放前保存到文件（或反向）舍入错误，而不是真正的转换舍入错误。

但是，我很好奇：傅里叶变换最常用的实现是什么：对称还是非对称？哪个值范围通常用于输入：[0,1] 或 [0,255]？通常以对数刻度显示光谱：例如[0,1]输入的非对称变换后的[0,M*N]直接对数缩放到[0,255]还是线性缩放到[0,255*M*N]之前？

【问题讨论】：

如果您使用相同的核心算法并且只是对结果进行不同的缩放，那么错误的差异应该是偶然的，本质上只是关于舍入如何发生的随机机会。我们需要您正在做的调查的详细信息。
@EricPostpischil：查看我的编辑

标签： c++ math floating-point fft

【解决方案1】：

您报告的错误很小、很正常，通常可以忽略。只需缩放您的结果并将目标区间之外的任何结果限制到端点。

在 FFT 的库实现中（即，编写的 FFT 例程通常由不同的应用程序使用，而不是为单个应用程序定制设计），很少考虑缩放；该例程通常只是简单地返回已通过算术自然缩放的数据，而无需使用额外的乘法运算来调整比例。这是因为比例通常与应用程序无关（例如，无论比例是多少，找到具有最大能量的频率都有效）或者比例可以通过乘法运算分布并且只执行一次（例如，而不是缩放在正向变换和逆变换中，应用程序可以通过显式缩放一次来获得相同的效果）。因此，由于通常不需要缩放，因此将其包含在库例程中是没有意义的。

数据缩放到的目标间隔取决于应用程序。

关于使用什么变换（对数或线性）来显示光谱的问题，我无法给出建议；我不使用可视化光谱。

【讨论】：

【解决方案2】：

缩放会导致舍入错误。因此，解决方案 1（缩放一次）优于解决方案 2（缩放两次）。同样，在求和之后缩放一次比在求和之前缩放所有内容要好。

您是从 0 运行 y 到 2*N 还是从 -N 到 +N ？在数学上它是相同的，但在后一种情况下你有额外的精度。

顺便说一句，mode 在 cos(-mode * stuff) 中做什么？

【讨论】：

y 从 0 运行到 N（不是 2*N）。是的，确实，我只能从 0 运行到 N/2（如果 N 偶数，否则 floor(N/2) 需要额外处理 ceil(N/2)），而不是 compute 而只是共轭剩余的y。也许这会更准确，我会检查一下。 mode 是变换的模式：如果 1 则为正向，如果 -1 则为逆向（正向和逆变换的公式不同（除了非对称情况的因子）仅由那个符号）。
在 sum 中进行缩放：我可以想象，即使存在绝对更多的舍入误差，它们也会补偿，因为 sin/cos 从 -1 波动到 1。
@mb84：这是不平凡的数学，但我相信舍入误差不会在这里取消。基本上，它们彼此不相关，但每个都与它们各自术语的大小相关。由于这些项本身取消，因此与单个项相比，总和相对较小（DC 除外，当然 v=0）。这意味着您要添加大量相当大的错误。虽然其预期误差将为0，但标准差将类似于O(sqrt(N))更大。
@mb84：模式注释是数学的：cos 甚至如此 cos(-mode*stuff) 只是 cos(stuff)。再说一次，你的编译器可能足够聪明，可以在适当的地方使用融合的 sincos 函数。
啊好的，现在我明白你对mode的评论了。由于可读性，我喜欢它所写的内容。而且我注意到我不能使用 F(u,v)==conjugate(F(-u,-v)) 属性，因为我想允许整个复杂的输入。关于roundig错误，我将尽可能少地乘以常数。与 [0,255] 相比，错误确实微不足道。它仍然是我感兴趣的频谱的对数标度：一般来说它是如何完成的（请参阅我第一篇文章中的最后一个问题）？