使用FFT反复计算函数导数时，会积累数值误差。

Question

我编写了一个C程序，使用FFTW计算函数的导数（重复计算）。我正在测试简单的sin(x)函数。每一步计算前一步答案的导数。我注意到误差会逐渐增大，20步后的结果是完全错误的。下面是示例输出。在特定点处的答案应该为0、+1或-1，但并不是这样。

---- out ----  data '(0) = 1.000000 -0.000000 
---- out ----  data '(1) = 0.000000 -0.000000 
---- out ----  data '(2) = -1.000000 0.000000 
---- out ----  data '(3) = -0.000000 0.000000 
---- out ----  data '(4) = 1.000000 -0.000000 
---- out ----  data '(5) = 0.000000 -0.000000 
---- out ----  data '(6) = -1.000000 0.000000 
---- out ----  data '(7) = -0.000000 0.000000 
---- out ----  data '(8) = 1.000000 -0.000000 
---- out ----  data '(9) = 0.000000 -0.000000 
---- out ----  data '(10) = -1.000000 0.000000 
---- out ----  data '(11) = -0.000000 0.000000 
---- out ----  data '(12) = 1.000000 -0.000002 
---- out ----  data '(13) = 0.000007 -0.000000 
---- out ----  data '(14) = -1.000000 0.000028 
---- out ----  data '(15) = -0.000113 0.000000 
---- out ----  data '(16) = 0.999997 -0.000444 
---- out ----  data '(17) = 0.001798 -0.000000 
---- out ----  data '(18) = -0.999969 0.007110 
---- out ----  data '(19) = -0.028621 0.000004  

我无法弄清楚为什么错误一直在增加。非常感谢您的任何建议。我将真实函数封装到复杂数据类型中，并将虚部设置为零。

以下是代码：

# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include <math.h>
# include <complex.h>
# include <fftw3.h>

