utils.c

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

0x88c0,0x77ac,0xf7ac,0x3fc5,
0xcd8d,0xc057,0x7feb,0xbfd3,
0xa22a,0x9035,0xa84e,0x3fe5,
};
#endif

#ifdef MIEEE
static unsigned short A[] = {
0xbc54,0x5cb7,0x2134,0xd0ef,
0x3c83,0x3362,0x977d,0xa589,
0xbcb1,0x84eb,0x721e,0xbbb4,
0x3cde,0xe6d8,0x93f6,0x5eba,
0xbd0a,0x5022,0xc297,0xfbeb,
0x3d35,0x9b46,0x4b26,0x2627,
0xbd61,0x164c,0x62ee,0x1af0,
0x3d89,0xfe2f,0xe19b,0xd324,
0xbdb2,0xfc95,0x7a94,0x6abc,
0x3dda,0x98be,0xcc74,0x3c10,
0xbe01,0xd4fe,0x13ae,0x9556,
0x3e26,0xd903,0xa454,0xcb34,
0xbe4b,0xeaf6,0x8c0b,0x30ab,
0x3e70,0x3b76,0x9d4d,0x6435,
0xbe91,0xec63,0x8f22,0x7f8d,
0x3eb2,0xbf24,0x978c,0xf4ac,
0xbed2,0x866f,0xcba5,0x6427,
0x3ef1,0x3f58,0xbe9a,0x2859,
0xbf0e,0x2b26,0x59c4,0x1d5a,
0x3f28,0xb51b,0x7410,0x7cab,
0xbf42,0xe2fd,0x1f15,0xeb52,
0x3f5a,0xdc75,0x8a12,0x100e,
0xbf71,0xb65e,0x201a,0xa849,
0x3f85,0x9961,0xf3dd,0xe3dd,
0xbf98,0x4e9e,0xf121,0xb6f0,
0x3fa9,0x3e8a,0xcea8,0xa32d,
0xbfb8,0x4b70,0x342d,0x06ea,
0x3fc5,0xf7ac,0x77ac,0x88c0,
0xbfd3,0x7feb,0xc057,0xcd8d,
0x3fe5,0xa84e,0x9035,0xa22a
};
#endif


/* Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
 * in the inverted interval [8,infinity].
 *
 * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
 */

#ifdef UNK
static double B[] =
{
-7.23318048787475395456E-18,
-4.83050448594418207126E-18,
 4.46562142029675999901E-17,
 3.46122286769746109310E-17,
-2.82762398051658348494E-16,
-3.42548561967721913462E-16,
 1.77256013305652638360E-15,
 3.81168066935262242075E-15,
-9.55484669882830764870E-15,
-4.15056934728722208663E-14,
 1.54008621752140982691E-14,
 3.85277838274214270114E-13,
 7.18012445138366623367E-13,
-1.79417853150680611778E-12,
-1.32158118404477131188E-11,
-3.14991652796324136454E-11,
 1.18891471078464383424E-11,
 4.94060238822496958910E-10,
 3.39623202570838634515E-9,
 2.26666899049817806459E-8,
 2.04891858946906374183E-7,
 2.89137052083475648297E-6,
 6.88975834691682398426E-5,
 3.36911647825569408990E-3,
 8.04490411014108831608E-1
};
#endif

#ifdef DEC
static unsigned short B[] = {
0122005,0066672,0123124,0054311,
0121662,0033323,0030214,0104602,
0022515,0170300,0113314,0020413,
0022437,0117350,0035402,0007146,
0123243,0000135,0057220,0177435,
0123305,0073476,0144106,0170702,
0023777,0071755,0017527,0154373,
0024211,0052214,0102247,0033270,
0124454,0017763,0171453,0012322,
0125072,0166316,0075505,0154616,
0024612,0133770,0065376,0025045,
0025730,0162143,0056036,0001632,
0026112,0015077,0150464,0063542,
0126374,0101030,0014274,0065457,
0127150,0077271,0125763,0157617,
0127412,0104350,0040713,0120445,
0027121,0023765,0057500,0001165,
0030407,0147146,0003643,0075644,
0031151,0061445,0044422,0156065,
0031702,0132224,0003266,0125551,
0032534,0000076,0147153,0005555,
0033502,0004536,0004016,0026055,
0034620,0076433,0142314,0171215,
0036134,0146145,0013454,0101104,
0040115,0171425,0062500,0047133
};
#endif

