FFT反复计算函数导数时出现的数值误差
Numerical error building up when computing derivative of function repeatedly with FFT
我编写了一个 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
更新：这是 M_PI 循环 25 次的输出。相同的错误累积。
---- 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
确实，细化 pi 的值并不能解决您的问题，即使它提高了 20 次导数的准确性。 问题是重复微分滤波器会放大小误差。要限制此问题，建议引入低通滤波器并使用四精度。引入 条件数 的概念有助于理解错误扩大的方式并相应地设置过滤器。 然而，计算 20 阶导数仍将是一场噩梦，因为 computing the 20th derivative is ill-conditionned：这根本不可能，即使对于最干净的实验输入...
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 的缺失数字引起的不连续性，所有频率都有小的非空值，这些值通过取导数膨胀。
由于使用了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 backward/forward 变换和缩放期间的有限精度计算，会出现小错误并被放大。再次。
2。为了限制错误的增长，可以引入过滤。
大多数实验输入都受到较大的高频噪声的干扰，可以通过应用低通滤波器（例如巴特沃斯滤波器）来降低这种噪声。有关详细信息和替代方案，请参阅 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。
另一个改进是通过不返回导数之间的真实 space 来实现的。这样，就避免了浮点计算过程中的大量舍入错误。在这种特殊情况下，将输入频率归零可确保误差保持为零，但它有点人为......从实际角度来看，如果输入信号以真实 space 形式提供，则使用滤波器计算几乎需要衍生品。
3。由于导数滤波器的条件数
，误差越来越大
导数是一个线性应用，它的特点是condition number。假设输入在所有频率上都受到错误 eps 的困扰。如果第一个频率被放大 alpha 倍，则频率 k 被放大 k*alpha 倍。因此，每次应用导数时，信噪比除以比率 kc（最大频率），称为条件数。如果滤波器重复20次，则信噪比除以kc^20.
双精度数大约 eps=10e-14 精确：这是您可以获得的最佳信噪比！大多数实验输入会比这差得多。例如，灰度图像通常使用 16bit=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
因此，如果输入是双精度的，最多只有 4 个频率的 20 阶导数是重要的... 所有高于 4 的频率都必须被过滤掉通过强大的低通滤波器。因此，可以建议使用四精度：参见 fftw's documentation to compile fftw for gcc's quad precision type __float128 linking against quadmath。 如果输入是图像，计算 20 阶导数就超出了范围：none 的频率将变得非常重要！
我编写了一个 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
更新：这是 M_PI 循环 25 次的输出。相同的错误累积。
---- 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
确实，细化 pi 的值并不能解决您的问题，即使它提高了 20 次导数的准确性。 问题是重复微分滤波器会放大小误差。要限制此问题，建议引入低通滤波器并使用四精度。引入 条件数 的概念有助于理解错误扩大的方式并相应地设置过滤器。 然而，计算 20 阶导数仍将是一场噩梦，因为 computing the 20th derivative is ill-conditionned：这根本不可能，即使对于最干净的实验输入...
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 的缺失数字引起的不连续性，所有频率都有小的非空值，这些值通过取导数膨胀。
由于使用了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 backward/forward 变换和缩放期间的有限精度计算，会出现小错误并被放大。再次。
2。为了限制错误的增长，可以引入过滤。
大多数实验输入都受到较大的高频噪声的干扰，可以通过应用低通滤波器（例如巴特沃斯滤波器）来降低这种噪声。有关详细信息和替代方案，请参阅 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。
另一个改进是通过不返回导数之间的真实 space 来实现的。这样，就避免了浮点计算过程中的大量舍入错误。在这种特殊情况下，将输入频率归零可确保误差保持为零，但它有点人为......从实际角度来看，如果输入信号以真实 space 形式提供，则使用滤波器计算几乎需要衍生品。
3。由于导数滤波器的条件数
，误差越来越大导数是一个线性应用，它的特点是condition number。假设输入在所有频率上都受到错误 eps 的困扰。如果第一个频率被放大 alpha 倍，则频率 k 被放大 k*alpha 倍。因此，每次应用导数时，信噪比除以比率 kc（最大频率），称为条件数。如果滤波器重复20次，则信噪比除以kc^20.
双精度数大约 eps=10e-14 精确：这是您可以获得的最佳信噪比！大多数实验输入会比这差得多。例如，灰度图像通常使用 16bit=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
因此，如果输入是双精度的，最多只有 4 个频率的 20 阶导数是重要的... 所有高于 4 的频率都必须被过滤掉通过强大的低通滤波器。因此，可以建议使用四精度：参见 fftw's documentation to compile fftw for gcc's quad precision type __float128 linking against quadmath。 如果输入是图像，计算 20 阶导数就超出了范围：none 的频率将变得非常重要！