计算科学 - c中的高精度离散傅里叶变换 - 吾爱随笔录

我正在尝试对信号进行高精度离散傅立叶变换。为了检查精度，我使用高斯函数作为信号，因为傅立叶变换也是高斯函数。

代码如下

#include <stdio.h>
#include <tgmath.h>

typedef long double complex cplx;
long double PI;
void _fft(cplx buf[], cplx out[], int n, int step)
{
if (step < n) {
    _fft(out, buf, n, step * 2);
    _fft(out + step, buf + step, n, step * 2);

    for (int i = 0; i < n; i += 2 * step) {
        cplx t = cexpl(- I * PI * i / n) * out[i + step];
        buf[i / 2]     = out[i] + t;
        buf[(i + n)/2] = out[i] - t;
    }
}
}

void fft(cplx buf[], int n)
{
cplx out[n];
for (int i = 0; i < n; i++) out[i] = buf[i];

_fft(buf, out, n, 1);
}


int main()
{
const int nPoints = pow(2, 12);
PI = atan2(1, 1) * 4;
cplx dt = 1e-3;
cplx dOmega = 1 / (dt * nPoints);
long double T[nPoints];
long double DOmega[nPoints];

cplx At[nPoints];
cplx tau = 28.4e-3;
for (int i = 0; i < nPoints; ++i)
{
    T[i] = dt*(i-nPoints/2);
    DOmega[i] = dOmega * (i - nPoints / 2);
    At[i] = cexpl(-T[i]*T[i]/2/(tau*tau));
}

fft(At, nPoints);
FILE* fw;
fw = fopen("fft_01.txt", "w+");

for (int i = 0; i < nPoints; ++i)
{
    fprintf(fw, "%.15Le, %.15Le\n", DOmega[i], fabs(At[i]) );
}

return 0;
}

当我将定义的类型从 float complex 更改为 double complex 时，结果确实显示了精度的提高。但是，当我将数据类型从双复数更改为长双复数时，虽然数字远离最小限制，但没有任何改进。

我认为峰值的两个翅膀应该像 Mathematica 的结果一样下降，但这里似乎卡在 10^(-35) 附近，无法进一步下降。

当我进一步将 iDFT 应用于 DFT 的结果时，更明显的是 long double 精度已丢失。蓝色的是原始信号，精度高。绿色的是iDFT(DFT(x))之后的结果，应该与蓝色的一致，但最小的数字要大得多。

谁能告诉我问题出在哪里，在此先感谢。

PS：我在 Mac OS 10.13.1 上，gcc用作编译器。我已经检查了我平台上的数字限制

storage size for float: 4
minimum float positive value: 1.175494e-38
Maximum float positive value: 3.402823e+38
precision value: 6

storage size for double: 8
minimum double positive value: 2.225074e-308
Maximum double positive value: 1.797693e+308
precision value: 15

storage size for long double: 16
minimum long double positive value: 3.362103e-4932
Maximum long double positive value: 1.189731e+4932
precision value: 18