#ifdef IBMPC
static unsigned short B[] = {
0x8b19,0x54ca,0xadb7,0xbc60,
0x9130,0x6611,0x46da,0xbc56,
0x8421,0x12d9,0xbe18,0x3c89,
0x41cd,0x0760,0xf3dd,0x3c83,
0x1fe4,0xabd2,0x600b,0xbcb4,
0xde38,0xd908,0xaee7,0xbcb8,
0xfb1f,0xa3ea,0xee7d,0x3cdf,
0xe6d7,0x9094,0x2a91,0x3cf1,
0x629a,0x7e65,0x83fe,0xbd05,
0xbb32,0xcf68,0x5d99,0xbd27,
0xc545,0x0d5f,0x56ff,0x3d11,
0xc073,0x6b83,0x1c8c,0x3d5b,
0x8cec,0xfa26,0x4347,0x3d69,
0x8d66,0x0317,0x9043,0xbd7f,
0x7bf2,0x357e,0x0fd7,0xbdad,
0x7425,0x0839,0x511d,0xbdc1,
0x004f,0xabe8,0x24fe,0x3daa,
0x6f75,0xc0f4,0xf9cc,0x3e00,
0x5b87,0xa922,0x2c64,0x3e2d,
0xd56d,0x80d6,0x5692,0x3e58,
0x616e,0xd9cd,0x8007,0x3e8b,
0xc586,0xc101,0x412b,0x3ec8,
0x9e52,0x7899,0x0fa3,0x3f12,
0x9049,0xa2e5,0x998c,0x3f6b,
0x09cb,0xaca8,0xbe62,0x3fe9
};
#endif

#ifdef MIEEE
static unsigned short B[] = {
0xbc60,0xadb7,0x54ca,0x8b19,
0xbc56,0x46da,0x6611,0x9130,
0x3c89,0xbe18,0x12d9,0x8421,
0x3c83,0xf3dd,0x0760,0x41cd,
0xbcb4,0x600b,0xabd2,0x1fe4,
0xbcb8,0xaee7,0xd908,0xde38,
0x3cdf,0xee7d,0xa3ea,0xfb1f,
0x3cf1,0x2a91,0x9094,0xe6d7,
0xbd05,0x83fe,0x7e65,0x629a,
0xbd27,0x5d99,0xcf68,0xbb32,
0x3d11,0x56ff,0x0d5f,0xc545,
0x3d5b,0x1c8c,0x6b83,0xc073,
0x3d69,0x4347,0xfa26,0x8cec,
0xbd7f,0x9043,0x0317,0x8d66,
0xbdad,0x0fd7,0x357e,0x7bf2,
0xbdc1,0x511d,0x0839,0x7425,
0x3daa,0x24fe,0xabe8,0x004f,
0x3e00,0xf9cc,0xc0f4,0x6f75,
0x3e2d,0x2c64,0xa922,0x5b87,
0x3e58,0x5692,0x80d6,0xd56d,
0x3e8b,0x8007,0xd9cd,0x616e,
0x3ec8,0x412b,0xc101,0xc586,
0x3f12,0x0fa3,0x7899,0x9e52,
0x3f6b,0x998c,0xa2e5,0x9049,
0x3fe9,0xbe62,0xaca8,0x09cb
};
#endif

#ifdef ANSIPROT
extern double exp ( double );
extern double sqrt ( double );
#else
double exp(), sqrt();
#endif

/** Modified Bessel function of order zero.
 *  Cephes Math Library Release 2.8:  June, 2000
 *  Copyright 1984, 1987, 2000 by Stephen L. Moshier
 */