int main ( int argc, char *argv[] ){
  int i;
  fftw_complex *in;
  fftw_complex *in2;
  fftw_complex *out;

  double pi = 3.14159265359;
  int nx = 8, k, t, tf = 20;
  double xi = 0, xf = 2*pi;
  double dx = (xf - xi)/((double)nx); // Step size

  complex double  *kx;

  fftw_plan plan_backward;
  fftw_plan plan_forward;

  in = fftw_malloc ( sizeof ( fftw_complex ) * nx );
  out = fftw_malloc ( sizeof ( fftw_complex ) * nx );
  in2 = fftw_malloc ( sizeof ( fftw_complex ) * nx );


  kx = malloc ( sizeof ( complex ) * nx );

  // only need it once, hence outside the loop
  for (k = 0; k < nx; k++){
     if (k < nx/2){
            kx[k] = I*2*pi*k/xf;
     } else if (k > nx/2){
           kx[k] = I*2*pi*(k-nx)/xf;
     } else if (k == nx/2){
        kx[k] = 0.0;
     }
  }

  // create plan outside the loop
  plan_forward = fftw_plan_dft_1d ( nx, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
  plan_backward = fftw_plan_dft_1d ( nx, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE );


  // initialize data
  for ( i = 0; i < nx; i++ )
  {
    in[i] = sin(i*dx) + I*0.0;  // note the complex notation.
  }
//-------------------- start time loop ---------------------------------------//

for (t = 0; t < tf; t++){
  // print input data 
  //for ( i = 0; i < nx; i++ ) { printf("initial data '(%f) = %f %f \n", i*dx, creal(in[i]), cimag(in[i]) ); }

  fftw_execute ( plan_forward );

  for ( i = 0; i < nx; i++ )
  {
    out[i] = out[i]*kx[i]; // for first order derivative
  }

  fftw_execute ( plan_backward );
  // normalize
  for ( i = 0; i < nx; i++ )
  {
    in2[i] = in2[i]/nx;
  }
    printf("---- out ----  data '(%d) = %f %f \n", t, creal(in2[0]), cimag(in2[0]) );
  // overwrite input array with this new output and loop over
  for ( i = 0; i < nx; i++ )
  {
    in[i] = in2[i];
  }
  // done with the curent loop.
}
//--------------------- end of loop ----------------------------------------//

  fftw_destroy_plan ( plan_forward );
  fftw_destroy_plan ( plan_backward );

  fftw_free ( in );
  fftw_free ( in2 );
  fftw_free ( out );

  return 0;
}

使用gcc source.c -lfftw3 -lm编译

更新：这是循环25次的M_PI输出。同样存在错误积累。

---- out ----  data '(0) = 1.000000 0.000000 
---- out ----  data '(1) = -0.000000 -0.000000 
---- out ----  data '(2) = -1.000000 -0.000000 
---- out ----  data '(3) = 0.000000 0.000000 
---- out ----  data '(4) = 1.000000 0.000000 
---- out ----  data '(5) = -0.000000 -0.000000 
---- out ----  data '(6) = -1.000000 -0.000000 
---- out ----  data '(7) = 0.000000 0.000000 
---- out ----  data '(8) = 1.000000 0.000000 
---- out ----  data '(9) = -0.000000 -0.000000 
---- out ----  data '(10) = -1.000000 -0.000000 
---- out ----  data '(11) = 0.000000 0.000000 
---- out ----  data '(12) = 1.000000 0.000000 
---- out ----  data '(13) = -0.000000 -0.000000 
---- out ----  data '(14) = -1.000000 -0.000000 
---- out ----  data '(15) = 0.000000 0.000000 
---- out ----  data '(16) = 1.000000 0.000001 
---- out ----  data '(17) = -0.000002 -0.000000 
---- out ----  data '(18) = -0.999999 -0.000008 
---- out ----  data '(19) = 0.000033 0.000004 
---- out ----  data '(20) = 0.999984 0.000132 
---- out ----  data '(21) = -0.000527 -0.000069 
---- out ----  data '(22) = -0.999735 -0.002104 
---- out ----  data '(23) = 0.008427 0.001099 
---- out ----  data '(24) = 0.995697 0.033667 

- Amit

---

你的圆周率值不够准确。你应该使用在math.h中定义的M_PI。 - FredK

是的，如果圆周率的值不够准确，信号就不再是连续周期性的，会出现杂散高频。然后，取导数（斜坡滤波器）会放大这些频率，使得信噪比显著降低。为了解决这个问题，可以将斜坡滤波器与低通滤波器结合使用，如http://www.owlnet.rice.edu/~elec539/Projects97/cult/node4.html所示。 - francis

更改为M_PI确实有所帮助，但效果不是很明显。误差确实减少了。以下是最后两次迭代使用M_PI的输出。---- 输出 ---- 数据'(18) = -0.999999 -0.000008 ---- 输出 ---- 数据'(19) = 0.000033 0.000004。如果超过20次迭代，则误差仍会继续累积。 - Amit

Answer 1

实际上，精确计算圆周率并不能解决问题，即使可以提高到20阶导数的精度。问题在于小误差会通过重复使用导数滤波器而被放大。为了限制这个问题，可以建议引入低通滤波器以及使用四倍精度。引入“条件数”概念有助于理解误差是如何被放大，并相应地设置滤波器。然而，计算20阶导数仍将是噩梦，因为计算20阶导数是病态的：即使对于最干净的实验输入，也根本不可能......

1.小误差总是存在的。

导数滤波器，也称为坡度滤波器，会比小频率更加放大高频率。最微小的误差在多次使用坡度滤波器时会被极大地放大。

通过打印频率来查看初始小误差：

printf("---- out ----  data '(%d) %d = %20g %20g \n", t,i, creal(out[i]), cimag(out[i]) );

使用 pi=3.14159265359，你得到：

---- out ----  data '(0) 0 =         -2.06712e-13                    0 
---- out ----  data '(0) 1 =          6.20699e-13                   -4 
---- out ----  data '(0) 2 =         -2.06823e-13          2.92322e-13 
---- out ----  data '(0) 3 =         -2.07053e-13          1.03695e-13 
---- out ----  data '(0) 4 =         -2.06934e-13                    0 
---- out ----  data '(0) 5 =         -2.07053e-13         -1.03695e-13 
---- out ----  data '(0) 6 =         -2.06823e-13         -2.92322e-13 
---- out ----  data '(0) 7 =          6.20699e-13                    4 

由于π的缺失数字引起的不连续性，所有频率都存在小的非空值，这些值通过取导数而被放大。

当使用pi=M_PI时，这些初始误差更小，但仍然非空。

---- out ----  data '(0) 0 =          1.14424e-17                    0 
---- out ----  data '(0) 1 =         -4.36483e-16                   -4 
---- out ----  data '(0) 2 =          1.22465e-16         -1.11022e-16 
---- out ----  data '(0) 3 =          1.91554e-16         -4.44089e-16 
---- out ----  data '(0) 4 =          2.33487e-16                    0 
---- out ----  data '(0) 5 =          1.91554e-16          4.44089e-16 
---- out ----  data '(0) 6 =          1.22465e-16          1.11022e-16 
---- out ----  data '(0) 7 =         -4.36483e-16                    4 

这些小错误与之前的错误一样被放大了，问题并没有完全解决。让我们在第一个循环中尝试将这些频率归零：

if(t==0){
     for (k = 0; k < nx; k++){
         if (k==1 || nx-k==1){
            out[k] = I*4.0;
         }else{
            out[k] =0.0;
         }
     }
}

这次在第一个循环t=0期间，唯一非空的频率是正确的频率。现在让我们看看第二个循环：

---- out ----  data '(1) 0 =                    0                    0 
---- out ----  data '(1) 1 =                   -4                    0 
---- out ----  data '(1) 2 =                    0                    0 
---- out ----  data '(1) 3 =         -4.44089e-16                    0 
---- out ----  data '(1) 4 =                    0                    0 
---- out ----  data '(1) 5 =          4.44089e-16                    0 
---- out ----  data '(1) 6 =                    0                    0 
---- out ----  data '(1) 7 =                    4                    0 

由于DFT反/正向变换和缩放过程中的有限精度计算，会出现小误差并被放大。再次出现。

2. 为了限制误差增长，可以引入滤波。

大多数实验输入都受到高频噪声的干扰，可以通过应用低通滤波器（如Butterworth滤波器）来减少噪声。有关详细信息和替代方案，请参见https://www.hindawi.com/journals/ijbi/2011/693795/。此滤波器由一个截止频率kc和一个指数组成，斜坡滤波器的频率响应将被修改如下：

    //parameters of Butterworth Filter:
    double kc=3;
    double n=16;
    // only need it once, hence outside the loop
    for (k = 0; k < nx; k++){
      if (k < nx/2){
        // add low pass filter:
        kx[k] = I*2*pi*k/xf;
        kx[k]*=1./(1.+pow(k/kc,n));
      } else if (k > nx/2){
        kx[k] = I*2*pi*(k-nx)/xf;
        kx[k]*=1./(1.+pow((nx-k)/kc,n));
      } else if (k == nx/2){
        kx[k] = 0.0;
      }
    }

使用这些参数，第20个导数的误差从5.27e-7降至1.22e-12。

另一个改进是在导数之间不返回实际空间。这样，在浮点计算期间避免了许多舍入误差。在这种特殊情况下，将输入频率清零可以确保误差保持为零，但它略带人工的味道...从实用的角度来看，如果输入信号以实际空间的形式提供，则几乎需要使用滤波器来计算导数。

3. 由于导数滤波器的条件数增加，误差正在增加

导数是一种线性应用程序，并且其特征由条件数确定。假设输入在所有频率上都存在误差eps。如果第一个频率放大了因子alpha，则频率k放大了因子k*alpha。因此，每次应用导数时，信噪比被一个比率kc（最大频率）除。如果重复20次滤波器，则信噪比将被kc^20除。

双精度数值大约精确到eps=10e-14：这是您可以获得的最佳信噪比！大多数实验输入远不如此。例如，灰度图像通常使用16位= 65536级灰度进行采样。因此，灰度图像最多只能精确到 eps=1/65536。同样，典型的音频位深为24，对应于 eps=6e-8。对于几乎解析的输入，建议使用四重精度，其精度大约为 esp=1e-34... 让我们找到频率，使得 kc^20*eps<1:

eps=10e-14    kc=5
eps=1/65536   kc=1
eps=1/2^24    kc=2
esp=1e-34     kc=44

因此，如果输入是双精度，最多只有20阶导数的4个频率会显著影响......所有高于第4个频率的频率必须通过强低通滤波器进行过滤。因此，建议使用四重精度：请参见 fftw's documentation以将fftw编译为gcc的四重精度类型__float128，并链接到quadmath。如果输入是图像，则计算20阶导数根本不在范围内：没有任何频率会显著影响！

- francis

---

我同意你的答案。谢谢!!! 最终，我想在循环中使用FFTW计算一些导数或Laplacian的术语。术语本身在循环内更改，因此与计算起始函数的20次导数不同。我想通过这里发布的代码来测试最终实现，结果证明它是一个极端情况。 - Amit

不客气！实际上，对于一阶导数和拉普拉斯算子来说，这个问题并不是很重要，尽管在实验输入的情况下经常会应用滤波。确实有宝贵可靠的基于FFT的算法，可以计算各种线性微分问题的解，其中绿函数算子已知为正弦输入。 - francis