使用 NFFT3 库的简单示例

Question

我正在尝试使用 NFFT3 库 执行从非等距空间数据到相应等距傅里叶系数的 NFFT。我正在 C 的 Linux 环境中实现我的项目。

为了测试该库的使用，我正在实现从具有值

[-0.1, 0.1]

的节点

[1.0, 1.0]

到

N=10

傅立叶系数的转换。预期结果应该是

2.0*cos(pi*nu/5)

，其中

nu

是频率。

我对 FFTW 有一些经验，但我似乎不知道如何让 NFFT3 工作。任何有关如何改进我的 C 代码以使用 NFFT3 获得预期结果的建议将非常感激。

我到目前为止的C代码如下：

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <complex.h>
#include "nfft3.h"

static void test_nfft_1d(void)
{
  /** problem size */
  int N=10; // # of Fourier coefficient sample points - equispaced
  int M=2;  // # of non-equispaced nodes

  /** init one dimensional plan */
  nfft_plan p;
  nfft_init_1d(&p, N, M);

  /** init non-equidistant nodes */
  p.x[0] = - 0.1;
  p.x[1] = 0.1;

  /** precompute psi */
  if(p.flags & PRE_ONE_PSI)
      nfft_precompute_one_psi(&p);

  /** init spacial values */
  p.f[0] = 1.0;
  p.f[1] = 1.0;

  /** perform adjoint nfft */
  nfft_adjoint(&p);

  /** show the result */
  for (int idx_k = 0; idx_k < N; idx_k++) {
    printf("%#Zg\n",p.f_hat[idx_k]);
  }

  /** finalise the one dimensional plan */
  nfft_finalize(&p);
}

int main(void)
{
  // computing a one dimensional nfft
  test_nfft_1d();

  return 1;
}

作为输出，我得到：

-2.00000
-1.61803
-0.618034
0.618034
1.61803
2.00000
1.61803
0.618034
-0.618034
-1.61803

我期望到达的地方：

2.00000     
1.53208 
0.34729     
-0.99999    
-1.87938    
-1.87938    
-1.00000    
0.34729 
1.53208     
1.99999 

以下 Matlab 代码产生预期结果：

%% Define Signal Parameters
N = 10; % # of Fourier coefficient sample points - equispaced
nu_range = 40; 
nu = linspace(-nu_range,nu_range,N);

%% Define Non-Equidistant Spacial Signal
x_nu = [-0.1,0.1]; % non-equidistant nodes
sig_nu = [1.0,1.0]; % signal at given nodes
 
%% Compute the NFFT
NFFT_mat = nufft(double(sig_nu),double(x_nu),kx);
NFFT_an = ndft(double(sig_nu),double(x_nu),kx);
FFT_an = 2.0*cos(pi*kx/5);

% NFFT results
figure(1)
hold on;
plot(kx,real(NFFT_mat),"--");
plot(kx,real(NFFT_an),":");
plot(kx,real(FFT_an),".-.");
legend("Matlab NFFT","Analytic NFFT","Rect function FT");
xlabel("frequency");
ylabel("real part of Fourier coefficient");
title("NFFT Results [Real Part]")
hold off;

%% Analytical expression for NDFT
function [out] = ndft(X,t,f)
  out = zeros(size(f));
  for j = 1:length(X)
    out(:) = out(:) + X(j) * exp(-1j*2*pi*t(j)*f(:));
  end
end

Answer 1

每个 NUFFT 代码都有不同的定义，因此与您自己编写的直接变换（作为您的 ndft 函数）进行比较至关重要。 MathWorks nufft 有一个非常规的定义，而你的 matlab 代码没有定义 kx，所以不完整，所以我无法提供帮助。我建议我们的 FINUFFT 库可能比 NFFT3 更容易使用。它具有不依赖于其自身实用函数的示例代码。祝你好运，亚历克斯