一维非均匀FFT正反向测试
一维非均匀FFT正反向测试
提问于 2020-02-03 08:35:06
我正在学习如何使用c++库来执行非均匀FFT (NUFFT)。该库提供3种类型的NUFFT。
- 类型1:从非均匀x网格到均匀k空间网格的正向转换。
- 类型2:从均匀k空间网格向非均匀x网格的反向转换
- 第三种类型:从非制服到非均匀
通过对测试函数sin(x)从-pi执行NUFFT到使用Type1 NUFFT对pi进行测试，在1D中测试库，使用Type2 NUFFT将其转换回来，并将输出与sin(x)进行比较。首先，我在一个统一的x网格上测试了它，这显示了一个非常小的错误。不幸的是，当测试是在非均匀的x网格上进行时，错误是非常大的。
有两种可能性：
- 我的NUFFT实现是不正确的，但实现相当简单，所以我怀疑是否是这样。
- 作者提到Type2不是Type1的反面，所以我认为这可能是问题所在。由于我不是NUFFT方面的专家，我想知道是否有一种替代的方法来执行NUFFT的前后向测试？
- 我的目的是在不规则网格上开发一个FFT泊松解算器，因此我需要向前和向后执行NUFFT，因此克服这个问题非常重要。除了使用FINUFFT，任何其他建议也是受欢迎的。
感谢您的阅读。这是给那些有兴趣的人的密码。
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex>
#include <fftw3.h>
#include <functional>
#include "finufft/src/finufft.h"
using namespace std;
int main()
{
double pi = 3.14159265359;
int N = 128*2;
int i;
double X[N];
double L = 2*pi;
double dx = L/(N);
nufft_opts opts; finufft_default_opts(&opts);
complex<double> R = complex<double>(1.0,0.0); // the real unit
complex<double> in1[N], out1[N], out2[N];
for(i = 0; i < N; i++) {
//X[i] = -(L/2) + i*dx ; // uniform grid
X[i] = -(L/2) + pow(double(i)/N,7.0)*L; //non-uniform grid
in1[i] = sin(X[i])*R ;}
int ier = finufft1d1(N,X,in1,-1,1e-10,N,out1,opts); // type-1 NUFFT
int ier2 = finufft1d2(N,X,out2,+1,1e-10,N,out1,opts); // type-2 NUFFT
// checking the error
double erl1 = 0.;
for ( i = 0; i < N; i++) {
erl1 += fabs( in1[i].real() - out2[i].real()/(N))*dx;
}
std::cout<< erl1 <<" " << ier << " "<< ier2<< std::endl ; // error
return 0;
}回答 1
Stack Overflow用户
发布于 2020-03-19 06:52:45
出于某种原因，开发人员在他们的页面上做了一个更新，准确地回答了我的问题。https://finufft.readthedocs.io/en/latest/examples.html#periodic-poisson-solve-on-non-cartesian-quadrature-grid。简单地说，他们的NUFFT代码在完全自适应的情况下是不好的，但是我仍然会在这里提供一个完整的答案和代码。
我的代码中缺少两种成分。
(1)在使用NUFFT之前，我需要用权重将函数sin(x)乘以。权重来自于二维例子中雅可比行列式的行列式，因此对于一个一维例子来说，权值只是非均匀坐标的导数。
(2) Nk小于N。
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex>
#include <fftw3.h>
#include <functional>
#include "finufft/src/finufft.h"
using namespace std;
int main()
{
double pi = 3.14159265359;
int N = 128*2;
int Nk = 32; // smaller than N
int i;
double X[N];
double L = 2*pi;
double dx = L/(N);
nufft_opts opts; finufft_default_opts(&opts);
complex<double> R = complex<double>(1.0,0.0); // the real unit
complex<double> in1[N], out1[N], out2[N];
for(i = 0; i < N; i++) {
ksi[i] = -(L/2) + i*dx ; //uniform grid
X[i] = -(L/2) + pow(double(i)/(N-1),6)*L; //nonuniform grid
}
dX = der(ksi,X,1); // your own derivative code
for(i = 0; i < N; i++) {
in1[i] = sin(X[i]) * dX[i] * R ; // applying weight
}
int ier = finufft1d1(N,X,in1,-1,1e-10,Nk,out1,opts); // type-1 NUFFT
int ier2 = finufft1d2(N,X,out2,+1,1e-10,Nk,out1,opts); // type-2 NUFFT
// checking the error
double erl1 = 0.;
for ( i = 0; i < N; i++) {
erl1 += fabs( in1[i].real() - out2[i].real()/(N))*dx;
}
std::cout<< erl1 <<" " << ier << " "<< ier2<< std::endl ; // error
return 0;
}页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/60035611
复制相关文章
相似问题
腾讯云开发者
