对 FFTW3 进行二阶导数
我已经使用前向 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 是一个整数，但如果我将 sin(3*pi*t) 作为测试函数，算法就会失败。我看不到我的编码中的错误。
请注意，因为该函数是真实的，所以我只需要取 k 值的一半。这不是问题。
I have tested my code for some real functions using a forward FFT and IFFT (normalized the result), this works fine.
I would, however, like to take a second derivative of a real function. For simplicity sake, I take sin(2*pi*t) as a test case. Here is the relevant code I use (FFT functions in a library):
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);
}
}
Now here is the problem: This algorithm works perfectly for any function of the form sin( 2*pi*t*K), where K is an integer, but if I take as a test function sin(3*pi*t), the algorithm fails. I am not able to see the mistake in my coding.
Please note that because the function is real, I only need to take half of the k values. This is not the problem.
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（2）
我的猜测是，
sin(3*pi*t)引入了不连续性，因为它没有给出采样间隔中的整数个周期。对于大多数与 FFT 相关的应用程序，您将应用窗口函数来处理此类不连续性，但显然这会在导数中引入错误项，我不确定您是否能够纠正这一点。My guess is that
sin(3*pi*t)introduces a discontinuity, since it does not give an integer number of cycles in your sample interval. For most FFT-related applications you would apply a window function to deal with such discontinuities, but obviously this will introduce an error term into your derivative and I'm not sure whether you will be able to correct for this.我不知道你是否已经解决了这个问题...但我猜主要问题是 sin(3 Pi t) 在域 [0,1](sin(0) != sin (3 *圆周率））。
FFT 无法正常工作...
I don't know if you have fixed this problem... But I guess the major problem is that sin(3 Pi t) is not periodic in the domain [0,1](sin(0) != sin (3 * Pi)).
FFT could not work properly...