nfft.c

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

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

#define MACRO_nfft_B(which_one)                                                \
inline void nfft_B_ ## which_one (nfft_plan *this)                             \
{                                                                              \
  int lprod;                            /**< 'regular bandwidth' of matrix B */\
  int u[this->d], o[this->d];           /**< multi band with respect to x_j  */\
  int t, t2;                            /**< index dimensions                */\
  int j;                                /**< index nodes                     */\
  int l_L, ix;                          /**< index one row of B              */\
  int l[this->d];                       /**< multi index u<=l<=o             */\
  int lj[this->d];                      /**< multi index 0<=lj<u+o+1         */\
  int ll_plain[this->d+1];              /**< postfix plain index in g        */\
  double phi_prod[this->d+1];           /**< postfix product of PHI          */\
  fftw_complex *f, *g;                  /**< local copy                      */\
  double *fj0,*fj1;                     /**< local copy                      */\
                                                                               \
  f=this->f; g=this->g;                                                        \
                                                                               \
  MACRO_nfft_B_init_result_ ## which_one;                                      \
                                                                               \
  if((this->nfft_flags & PRE_PSI)&&(this->nfft_flags & PRE_FULL_PSI))          \
    {                                                                          \
      for(ix=0, j=0, fj0=&f[0][0], fj1=&f[0][1]; j<this->M; j++,fj0+=2,fj1+=2) \
        for(l_L=0; l_L<this->psi_index_f[j]; l_L++, ix++)                      \
	  MACRO_nfft_B_PRE_FULL_PSI_compute_ ## which_one;                     \
      return;                                                                  \
    }                                                                          \
                                                                               \
  phi_prod[0]=1;                                                               \
  ll_plain[0]=0;                                                               \
                                                                               \
  for(t=0,lprod = 1; t<this->d; t++)                                           \
    lprod *= (2*this->m+2);                                                    \
                                                                               \
  if(this->nfft_flags & PRE_PSI)                                               \
    {                                                                          \
      for(j=0, fj0=&f[0][0], fj1=&f[0][1]; j<this->M; j++, fj0+=2, fj1+=2)     \
	{                                                                      \
          MACRO_init_uo_l_lj_t;                                                \
                                                                               \
	  for(l_L=0; l_L<lprod; l_L++)                                         \
	    {                                                                  \
              MACRO_update_phi_prod_ll_plain(with_PRE_PSI);                    \
                                                                               \
	      MACRO_nfft_B_compute_ ## which_one;                              \
		                                                               \
	      MACRO_count_uo_l_lj_t;                                           \
            } /* for(l_L) */                                                   \
	} /* for(j) */                                                         \
      return;                                                                  \
    } /* if(PRE_PSI) */                                                        \
                                                                               \
  /* no precomputed psi at all */                                              \
  for(j=0, fj0=&f[0][0], fj1=&f[0][1]; j<this->M; j++, fj0+=2, fj1+=2)         \
    {                                                                          \
      MACRO_init_uo_l_lj_t;                                                    \
	                                                                       \
      for(l_L=0; l_L<lprod; l_L++)                                             \
     	{                                                                      \
          MACRO_update_phi_prod_ll_plain(without_PRE_PSI);                     \
                                                                               \
          MACRO_nfft_B_compute_ ## which_one;                                  \
		                                                               \
          MACRO_count_uo_l_lj_t;                                               \
	} /* for(l_L) */                                                       \
    } /* for(j) */                                                             \
} /* nfft_B */                                                                 \

MACRO_nfft_B(A)
MACRO_nfft_B(T)


/** user routines
 */
void nfft_trafo(nfft_plan *this)
{
  /* use this->my_fftw_plan1 */
  this->g_hat=this->g1;
  this->g=this->g2;
 
  /** form \f$ \hat g_k = \frac{\hat f_k}{c_k\left(\phi\right) \text{ for }
   *  k \in I_N \f$
   */ 
  T1;
  nfft_D_A(this);
  T2(1);

  /** compute by d-variate discrete Fourier transform
   *  \f$ g_l = \sum_{k \in I_N} \hat g_k {\rm e}^{-2\pi {\rm i} \frac{kl}{n}}
   *  \text{ for } l \in I_n \f$
   */
  T1;
  fftw_execute(this->my_fftw_plan1);
  T2(2);

  /** set \f$ f_j = \sum_{l \in I_n,m(x_j)} g_l \psi\left(x_j-\frac{l}{n}\right)
   *  \text{ for } j=0,\hdots,M-1 \f$
   */
  T1;
  nfft_B_A(this);
  T2(3);
} /* nfft_trafo */

