0

好的，我要实现的目标很简单。对一些随机数据应用 FFT，然后在输出上应用反向算法以取回输入。我正在kissFFT为此使用图书馆。

代码：

const int fft_siz = 512;
const int inverse = 1;

kiss_fft_cpx* in = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * fft_siz);
kiss_fft_cpx* out = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * fft_siz);
kiss_fft_cpx* rec = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * fft_siz);
kiss_fft_cfg cfg = kiss_fft_alloc(fft_siz, !inverse, NULL, NULL);
kiss_fft_cfg icfg = kiss_fft_alloc(fft_siz, inverse, NULL, NULL);

srand((unsigned int)time(NULL));
for(int i = 0; i < fft_siz; i++)
{
    in[i].r = rand() % 256;
    in[i].i = rand() % 256;
}

kiss_fft(cfg, in, out);

// scaling
for(int i = 0; i < fft_siz; i++)
{
    out[i].r /= fft_siz;
    out[i].i /= fft_siz;
}

kiss_fft(icfg, out, rec);

unsigned int count = 0;
for(int i = 0; i < fft_siz; i++)
    if(in[i].r != rec[i].r)
    {
        count++;
        printf( "in[%3d].r does not match rec[%3d].r :: %3d :: %f\n",
                i, i, count, in[i].r - rec[i].r);
    }
    else if(in[i].i != rec[i].i)
    {
        count++;
        printf( "in[%3d].i does not match rec[%3d].i :: %3d :: %f\n",
                i, i, count, in[i].i - rec[i].i);
    }

free(in);
free(out);
free(rec);
free(cfg);
free(icfg);

kiss_fft_cleanup();

输出：

in[  0]:     71.000000       85.000000 -- out[  0]:    127.095703      124.541016
in[  1]:    248.000000       27.000000 -- out[  1]:     -7.083314        0.072701
in[  2]:     64.000000       18.000000 -- out[  2]:     -3.770610        2.682554
in[  3]:      6.000000       96.000000 -- out[  3]:     -7.929140       -2.897723
in[  4]:     98.000000       23.000000 -- out[  4]:     -0.719621       -5.854260
in[  5]:    250.000000      188.000000 -- out[  5]:      0.397226       -1.248124
in[  6]:    231.000000        3.000000 -- out[  6]:     -7.934285       -2.367196
in[  7]:      6.000000      105.000000 -- out[  7]:     -0.317480       -2.955601
in[  8]:    172.000000      143.000000 -- out[  8]:     -4.236186        3.911616
in[  9]:     16.000000      134.000000 -- out[  9]:     -0.162577       -5.353521
in[ 10]:    230.000000      112.000000 -- out[ 10]:     -4.703711        7.791993
in[ 11]:      5.000000       26.000000 -- out[ 11]:     -2.636305        0.188381
in[ 12]:     16.000000      127.000000 -- out[ 12]:      1.137413        4.576081
in[ 13]:    112.000000       86.000000 -- out[ 13]:      0.978051       -0.408992
in[ 14]:     40.000000       23.000000 -- out[ 14]:      5.231920       -2.347566
in[ 15]:     75.000000       26.000000 -- out[ 15]:      0.009981       -2.091559

note                                ::count::difference
--------------------------------------------------------
in[  1].r does not match rec[  1].r ::   1 :: -0.000031
in[  3].r does not match rec[  3].r ::   2 :: -0.000015
in[  4].i does not match rec[  4].i ::   3 :: -0.000004
in[  6].i does not match rec[  6].i ::   4 :: -0.000008
in[  7].r does not match rec[  7].r ::   5 :: -0.000002
in[  9].r does not match rec[  9].r ::   6 :: -0.000015
in[ 11].r does not match rec[ 11].r ::   7 :: -0.000015
in[ 12].r does not match rec[ 12].r ::   8 :: -0.000015
in[ 13].i does not match rec[ 13].i ::   9 :: -0.000008
in[ 14].i does not match rec[ 14].i ::  10 :: 0.000008
in[ 15].r does not match rec[ 15].r ::  11 :: -0.000015

调试。如果你走到底部，你会看到有 317 个不匹配。我还输出值之间的差异，即(in[].r - rec[].r)或(in[].i - rec[].i)。

接下来我要展示的是输入数据，其中白点代表实部，红点代表虚部。

在此处输入图像描述

这是用紫色表示的 FFT 的输出数据以及用白色和红色表示的重建数据。

在此处输入图像描述

在此处输入图像描述

注意到细微的差别了吗？我猜这与浮点精度有关。我怎样才能克服这个问题来获得与我使用 FFT 完全相同的输入数据？

编辑：

我注意到在我的情况下，范围关闭错误类似于 ]0 , 0.002]。因此，作为一种解决方法，我对重建的数据进行了四舍五入，并得到了一个很好的结果。但是...这仅在我的数字的小数部分为 0.0 时才有效。

4

1 回答 1

3

You're guess is correct; floats have a limited precision, and some of the input precision will be lost in the process. If you want the exact input, you'll need some sort of arbitrary-precision floating-point library (for example, the GNU MP library).

于 2013-06-30T04:10:36.557 回答