我已经使用前向 FFT 和 IFFT(对结果进行归一化)测试了我的代码的一些实际功能,这工作正常。
不过,我想对实函数求二阶导数。为了简单起见,我将 sin(2*pi*t)
作为测试用例。这是我使用的相关代码(库中的 FFT 函数):
int main(void)
{
int i;
int nyh = (N/2) + 1;
double result_array[nyh][2];
double x_k[nyh][2];
double x_r[N];
FILE* psit;
psit=fopen("psitest.txt","w");
init();
fft(x, result_array); //function in a library, this has been tested
psi(result_array, x_k);
ifft(x_k, x_r); //function in a library, this has been tested
for(i=0;i<N;i++)
{
fprintf(psit, "%f\n", x_r[i]);
}
fclose(psit);
return 0;
}
void psi(double array[nyh][2], double out[nyh][2])
{
int i;
for ( i = 0; i < N/2; i++ )
{
out[i][0] = -4.0*pi*pi*i*i*array[i][0];
out[i][1] = -4.0*pi*pi*i*i*array[i][1];
}
out[N/2][0]=0.0;
out[N/2][1]=0.0;
}
void init()
{
int i;
for(i=0;i<N;i++)
{
x[i] = sin(2.0*pi*i/N);
}
}
现在问题来了:该算法对 sin( 2*pi*t*K)
形式的任何函数都能完美运行,其中 K 是一个整数,但如果我将其作为测试function sin(3*pi*t)
,算法失败。我看不到编码中的错误。
请注意,因为函数是实数,所以我只需要取k
值的一半。这不是问题。
最佳答案
我的猜测是 sin(3*pi*t)
引入了不连续性,因为它不会在您的样本间隔中给出整数个周期。对于大多数与 FFT 相关的应用程序,您将应用一个窗口函数来处理此类不连续性,但显然这会在您的导数中引入一个误差项,我不确定您是否能够纠正这一点。
关于c - 在 FFTW3 中进行二阶导数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/6451997/