double i0(double x)
{
  double y;
  
  if( x < 0 )
    x = -x;
  if( x <= 8.0 )
    {
      y = (x/2.0) - 2.0;
      return( exp(x) * chbevl( y, A, 30 ) );
    }
  
  return(  exp(x) * chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) ); 
}

/** Computes the inner/dot product \f$x^H x\f$.
 */
double dotproductc(fftw_complex* x, int n)
{
  int k;
  double result=0.0;

  for(k=0;k<n;k++)
    result+=(x[k][0]*x[k][0]+x[k][1]*x[k][1]);

  return result;
}

/** Computes the weighted inner/dot product \f$x^H (w \odot x)\f$.
 */
double dotproductc_w(fftw_complex* x, double* w, int n)
{
  int k;
  double result=0.0;

  for(k=0;k<n;k++)
    result+=w[k]*(x[k][0]*x[k][0]+x[k][1]*x[k][1]);

  return result;
}

/** Copies \f$x \leftarrow y\f$.
 */
void copyc(fftw_complex* x,fftw_complex* y, int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]=y[l][0];
      x[l][1]=y[l][1];
    }
}

/** Copies \f$x \leftarrow a y\f$.
 */
void copyc_a(fftw_complex* x,double a, fftw_complex* y, int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]=a*y[l][0];
      x[l][1]=a*y[l][1];
    }
}

/** Copies \f$x \leftarrow w\odot y\f$.
 */
void copyc_w(fftw_complex* x,double* w,fftw_complex* y, int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]=w[l]*y[l][0];
      x[l][1]=w[l]*y[l][1];
    }
}

/** Updates \f$x \leftarrow a x + y\f$.
 */
void updatec_axpy(fftw_complex* x,double a, fftw_complex* y, int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]=a*x[l][0]+y[l][0];
      x[l][1]=a*x[l][1]+y[l][1];
    }
}

/** Updates \f$x \leftarrow x + a y\f$.
 */
void updatec_xpay(fftw_complex* x,double a, fftw_complex* y, int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]+=a*y[l][0];
      x[l][1]+=a*y[l][1];
    }
}

/** Updates \f$x \leftarrow a x + b y\f$.
 */
void updatec_axpby(fftw_complex* x, double a, fftw_complex* y, double b,
		   int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]=a*x[l][0]+b*y[l][0];
      x[l][1]=a*x[l][1]+b*y[l][1];
    }
}

/** Updates \f$x \leftarrow x + a w\odot y\f$.
 */
void updatec_xpawy(fftw_complex* x,double a, double* w, fftw_complex* y,
		   int n)
{
  int l;

  for(l=0;l<n;l++)
    {
      x[l][0]+=a*w[l]*y[l][0];
      x[l][1]+=a*w[l]*y[l][1];
    }
}

/** computes \f$\frac{\|x_0-x\|_{\infty}}{\|x\|_{\infty}} \f$
 */
double error(double *x0, double *x, int n)
{
  double maximum=0.0, maxdiff=0.0, xd, bx;
  int k;
  
  for (k=0; k<n; k++)
    {
      xd=(double)fabs(x0[k]-x[k]);
      
      bx=(double)fabs(x[k]);

      if(xd>maxdiff)
	maxdiff=xd;
      if(bx>maximum)
	maximum=bx;
    }
  
  return (double)maxdiff/maximum;
}

/** computes \f$\frac{\|x_0-x\|_{\infty}}{\|x\|_{\infty}} \f$
 */
double errorc(fftw_complex *x0, fftw_complex *x, int n)
{
  double maximum=0.0, maxdiff=0.0, xd, bx;
  fftw_complex xda;
  int k;
  
  for (k=0; k<n; k++)
    {
      xda[0]=(double)fabs(x0[k][0]-x[k][0]);
      xda[1]=(double)fabs(x0[k][1]-x[k][1]);
      xd=sqrt(xda[0]*xda[0]+xda[1]*xda[1]);
      
      bx=sqrt(x[k][0]*x[k][0]+x[k][1]*x[k][1]);

      if(xd>maxdiff)
	maxdiff=xd;
      if(bx>maximum)
	maximum=bx;
    }
  
  return (double)maxdiff/maximum;
}

