nfft.c

用于实现非均匀采样FFT（NDFT）

/**
 * Library.
 * Includes simple and fast computation of the NDFT (direct problem) as well
 * as (iterative) solution to the inverse problem.
 * authors: D. Potts, S. Kunis (c) 2002-2005
 */

#include "nfft.h"
#include "window_defines.h"

/** direct computation of non equispaced fourier transforms
 *  ndft_trafo, ndft_conjugated, ndft_adjoint, ndft_transposed
 *  require O(M N^d) arithemtical operations
 *
 * direct computation of the ndft_trafo and ndft_conjugated, formula (1.1)
 * ndft_trafo:
 * for j=0,...,M-1                                                             
 *  f[j] = sum_{k in I_N^d} f_hat[k] * exp(-2 (pi) k x[j])
 * ndft_conjugated:
 * for j=0,...,M-1                                                             
 *  f[j] = sum_{k in I_N^d} f_hat[k] * exp(+2 (pi) k x[j])
 *
 * direct computation of the ndft_adjoint and ndft_transposed, formula (1.2)
 * ndft_adjoint:
 * for k in I_N^d
 *  f_hat[k] = sum_{j=0}^{M-1} f[j] * exp(+2(pi) k x[j])
 * ndft_transposed:
 * for k in I_N^d
 *  f_hat[k] = sum_{j=0}^{M-1} f[j] * exp(-2(pi) k x[j])
 */

/** macros and small sub routines for the direct transforms
 */
#define MACRO_ndft_init_result_trafo memset(f,0,this->M*sizeof(fftw_complex));
#define MACRO_ndft_init_result_conjugated MACRO_ndft_init_result_trafo
#define MACRO_ndft_init_result_adjoint memset(f_hat,0,this->N_L*               \
                                              sizeof(fftw_complex));
#define MACRO_ndft_init_result_transposed MACRO_ndft_init_result_adjoint

#define MACRO_ndft_sign_trafo      +2*PI*this->x[j*this->d+t]
#define MACRO_ndft_sign_conjugated -2*PI*this->x[j*this->d+t]
#define MACRO_ndft_sign_adjoint    +2*PI*this->x[j*this->d+t]
#define MACRO_ndft_sign_transposed -2*PI*this->x[j*this->d+t]

#define MACRO_init_k_N_Omega_x(which_one) {                                    \
for(t=0; t<this->d; t++)                                                       \
  {                                                                            \
    k[t]=-this->N[t]/2;                                                        \
    x[t]= MACRO_ndft_sign_ ## which_one;                                       \
    Omega[t+1]=k[t]*x[t]+Omega[t];                                             \
  }                                                                            \
omega=Omega[this->d];                                                          \
}                                                                              \

#define MACRO_count_k_N_Omega {                                                \
for(t = this->d-1; (t >= 1) && (k[t] == this->N[t]/2-1); t--)                  \
  k[t]-= this->N[t]-1;                                                         \
                                                                               \
k[t]++;                                                                        \
                                                                               \
for(t2 = t; t2<this->d; t2++)                                                  \
  Omega[t2+1]=k[t2]*x[t2]+Omega[t2];                                           \
                                                                               \
omega=Omega[this->d];                                                          \
}

#define MACRO_ndft_compute_trafo {                                             \
  (*fj0)+=(*f_hat_k0)*cos_omega+(*f_hat_k1)*sin_omega;                         \
  (*fj1)+=(*f_hat_k1)*cos_omega-(*f_hat_k0)*sin_omega;                         \
}
#define MACRO_ndft_compute_conjugated MACRO_ndft_compute_trafo

#define MACRO_ndft_compute_adjoint {                                           \
  (*f_hat_k0)+= (*fj0)*cos_omega-(*fj1)*sin_omega;                             \
  (*f_hat_k1)+= (*fj0)*sin_omega+(*fj1)*cos_omega;                             \
}
#define MACRO_ndft_compute_transposed MACRO_ndft_compute_adjoint