void nfft_conjugated(nfft_plan *this)
{
  /* use this->my_fftw_plan2 */
  this->g_hat=this->g2;
  this->g=this->g1;
 
  /** form \f$ \hat g_k = \frac{\hat f_k}{c_k\left(\phi\right) \text{ for }
   *  k \in I_N \f$
   */  
  T1;
  nfft_D_A(this);
  T2(1);

  /** compute by d-variate discrete Fourier transform
   *  \f$ g_l = \sum_{k \in I_N} \hat g_k {\rm e}^{+2\pi {\rm i} \frac{kl}{n}}
   *  \text{ for } l \in I_n \f$
   */
  T1;
  fftw_execute(this->my_fftw_plan2);
  T2(2);

  /** set \f$ f_j = \sum_{l \in I_n,m(x_j)} g_l \psi\left(x_j-\frac{l}{n}\right)
   *  \text{ for } j=0,\hdots,M-1 \f$
   */
  T1;
  nfft_B_A(this);
  T2(3);
} /* nfft_conjugated */

void nfft_adjoint(nfft_plan *this)
{
  /* use this->my_fftw_plan2 */
  this->g_hat=this->g1;
  this->g=this->g2;
  
  /** set \f$ g_l = \sum_{j=0}^{M-1} f_j \psi\left(x_j-\frac{l}{n}\right)
   *  \text{ for } l \in I_n,m(x_j) \f$
   */
  T1;
  nfft_B_T(this);
  T2(1);
 
  /** compute by d-variate discrete Fourier transform
   *  \f$ \hat g_k = \sum_{l \in I_n} g_l {\rm e}^{+2\pi {\rm i} \frac{kl}{n}}
   *  \text{ for }  k \in I_N\f$
   */
  T1;
  fftw_execute(this->my_fftw_plan2);
  T2(2);
  
  /** form \f$ \hat f_k = \frac{\hat g_k}{c_k\left(\phi\right) \text{ for }
   *  k \in I_N \f$
   */
  T1;
  nfft_D_T(this);
  T2(3);
} /* nfft_adjoint */

void nfft_transposed(nfft_plan *this)
{
  /* use this->my_fftw_plan1 */
  this->g_hat=this->g2;
  this->g=this->g1;

  /** set \f$ g_l = \sum_{j=0}^{M-1} f_j \psi\left(x_j-\frac{l}{n}\right)
   *  \text{ for } l \in I_n,m(x_j) \f$
   */
  T1;
  nfft_B_T(this);
  T2(1);
 
  /** compute by d-variate discrete Fourier transform
   *  \f$ \hat g_k = \sum_{l \in I_n} g_l {\rm e}^{-2\pi {\rm i} \frac{kl}{n}}
   *  \text{ for }  k \in I_N\f$
   */ 
  T1;
  fftw_execute(this->my_fftw_plan1);
  T2(2);
  
  /** form \f$ \hat f_k = \frac{\hat g_k}{c_k\left(\phi\right) \text{ for }
   *  k \in I_N \f$
   */
  T1;
  nfft_D_T(this);
  T2(3);
} /* nfft_transposed */


/** initialisation of direct transform 
 */
void nfft_precompute_phi_hut(nfft_plan *this)
{
  int ks[this->d];                      /**< index over all frequencies       */
  int t;                                /**< index over all dimensions        */

  this->c_phi_inv = (double**) fftw_malloc(this->d*sizeof(double*));

  for(t=0; t<this->d; t++)
    {
      this->c_phi_inv[t]= (double*)fftw_malloc(this->N[t]*sizeof(double));
      for(ks[t]=0; ks[t]<this->N[t]; ks[t]++)  
	this->c_phi_inv[t][ks[t]]= 1.0/(PHI_HUT(ks[t]-this->N[t]/2,t));
    }
} /* nfft_phi_hut */

