时间戳

首页 > 解决方案 > 1D 中的非均匀 FFT 前向和后向测试

c++ - 1D 中的非均匀 FFT 前向和后向测试

问题描述

我正在学习使用 c++ 库来执行非均匀 FFT ( NUFFT )。该库提供了 3 种类型的 NUFFT。

我通过使用 Type1 NUFFT 对从 -pi 到 pi 的测试函数 sin(x) 执行 NUFFT 来测试一维库，使用 Type2 NUFFT 将其转换回来，并将输出与 sin(x) 进行比较。起初，我在一个统一的 x 网格上测试它，显示出非常小的误差。不幸的是，当在非均匀 x 网格上进行测试时，错误非常大。

两种可能：

感谢您的阅读。代码在这里，有兴趣的人可以看看。

#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;
}

标签： c++fft

解决方案

出于某种原因，开发人员在他们的页面上进行了更新，这完全回答了我的问题。https://finufft.readthedocs.io/en/latest/examples.html#periodic-poisson-solve-on-non-cartesian-quadrature-grid。简而言之，他们的 NUFFT 代码在完全自适应方案的情况下并不好，但为了完整起见，我仍然会在此处提供答案和代码。

我的代码中缺少两种成分。

(1) 在使用 NUFFT 之前，我需要将函数 sin(x) 与权重相乘。权重来自二维示例中雅可比行列式的行列式，因此对于一维示例，权重只是非均匀坐标相对于统一坐标 dx/dksi 的导数。

(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;
}

推荐阅读