#define MACRO_ndft(which_one)                                                  \
void ndft_ ## which_one (nfft_plan *this)                                      \
{                                                                              \
  int j;                                /**< index over all nodes            */\
  int t,t2;                             /**< index for dimensions            */\
  int k_L;                              /**< plain index for summation       */\
  fftw_complex *f_hat, *f;              /**< dito                            */\
  double *f_hat_k0, *f_hat_k1;          /**< actual Fourier coefficient      */\
  double *fj0,*fj1;                     /**< actual sample                   */\
  double x[this->d];                    /**< actual node x[d*j+t]            */\
  int k[this->d];                       /**< multi index for summation       */\
  double omega, Omega[this->d+1];       /**< sign times 2*pi*k*x             */\
  double sin_omega,cos_omega;           /**< just for better reading         */\
                                                                               \
  f_hat=this->f_hat; f=this->f;                                                \
                                                                               \
  MACRO_ndft_init_result_ ## which_one                                         \
                                                                               \
  if(this->d==1)                                                               \
    {                                                                          \
  /* univariate case (due to performance) */                                   \
      t=0;                                                                     \
      for(j=0, fj0= &f[0][0], fj1= &f[0][1];                                   \
          j<this->M; j++, fj0+=2, fj1+=2)                                      \
        {                                                                      \
	  for(k_L=0, f_hat_k0= &f_hat[0][0], f_hat_k1= &f_hat[0][1];           \
              k_L<this->N_L; k_L++, f_hat_k0+=2, f_hat_k1+=2)                  \
	    {                                                                  \
	      omega=(k_L-this->N_L/2)* MACRO_ndft_sign_ ## which_one;          \
	      cos_omega= cos(omega);                                           \
	      sin_omega= sin(omega);                                           \
              MACRO_ndft_compute_ ## which_one;                                \
	    }                                                                  \
        }                                                                      \
    }                                                                          \
  else                                                                         \
    {                                                                          \
  /* multivariate case */                                                      \
      Omega[0]=0;                                                              \
      for(j=0, fj0=&f[0][0], fj1=&f[0][1];                                     \
          j<this->M; j++, fj0+=2, fj1+=2)                                      \
        {                                                                      \
          MACRO_init_k_N_Omega_x(which_one);                                   \
          for(k_L=0, f_hat_k0= &f_hat[0][0], f_hat_k1= &f_hat[0][1];           \
              k_L<this->N_L; k_L++, f_hat_k0+=2, f_hat_k1+=2)                  \
	    {                                                                  \
              cos_omega= cos(omega);                                           \
              sin_omega= sin(omega);                                           \
                                                                               \
              MACRO_ndft_compute_ ## which_one;                                \
                                                                               \
	      MACRO_count_k_N_Omega;                                           \
	    } /* for(k_L) */                                                   \
        } /* for(j) */                                                         \
    } /* else */                                                               \
} /* ndft_trafo */


/** user routines
 */
MACRO_ndft(trafo)
MACRO_ndft(conjugated)
MACRO_ndft(adjoint)
MACRO_ndft(transposed)


/** fast computation of non equispaced fourier transforms
 *  require O(N^d log(N) + M) arithemtical operations
 *
 * fast computation of the nfft_trafo and nfft_conjugated, formula (1.1)
 * nfft_trafo:
 * for j=0,...,M-1                                                             
 *  f[j] = sum_{k in I_N^d} f_hat[k] * exp(-2 (pi) k x[j])
 * nfft_conjugated:
 * for j=0,...,M-1                                                             
 *  f[j] = sum_{k in I_N^d} f_hat[k] * exp(+2 (pi) k x[j])
 *
 * direct computation of the nfft_adjoint and nfft_transposed, formula (1.2)
 * nfft_adjoint:
 * for k in I_N^d
 *  f_hat[k] = sum_{j=0}^{M-1} f[j] * exp(+2(pi) k x[j])
 * nfft_transposed:
 * for k in I_N^d
 *  f_hat[k] = sum_{j=0}^{M-1} f[j] * exp(-2(pi) k x[j])
 */

/** macros and small sub routines for the fast transforms
 */

/** computes 2m+2 indices for the matrix B
 */
