傅里叶变换浮点问题

Question

我正在实现一个传统的(这意味着不快)，分离傅里叶变换的图像。我知道，在浮点上，等间距样本中超过一个周期的sin或cos的和不是完全零的，这更像是常规变换的问题，而不是快速的问题。

该算法适用于二维双阵列，是正确的。逆在内部完成(在使用非对称公式时，通过双符号标志和条件检查)，而不是在外边使用共轭。结果与预期接近100%，所以这是一个关于细节的问题：

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

1. F(u，v) =求和(x=0->M-1)和(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 *(如上文所示)
2. F(u，v) = 1/sqrt(M*N) *(与上文类似) f(x，y) = 1/sqrt(M*N) *(如上文所示)

第一个是非对称变换对，第二个是对称变换对。对于非对称对，舍入误差更多地出现在图像的亮点中(一些像素在值范围(例如256)外稍微四舍五入)。对于对称对，误差更多出现在图像的恒定中程区域(不超过值范围！)。总的来说，对称对似乎会产生更多的舍入错误。

然后，它还依赖于输入:当图像存储在0,255中时，舍入错误不同于0,1。

那么，我的问题是:如何实现最优、最精确的算法(理论上，没有代码)：非对称/对称对？输入值范围为0,255或0,1？在将对数保存到文件中之前，线性提升结果如何？

编辑

我的算法只计算分离的非对称或对称的DFT公式。利用Eulers恒等式将因子分解为实部和虚部，然后分别展开和归纳为实部和虚部：

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 -

在我看来，这个值分组indside cos和sin应该给出最小的舍入误差(与例如cos(-mode*2*pi*v*y/N)相比)，因为不会多次增加/除以显著的假圆周率，但只有一次。难到不是么？

标度因子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 )。
• 输入inde0，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，M*N，在0,1输入不对称变换后，直接标度到0,255或线性缩放到0,255*M*N。

Answer 1

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

在FFT的库实现中(即编写的FFT例程通常由不同的应用程序使用，而不是为单个应用程序设计的自定义)，很少考虑缩放；该例程通常只是返回算法自然缩放的数据，不需要额外的乘法操作来调整比例。这是因为标度通常与应用程序无关(例如，找到具有最大能量的频率，无论尺度是什么)，或者标度可以通过乘运算分配，并且只执行一次(例如，应用程序可以通过显式缩放一次，而不是在正向变换和逆变换中进行缩放，从而获得同样的效果)。因此，由于通常不需要扩展，所以没有必要将其包含在库例程中。

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

关于使用什么变换(对数或线性)来显示光谱的问题，我不能建议；我不使用可视化光谱。

Answer 2

缩放会导致舍入错误。因此，解决方案1(只扩展一次)比解决方案2(实现两次)要好。类似地，在求和后一次缩放比在求和前缩放一切都要好。

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

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