/** computes \f$\frac{\|x_0-x\|_{\infty}}{\|\hat f\|_1} \f$
 */
double E_infty_error_c(fftw_complex *x0, fftw_complex *x, int M,
			  fftw_complex *f_hat, int N_d)
{
  double maxdiff=0.0, xd, bx;
  fftw_complex xda;
  int k,j;

  for (j=0; j<M; j++)
    {
      xda[0]=(double)fabs(x0[j][0]-x[j][0]);
      xda[1]=(double)fabs(x0[j][1]-x[j][1]);
      xd=sqrt(xda[0]*xda[0]+xda[1]*xda[1]);

      if((xd>maxdiff)||isnan(xd))
	maxdiff=xd;
    }

  bx=0.0;
  for (k=0; k<N_d; k++)
    {
      bx=bx+sqrt(f_hat[k][0]*f_hat[k][0]+f_hat[k][1]*f_hat[k][1]);
    }

  return ((double)(maxdiff/bx));
}

/** computes \f$\frac{\|x_0-x\|_2}{\|x_0\|_2} \f$
 */
double E_2_error_c(fftw_complex *x0, fftw_complex *x, int M)
{
  double sumdiff=0.0, xd, bx;
  fftw_complex xda;
  int j;

  for (j=0; j<M; j++)
    {
      xda[0]=(double)fabs(x0[j][0]-x[j][0]);
      xda[1]=(double)fabs(x0[j][1]-x[j][1]);
      xd=xda[0]*xda[0]+xda[1]*xda[1];

      sumdiff+=xd;
    }

  bx=0.0;
  for (j=0; j<M; j++)
    {
      bx=bx+x0[j][0]*x0[j][0]+x0[j][1]*x0[j][1];
    }

  return ((double)(sqrt(sumdiff/bx)));
}

/** vector print
 */
void vpr(double *x, int n, char *text)
{
  int k;

  printf ("\n %s, adr=%p\n", text, (void*)x); 
  for (k=0; k<n; k++)
    {
      if (k%5==0) 
	printf("%d.", k);
      printf("%14.7e,", x[k]);
      if (k%5==4) 
	printf("\n");
      fflush(stdout);
    }
  printf("\n");
}

/** vector print
 */
void vpr_c(fftw_complex *x, int n, char *text)
{
  int k;

  printf ("\n %s, adr=%p\n", text, (void*)x);
  for (k=0; k<n; k++)
    {
      if (k%4==0) 
	printf("%d.", k);
      printf("%5.4e+%5.4e i,", x[k][0],x[k][1]);
      if (k%4==3) 
	printf("\n");
      fflush(stdout);
    }
  printf("\n");
}

/** vector print
 */
void vpr_int(int *x, int n, char *text)
{
  int k;

  printf ("\n %s, adr=%p\n", text, (void*)x);
  for (k=0; k<n; k++)
    {
      if (k%5==0) 
	printf("%d.", k);
      printf("%d,", x[k]);
      if (k%5==4) 
	printf("\n");
      fflush(stdout);
    }
}

/** Computes non periodic voronoi weights 
 *  assumes ordered x_j */
void voronoi_weights_1d(double *w, double *x, int M)
{
  int j;
  
  w[0]=(x[1]-x[0])/2;
  for(j=1;j<M-1;j++)
    w[j]=(x[j+1]-x[j-1])/2;
  w[M-1]=(x[M-1]-x[M-2])/2;
}

/** Computes the damping factor for the modified Fejer kernel.
 *  /f$\frac{2}{N}\left(1-\frac{\left|2k+1\right|}{N}\right)/f$
 */
