【发布时间】：2018-07-05 05:56:24
【问题描述】：
我编写了一个 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 值不够准确。您应该使用 math.h 中定义的 M_PI
-
是的，确实，如果 pi 的值不够准确，信号就不再是连续周期性的，会出现杂散的高频。然后，采用导数（斜坡滤波器）放大这些频率，使信噪比显着降低。要解决此问题，请将斜坡滤波器与低通滤波器结合使用，如 owlnet.rice.edu/~elec539/Projects97/cult/node4.html 所示
-
更改为 M_PI 确实有帮助，但作用不大。错误确实减少了。这是最后两次迭代的 M_PI 输出。 ----出----数据'（18）= -0.999999 -0.000008 ----出----数据'（19）= 0.000033 0.000004。如果我超过 20，那么错误仍然会不断累积。