CODEWITHSUNDEEP
cfftfftwnfft
pieterinjapan
pieterinjapan

Reputation: 13

Simpe example if using NFFT3 library

I am trying to use the NFFT3 library to perform a NFFT from non-equidistant spatial data to the corresponding equidistant Fourier coefficients. I'm implementing my project in C, in a Linux environment.

To test the use of the library I am implementing a transform from nodes [-0.1, 0.1] with values [1.0, 1.0] to N=10 Fourier coefficients. The expected result should be 2.0*cos(pi*nu/5), where nu is the frequency.

I have some experience with FFTW, but I can't seem to figure out how to get NFFT3 to work. Any advice on how to improve my C-code to get the expected results using NFFT3 would be very much appreciated.

The C-code I have so far is as follows:

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

As output I get:

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

Where I expect to get:

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

The following Matlab code produces the expected results:

%% 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

Upvotes: 1

Views: 84

Answers (1)

Alex Barnett
Alex Barnett

Reputation: 1

Each NUFFT code has a different definition, so comparison against a direct transform that you write yourself (as your ndft function) is essential. The MathWorks nufft has an unconventional definition, and your matlab code does not define kx, so is incomplete, so I can't help there. I suggest that our FINUFFT library may be easier to use than NFFT3. Is has example codes that do not depend on its own utility functions. Good luck, Alex

Upvotes: 0

Related Questions