void nfft_precompute_psi(nfft_plan *this)
{
  int t;                                /**< index over all dimensions        */
  int j;                                /**< index over all nodes             */
  int l;                                /**< index u<=l<=o                    */
  int lj;                               /**< index 0<=lj<u+o+1                */
  int u, o;                             /**< depends on x_j                   */
  
  for (t=0; t<this->d; t++)
    for(j=0;j<this->M;j++)
      {
	nfft_uo(this,j,&u,&o,t);
	
	for(l=u, lj=0; l <= o; l++, lj++)
	  this->psi[(j*this->d+t)*(2*this->m+2)+lj]=
	    (PHI((this->x[j*this->d+t]-((double)l)/this->n[t]),t));
      } /* for(j) */
  /* for(t) */
} /* nfft_precompute_psi */

/** more memory usage, a bit faster */
void nfft_full_psi(nfft_plan *this, double eps)
{
  int t,t2;                             /**< index over all dimensions        */
  int j;                                /**< index over all nodes             */
  int l_L;                              /**< plain index 0<=l_L<lprod         */
  int l[this->d];                       /**< multi index u<=l<=o              */
  int lj[this->d];                      /**< multi index 0<=lj<u+o+1          */
  int ll_plain[this->d+1];              /**< postfix plain index              */
  int lprod;                            /**< 'bandwidth' of matrix B          */
  int u[this->d], o[this->d];           /**< depends on x_j                   */
  
  double phi_prod[this->d+1];

  int *index_g, *index_f;
  double *new_psi;
  int ix,ix_old,size_psi;

  phi_prod[0]=1;
  ll_plain[0]=0;

  if(this->nfft_flags & PRE_PSI)
    {
      size_psi=this->M;
      index_f = (int*) malloc(this->M*sizeof(int));
      index_g = (int*) malloc(size_psi*sizeof(int));
      new_psi = (double*) malloc(size_psi*sizeof(double));
      
      for(t=0,lprod = 1; t<this->d; t++)
	{
	  lprod *= 2*this->m+2;
	  eps=eps*(PHI(0,t));
	}
      
      for(ix=0,ix_old=0,j=0; j<this->M; j++)
	{
	  MACRO_init_uo_l_lj_t;
	  
	  for(l_L=0; l_L<lprod; l_L++)
            {
	      MACRO_update_phi_prod_ll_plain(with_PRE_PSI);

              if(phi_prod[this->d]>eps)
		{
		  index_g[ix]=ll_plain[this->d];
		  new_psi[ix]=phi_prod[this->d];
		  ix++;
		  if(ix==size_psi)
		    {
		      size_psi+=this->M;
		      index_g=(int*)realloc(index_g,size_psi*sizeof(int));
		      new_psi=(double*)realloc(new_psi,size_psi*sizeof(double));
		    }
		}
	      MACRO_count_uo_l_lj_t;
            } /* for(l_L) */

	  
	  index_f[j]=ix-ix_old;
	  ix_old=ix;
	} /* for(j) */

      free(this->psi);
      size_psi=ix;
      this->size_psi=size_psi;
      index_g=(int*)realloc(index_g,size_psi*sizeof(int));
      new_psi=(double*)realloc(new_psi,size_psi*sizeof(double));

      this->psi         =new_psi;
      this->psi_index_g =index_g; 
      this->psi_index_f=index_f;
      
      //printf("size_psi / (lprod*M)=%e\n",size_psi / ((double)lprod*this->M));
    } /* if(PRE_PSI) */
}


void nfft_init_help(nfft_plan *this)
{
  int t;                                /**< index over all dimensions        */

  this->N_L=nfft_prod_int(this->N, this->d);
  this->n_L=nfft_prod_int(this->n, this->d);

  this->sigma = (double*) fftw_malloc(this->d*sizeof(double));
  for(t = 0;t < this->d; t++)
    this->sigma[t] = ((double)this->n[t])/this->N[t];
  
  WINDOW_HELP_INIT;

  if(this->nfft_flags & MALLOC_X)
    this->x = (double*)fftw_malloc(this->d*this->M*
					sizeof(double));

  if(this->nfft_flags & MALLOC_F_HAT)
    this->f_hat = (fftw_complex*)fftw_malloc(this->N_L*
						  sizeof(fftw_complex));
  if(this->nfft_flags & MALLOC_F)
    this->f=(fftw_complex*)fftw_malloc(this->M*sizeof(fftw_complex));

  if(this->nfft_flags & PRE_PHI_HUT)
    nfft_precompute_phi_hut(this);

  /* NO FFTW_MALLOC HERE */
  if(this->nfft_flags & PRE_PSI)
    this->psi = (double*) malloc(this->M*this->d*