void nfft_uo(nfft_plan *this,int j,int *up,int *op,int act_dim)
{
  double c;
  int u,o;

  c = this->x[j*this->d+act_dim] * this->n[act_dim];
  u = c; o = c;
  if(c < 0)                  
    u = u-1;                  
  else
    o = o+1;
  
  u = u - (this->m); o = o + (this->m);

  up[0]=u; op[0]=o;
}

#define MACRO_nfft_D_compute_A {                                               \
 g_hat[k_plain[this->d]][0]= f_hat[ks_plain[this->d]][0]*c_phi_inv_k[this->d]; \
 g_hat[k_plain[this->d]][1]= f_hat[ks_plain[this->d]][1]*c_phi_inv_k[this->d]; \
}

#define MACRO_nfft_D_compute_T {                                               \
 f_hat[ks_plain[this->d]][0]= g_hat[k_plain[this->d]][0]*c_phi_inv_k[this->d]; \
 f_hat[ks_plain[this->d]][1]= g_hat[k_plain[this->d]][1]*c_phi_inv_k[this->d]; \
}

#define MACRO_nfft_D_init_result_A  memset(g_hat,0,this->n_L*                  \
                                          sizeof(fftw_complex));
#define MACRO_nfft_D_init_result_T memset(f_hat,0,this->N_L*                   \
                                          sizeof(fftw_complex));

#define MACRO_with_PRE_PHI_HUT * this->c_phi_inv[t2][ks[t2]];
#define MACRO_without_PRE_PHI_HUT / (PHI_HUT(ks[t2]-this->N[t2]/2,t2));

#define MACRO_init_k_ks {                                                      \
  for(t = this->d-1; t>=0; t--)                                                \
    {                                                                          \
      kp[t]= 0;                                                                \
      k[t] = 0;                                                                \
      ks[t] = this->N[t]/2;                                                    \
    }                                                                          \
  t++;                                                                         \
}

#define MACRO_update_c_phi_inv_k(which_one) {                                  \
  for(t2=t; t2<this->d; t2++)                                                  \
    {                                                                          \
      c_phi_inv_k[t2+1]= c_phi_inv_k[t2] MACRO_ ##which_one;                   \
      ks_plain[t2+1]= ks_plain[t2]*this->N[t2]+ks[t2];                         \
      k_plain[t2+1]= k_plain[t2]*this->n[t2]+k[t2];                            \
    }                                                                          \
}

#define MACRO_count_k_ks {                                                     \
  for(t=this->d-1; (t>0)&& (kp[t]==this->N[t]-1); t--)                         \
    {                                                                          \
      kp[t]= 0;                                                                \
      k[t]= 0;                                                                 \
      ks[t]= this->N[t]/2;                                                     \
    }                                                                          \
                                                                               \
  kp[t]++; k[t]++; ks[t]++;                                                    \
  if(kp[t]==this->N[t]/2)                                                      \
    {                                                                          \
      k[t]= this->n[t]-this->N[t]/2;                                           \
      ks[t]= 0;                                                                \
    }                                                                          \
}                                                                              \


/** sub routines for the fast transforms
 *  matrix vector multiplication with \f$D, D^T\f$
 */