double modified_fejer(int N,int kk)
{
  double result;

  result=2.0/((double)N*N)*(1-fabs(2.0*kk+1)/((double)N));

  return result;
}

/** Computes the damping factor for the modified Jackson kernel.
 */
double modified_jackson2(int N,int kk)
{
  int kj;
  double result;
  double n=(N/2+1)/2;
  double k;
  
  for(result=0,kj=kk;kj<=kk+1;kj++)
    {
      k=fabs(kj);
      if(k/n<1)
	result+= 1 - (3.0*k + 6.0*n*pow(k,2) - 3.0*pow(k,3)) 
	  / (2.0*n*(2.0*pow(n,2)+1.0));
      else
	result+= (2*n-k)*(pow(2*n-k,2)-1) / (2.0*n*(2.0*pow(n,2)+1.0));
    }
      
  return result;
}

/** Computes the damping factor for the modified generalised Jackson kernel.
 */
double modified_jackson4(int N,int kk)
{
  int kj;
  double result;
  double n=(N/2+3)/4;
  double k;
  double normalisation=(2416*pow(n,7)+1120*pow(n,5)+784*pow(n,3)+720*n);
  
  for(result=0,kj=kk;kj<=kk+1;kj++)
    {
      k=fabs(kj);

      if(k/n<1)
	result+= 1 - (1260*k + (1680*pow(n,5)+2240*pow(n,3)+2940*n)*pow(k,2) -
		      1715*pow(k,3) - (560*pow(n,3)+1400*n)*pow(k,4) + 490*
		      pow(k,5) + 140*n*pow(k,6) - 35*pow(k,7)) / normalisation;
      
      if((1<=k/n)&&(k/n<2))
  	result+= ((2472*pow(n,7)+336*pow(n,5)+3528*pow(n,3)-1296*n) - 
		  (392*pow(n,6)-3920*pow(n,4)+8232*pow(n,2)-756)*k -
		  (504*pow(n,5)+10080*pow(n,3)-5292*n)*pow(k,2) - 
		  (1960*pow(n,4)-7840*pow(n,2)+1029)*pow(k,3) +
		  (2520*pow(n,3)-2520*n)*pow(k,4) - (1176*pow(n,2)-294)*pow(k,5)
		  + 252*n*pow(k,6) - 21*pow(k,7)) / normalisation;
      
      if((2<=k/n)&&(k/n<3))
 	result+= (-(1112*pow(n,7)-12880*pow(n,5)+7448*pow(n,3)-720*n) +
		  (12152*pow(n,6)-27440*pow(n,4)+8232*pow(n,2)-252)*k -
		  (19320*pow(n,5)-21280*pow(n,3)+2940*n)*pow(k,2) +
		  (13720*pow(n,4)-7840*pow(n,2)+343)*pow(k,3) -
		  (5320*pow(n,3)-1400*n)*pow(k,4) + (1176*pow(n,2)-98)*pow(k,5)
		  - 140*n*pow(k,6) + 7*pow(k,7)) / normalisation;
      
      if((3<=k/n)&&(k/n<4))
	result+= ((4*n-k)*(pow(4*n-k,2)-1)*(pow(4*n-k,2)-4)*(pow(4*n-k,2)-9)) /
	  normalisation;
    }

  return result;
}

/** Computes the damping factor for the modified Sobolev kernel.
 */
double modified_sobolev(double mu,int kk)
{
  double result;
  int kj,k;
  
  for(result=0,kj=kk;kj<=kk+1;kj++)
    {
      k=fabs(kj);
      if(k==0)
	result+=1;
      else
	result+=pow(k,-2*mu);
    }

  return result;
}

/** Computes the damping factor for the modified multiquadric kernel.
 */
double modified_multiquadric(double mu,double c,int kk)
{
  double result;
  int kj,k;
  
  for(result=0,kj=kk;kj<=kk+1;kj++)
    {
      k=fabs(kj);
      result+=pow(k*k+c*c,-mu);
    }

  return result;
}