【问题标题】:Rounding errors giving incorrect tesults in DFT?在 DFT 中给出不正确结果的舍入错误?
【发布时间】:2013-10-09 16:13:01
【问题描述】:

我一直在这个 DFT 上撞墙。它应该打印出: 8,0,0,0,0,0,0,0 但我得到 8 然后是非常非常小的数字。这些是舍入错误吗?有什么我可以做的吗?我的 Radix2 FFT 给出了正确的结果,看起来愚蠢的 DFT 也不能工作。

我从复数开始,所以我知道有一点遗漏,我试图将其剥离以说明问题。

#include <cstdlib>
#include <math.h>
#include <iostream>
#include <complex>
#include <cassert>

#define SIZE 8
#define M_PI 3.14159265358979323846

void fft(const double src[], double dst[], const unsigned int n) 
{
    for(int i=0; i < SIZE; i++)
    {
        const double ph = -(2*M_PI) / n;
        const int gid = i;

        double res = 0.0f;
        for (int k = 0; k < n; k++) {

            double t = src[k];

            const double val = ph * k * gid;
            double cs = cos(val);
            double sn = sin(val);

            res += ((t * cs) - (t * sn));
            int a = 1;
        }

        dst[i] = res;
        std::cout << dst[i] << std::endl;
    }
}

int main(void)
{
    double array1[SIZE];
    double array2[SIZE];

    for(int i=0; i < SIZE; i++){
        array1[i] = 1;
        array2[i] = 0;
    }

    fft(array1, array2, SIZE);

    return 666;
}

【问题讨论】:

  • 小数字有多小?
  • @PaulR 在 1e-16 的顺序上
  • 这听起来很正确 - 16 或 17 位小数大约是双精度的限制。请注意,您应该真正使用 中的 M_PI 而不是自己滚动,但我认为在这种情况下并没有太大区别。
  • 是的,这真的很烦人。我可以使用 std::cout.setf ( std::ios::fixed, std::ios::floatfield ) 修复输出,但这只能修复显示。在大型数据集上,累积误差实际上加起来是有意义的。

标签: c++ fft rounding dft


【解决方案1】:

FFT 实际上可以产生比直接 DFT 计算更准确的结果,因为较少的算术运算通常允许较少的算术量化误差累积机会。 FFTW 的一位作者有一篇关于此主题的论文。

由于 DFT/FFT 处理超越基函数,因此使用任何非符号和有限计算机数字格式的结果将永远不会完全正确(可能在少数特殊情况下,或幸运的意外除外)。因此,与零非常接近(在几个 LSB 内)的值应该被简单地作为噪声忽略,或者认为与零相同。

【讨论】:

  • 是的,我现在正在做:int tmp = res*100; dst[i] = tmp /100.0f;有什么更有效的吗?我也有一个基数 2 FFT,但它们只适用于 2^N 大小的数据集,这对我没有多大好处。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2012-08-27
  • 2017-12-14
  • 2019-03-12
  • 1970-01-01
  • 1970-01-01
  • 2023-03-29
相关资源
最近更新 更多