#define MACRO_nfft_D(which_one)                                                \
inline void nfft_D_ ## which_one (nfft_plan *this)                             \
{                                                                              \
  int t, t2;                            /**< index dimensions                */\
  int k_L;                              /**< plain index                     */\
  int kp[this->d];                      /**< multi index (simple)            */\
  int k[this->d];                       /**< multi index in g_hat            */\
  int ks[this->d];                      /**< multi index in f_hat, c_phi_inv */\
  double c_phi_inv_k[this->d+1];        /**< postfix product of PHI_HUT      */\
  int k_plain[this->d+1];               /**< postfix plain index             */\
  int ks_plain[this->d+1];              /**< postfix plain index             */\
  fftw_complex *f_hat, *g_hat;          /**< local copy                      */\
                                                                               \
  f_hat=this->f_hat; g_hat=this->g_hat;                                        \
  MACRO_nfft_D_init_result_ ## which_one;                                      \
                                                                               \
  c_phi_inv_k[0]=1;                                                            \
  k_plain[0]=0;                                                                \
  ks_plain[0]=0;                                                               \
                                                                               \
  if(this->nfft_flags & PRE_PHI_HUT)                                           \
    {                                                                          \
      MACRO_init_k_ks;                                                         \
                                                                               \
      for(k_L=0; k_L<this->N_L; k_L++)                                         \
	{                                                                      \
          MACRO_update_c_phi_inv_k(with_PRE_PHI_HUT);                          \
                                                                               \
	  MACRO_nfft_D_compute_ ## which_one;                                  \
	                                                                       \
	  MACRO_count_k_ks;                                                    \
	} /* for(k_L) */                                                       \
    } /* if(PRE_PHI_HUT) */                                                    \
  else                                                                         \
    {                                                                          \
      MACRO_init_k_ks;                                                         \
                                                                               \
      for(k_L=0; k_L<this->N_L; k_L++)                                         \
	{                                                                      \
          MACRO_update_c_phi_inv_k(without_PRE_PHI_HUT);                       \
                                                                               \
	  MACRO_nfft_D_compute_ ## which_one;                                  \
	                                                                       \
	  MACRO_count_k_ks;                                                    \
	} /* for(k_L) */                                                       \
    } /* else(PRE_PHI_HUT) */                                                  \
} /* nfft_D */

MACRO_nfft_D(A)
MACRO_nfft_D(T)

/** sub routines for the fast transforms
 *  matrix vector multiplication with \f$B, B^{\rm T}\f$
 */ 
#define MACRO_nfft_B_init_result_A  memset(f,0,this->M*sizeof(fftw_complex));
#define MACRO_nfft_B_init_result_T memset(g,0,this->n_L*sizeof(fftw_complex));

#define MACRO_nfft_B_PRE_FULL_PSI_compute_A {                                  \
  (*fj0)+=this->psi[ix]*g[this->psi_index_g[ix]][0];                           \
  (*fj1)+=this->psi[ix]*g[this->psi_index_g[ix]][1];                           \
}

#define MACRO_nfft_B_PRE_FULL_PSI_compute_T {                                  \
  g[this->psi_index_g[ix]][0]+=this->psi[ix]*(*fj0);                           \
  g[this->psi_index_g[ix]][1]+=this->psi[ix]*(*fj1);                           \
}

#define MACRO_nfft_B_compute_A {                                               \
  (*fj0) += phi_prod[this->d]* g[ll_plain[this->d]][0];                        \
  (*fj1) += phi_prod[this->d]* g[ll_plain[this->d]][1];                        \
}

#define MACRO_nfft_B_compute_T {                                               \
  g[ll_plain[this->d]][0] += phi_prod[this->d]* (*fj0);                        \
  g[ll_plain[this->d]][1] += phi_prod[this->d]* (*fj1);                        \
}

#define MACRO_with_PRE_PSI     this->psi[(j*this->d+t2)*(2*this->m+2)+lj[t2]]
#define MACRO_without_PRE_PSI  PHI(this->x[j*this->d+t2]-                      \
                               ((double)l[t2])/this->n[t2], t2)

#define MACRO_init_uo_l_lj_t {                                                 \
  for(t = this->d-1; t>=0; t--)                                                \
    {                                                                          \
      nfft_uo(this,j,&u[t],&o[t],t);                                           \
      l[t] = u[t];                                                             \
      lj[t] = 0;                                                               \
    } /* for(t) */                                                             \
  t++;                                                                         \
}

#define MACRO_update_phi_prod_ll_plain(which_one) {                            \
  for(t2=t; t2<this->d; t2++)                                                  \
    {                                                                          \
      phi_prod[t2+1]=phi_prod[t2]* MACRO_ ## which_one;                        \
      ll_plain[t2+1]=ll_plain[t2]*this->n[t2] +(l[t2]+this->n[t2])%this->n[t2];\
    } /* for(t2) */                                                            \
}

#define MACRO_count_uo_l_lj_t {                                                \