00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef __INFFT_H__
00023 #define __INFFT_H__
00024
00025 #include "config.h"
00026
00027 #include <math.h>
00028 #include <float.h>
00029
00030 #include <stdlib.h>
00031 #include <stdarg.h>
00032 #include <stddef.h>
00033
00034 #if HAVE_SYS_TYPES_H
00035 #include <sys/types.h>
00036 #endif
00037
00038 #if HAVE_STDINT_H
00039 #include <stdint.h>
00040 #endif
00041
00042 #if HAVE_INTTYPES_H
00043 #include <inttypes.h>
00044 #endif
00045
00046 #if HAVE_DISPATCH_DISPATCH_H
00047 #include <dispatch/dispatch.h>
00048 #endif
00049
00050 #if defined(DIRAC_DELTA)
00051 #define PHI_HUT(k,d) 1.0
00052 #define PHI(x,d) (fabs((x))<10e-8)? 1.0 : 0.0
00053 #define WINDOW_HELP_INIT(d)
00054 #define WINDOW_HELP_FINALIZE
00055 #define WINDOW_HELP_ESTIMATE_m {ths->m = 0;}
00056 #elif defined(GAUSSIAN)
00057 #define PHI_HUT(k,d) ((double)exp(-(pow(PI*(k)/ths->n[d],2.0)*ths->b[d])))
00058 #define PHI(x,d) ((double)exp(-pow((x)*ths->n[d],2.0)/ ths->b[d])/sqrt(PI*ths->b[d]))
00059 #define WINDOW_HELP_INIT \
00060 { \
00061 int WINDOW_idx; \
00062 ths->b = (double*) nfft_malloc(ths->d*sizeof(double)); \
00063 for(WINDOW_idx=0; WINDOW_idx<ths->d; WINDOW_idx++) \
00064 ths->b[WINDOW_idx]=((double)2*ths->sigma[WINDOW_idx])/ \
00065 (2*ths->sigma[WINDOW_idx]-1)*(((double)ths->m) / PI); \
00066 }
00067 #define WINDOW_HELP_FINALIZE {nfft_free(ths->b);}
00068 #define WINDOW_HELP_ESTIMATE_m {ths->m =12;}
00069 #elif defined(B_SPLINE)
00070 #define PHI_HUT(k,d) ((double)(((k)==0)? 1.0/ths->n[(d)] : \
00071 pow(sin((k)*PI/ths->n[(d)])/((k)*PI/ths->n[(d)]),2*ths->m)/ths->n[(d)]))
00072 #define PHI(x,d) (nfft_bspline(2*ths->m,((x)*ths->n[(d)])+ \
00073 (double)ths->m,ths->spline_coeffs)/ths->n[(d)])
00074 #define WINDOW_HELP_INIT \
00075 { \
00076 ths->spline_coeffs= (double*)nfft_malloc(2*ths->m*sizeof(double)); \
00077 }
00078 #define WINDOW_HELP_FINALIZE {nfft_free(ths->spline_coeffs);}
00079 #define WINDOW_HELP_ESTIMATE_m {ths->m =11;}
00080 #elif defined(SINC_POWER)
00081 #define PHI_HUT(k,d) (nfft_bspline(2*ths->m,((double)2*ths->m*(k))/ \
00082 ((2*ths->sigma[(d)]-1)*ths->n[(d)]/ths->sigma[(d)])+ (double)ths->m, \
00083 ths->spline_coeffs))
00084 #define PHI(x,d) ((double)(ths->n[(d)]/ths->sigma[(d)]*(2*ths->sigma[(d)]-1)/\
00085 (2*ths->m)*pow(nfft_sinc(PI*ths->n[(d)]/ths->sigma[(d)]*(x)* \
00086 (2*ths->sigma[(d)]-1)/(2*ths->m)),2*ths->m)/ths->n[(d)]))
00087 #define WINDOW_HELP_INIT \
00088 { \
00089 ths->spline_coeffs= (double*)nfft_malloc(2*ths->m*sizeof(double)); \
00090 }
00091 #define WINDOW_HELP_FINALIZE {nfft_free(ths->spline_coeffs);}
00092 #define WINDOW_HELP_ESTIMATE_m {ths->m = 9;}
00093 #else
00094 #define PHI_HUT(k,d) ((double)nfft_i0( ths->m*sqrt(\
00095 pow((double)(ths->b[d]),2.0) - pow(2.0*PI*(k)/ths->n[d],2.0))))
00096 #define PHI(x,d) ((double)((pow((double)(ths->m),2.0)\
00097 -pow((x)*ths->n[d],2.0))>0)? \
00098 sinh(ths->b[d]*sqrt(pow((double)(ths->m),2.0)- \
00099 pow((x)*ths->n[d],2.0)))/(PI*sqrt(pow((double)(ths->m),2.0)- \
00100 pow((x)*ths->n[d],2.0))): (((pow((double)(ths->m),2.0)- \
00101 pow((x)*ths->n[d],2.0))<0)? sin(ths->b[d]* \
00102 sqrt(pow(ths->n[d]*(x),2.0)-pow((double)(ths->m),2.0)))/ \
00103 (PI*sqrt(pow(ths->n[d]*(x),2.0)-pow((double)(ths->m),2.0))):1.0))
00104 #define WINDOW_HELP_INIT \
00105 { \
00106 int WINDOW_idx; \
00107 ths->b = (double*) nfft_malloc(ths->d*sizeof(double)); \
00108 for(WINDOW_idx=0; WINDOW_idx<ths->d; WINDOW_idx++) \
00109 ths->b[WINDOW_idx] = ((double)PI*(2.0-1.0/ths->sigma[WINDOW_idx])); \
00110 }
00111 #define WINDOW_HELP_FINALIZE {nfft_free(ths->b);}
00112 #define WINDOW_HELP_ESTIMATE_m {ths->m = 6;}
00113 #endif
00114
00115
00116 #define CONCAT(prefix, name) prefix ## name
00117 #if defined(NFFT_SINGLE)
00118 typedef float R;
00119 typedef float _Complex C;
00120 #define X(name) CONCAT(nfftf_, name)
00121 #elif defined(NFFT_LDOUBLE)
00122 typedef long double R;
00123 typedef long double _Complex C;
00124 #define X(name) CONCAT(nfftl_, name)
00125 #else
00126 typedef double R;
00127 typedef double _Complex C;
00128 #define X(name) CONCAT(nfft_, name)
00129 #endif
00130
00131 #ifdef NFFT_LDOUBLE
00132 # define K(x) ((R) x##L)
00133 #else
00134 # define K(x) ((R) x)
00135 #endif
00136 #define DK(name, value) const R name = K(value)
00137
00142 typedef ptrdiff_t INT;
00143
00144 #define IF(x,a,b) ((x)?(a):(b))
00145
00146 #define MIN(a,b) (((a)<(b))?(a):(b))
00147 #define MAX(a,b) (((a)>(b))?(a):(b))
00148 #define ABS(x) (((x)>K(0.0))?(x):(-(x)))
00149 #define SIGN(a) (((a)>=0)?1:-1)
00150 #define SIGN(a) (((a)>=0)?1:-1)
00151 #define SIGNF(a) IF((a)<K(0.0),K(-1.0),K(1.0))
00152
00153 #if defined(NFFT_LDOUBLE)
00154 #define LOG1P log1pl
00155 #if HAVE_DECL_LOG1PL == 0
00156 extern long double log1pl(const long double);
00157 #endif
00158 #define LOG10 log10l
00159 #if HAVE_DECL_LOG10L == 0
00160 extern long double log10l(const long double);
00161 #endif
00162 #define SQRT sqrtl
00163 #if HAVE_DECL_SQRTL == 0
00164 extern long double sqrtl(const long double);
00165 #endif
00166 #define TGAMMA tgammal
00167 #if HAVE_DECL_TGAMMAL == 0
00168 extern long double tgammal(const long double);
00169 #endif
00170 #define LGAMMA lgammal
00171 #if HAVE_DECL_LGAMMAL == 0
00172 extern long double lgammal(const long double);
00173 #endif
00174 #define LOG logl
00175 #if HAVE_DECL_LOGL == 0
00176 extern long double logl(const long double);
00177 #endif
00178 #define LOG2 log2l
00179 #if HAVE_DECL_LOG2L == 0
00180 extern long double log2l(const long double);
00181 #endif
00182 #define SIN sinl
00183 #if HAVE_DECL_SINL == 0
00184 extern long double sinl(const long double);
00185 #endif
00186 #define COS cosl
00187 #if HAVE_DECL_COSL == 0
00188 extern long double cosl(const long double);
00189 #endif
00190 #define ACOS acosl
00191 #if HAVE_DECL_ACOSL == 0
00192 extern long double acosl(const long double);
00193 #endif
00194 #define POW powl
00195 #if HAVE_DECL_POWL == 0
00196 extern long double powl(const long double, const long double);
00197 #endif
00198 #define EXP expl
00199 #if HAVE_DECL_EXPL == 0
00200 extern long double expl(const long double);
00201 #endif
00202 #define COPYSIGN copysignl
00203 #if HAVE_DECL_COPYSIGNL == 0
00204 extern long double copysignl(const long double, const long double);
00205 #endif
00206 #define LRINT lrintl
00207 #if HAVE_DECL_LRINTL == 0
00208 extern long int lrintl(const long double);
00209 #endif
00210 #define FMAX fmaxl
00211 #if HAVE_DECL_FMAXL == 0
00212 extern double fmaxl(const long double x, const long double y);
00213 #endif
00214 #elif defined(NFFT_SINGLE)
00215 #define LOG1P log1pf
00216 #if HAVE_DECL_LOG1PF == 0
00217 extern float log1pf(const float);
00218 #endif
00219 #define LOG10 log10f
00220 #if HAVE_DECL_LOG10F == 0
00221 extern float log10f(const float);
00222 #endif
00223 #define SQRT sqrtf
00224 #if HAVE_DECL_SQRTF == 0
00225 extern float sqrtf(const float);
00226 #endif
00227 #define >TGAMMA tgammaf
00228 #if HAVE_DECL_TGAMMAF == 0
00229 extern float tgammaf(const float);
00230 #endif
00231 #define LGAMMA lgammaf
00232 #if HAVE_DECL_LGAMMAF == 0
00233 extern float lgammaf(const float);
00234 #endif
00235 #define LOG logf
00236 #if HAVE_DECL_LOGF == 0
00237 extern float logf(const float);
00238 #endif
00239 #define LOG2 log2f
00240 #if HAVE_DECL_LOG2F == 0
00241 extern float log2f(const float);
00242 #endif
00243 #define SIN sinf
00244 #if HAVE_DECL_SINF == 0
00245 extern float sinf(const float);
00246 #endif
00247 #define COS cosf
00248 #if HAVE_DECL_COSF == 0
00249 extern float cosf(const float);
00250 #endif
00251 #define ACOS acosf
00252 #if HAVE_DECL_ACOSF == 0
00253 extern float acosf(const float);
00254 #endif
00255 #define POW powf
00256 #if HAVE_DECL_POWF == 0
00257 extern float powf(const float, const float);
00258 #endif
00259 #define EXP expf
00260 #if HAVE_DECL_EXPF == 0
00261 extern float expf(const float);
00262 #endif
00263 #define COPYSIGN copysignf
00264 #if HAVE_DECL_COPYSIGNF == 0
00265 extern float copysignf(const float, const float);
00266 #endif
00267 #define LRINT lrintf
00268 #if HAVE_DECL_LRINTF == 0
00269 extern long int lrintf(const float);
00270 #endif
00271 #define FMAX fmaxf
00272 #if HAVE_DECL_FMAXF == 0
00273 extern double fmaxf(const float x, const float y);
00274 #endif
00275 #else
00276 #define LOG1P log1p
00277 #if HAVE_DECL_LOG1P == 0
00278 extern double log1p(const double);
00279 #endif
00280 #define LOG10 log10
00281 #if HAVE_DECL_LOG10 == 0
00282 extern double log10(const double);
00283 #endif
00284 #define SQRT sqrt
00285 #if HAVE_DECL_SQRT == 0
00286 extern double sqrt(const double);
00287 #endif
00288 #define TGAMMA tgamma
00289 #if HAVE_DECL_TGAMMA == 0
00290 extern double tgamma(const double);
00291 #endif
00292 #define LGAMMA lgamma
00293 #if HAVE_DECL_LGAMMA == 0
00294 extern double lgamma(const double);
00295 #endif
00296 #define LOG log
00297 #if HAVE_DECL_LOG == 0
00298 extern double log(const double);
00299 #endif
00300 #define LOG2 log2
00301 #if HAVE_DECL_LOG2 == 0
00302 extern double log2(const double);
00303 #endif
00304 #define SIN sin
00305 #if HAVE_DECL_SIN == 0
00306 extern double sin(const double);
00307 #endif
00308 #define COS cos
00309 #if HAVE_DECL_COS == 0
00310 extern double cos(const double);
00311 #endif
00312 #define ACOS acos
00313 #if HAVE_DECL_ACOS == 0
00314 extern double acos(const double);
00315 #endif
00316 #define POW pow
00317 #if HAVE_DECL_POW == 0
00318 extern double pow(const double, const double);
00319 #endif
00320 #define EXP exp
00321 #if HAVE_DECL_EXP == 0
00322 extern double exp(const double);
00323 #endif
00324 #define COPYSIGN copysign
00325 #if HAVE_DECL_COPYSIGN == 0
00326 extern double copysign(const double, const double);
00327 #endif
00328 #define LRINT lrint
00329 #if HAVE_DECL_LRINT == 0
00330 extern long int lrint(const double);
00331 #endif
00332 #define FMAX fmax
00333 #if HAVE_DECL_FMAX == 0
00334 extern double fmax(const double x, const double y);
00335 #endif
00336 #endif
00337
00338 #if defined(NFFT_LDOUBLE)
00339 # define CEIL ceill
00340 # define FLOOR floorl
00341 # define ROUND roundl
00342 # define FABS fabsl
00343 # define FMAX fmaxl
00344 # define FMIN fminl
00345 # define FREXP frexpl
00346 # define LDEXP ldexpl
00347 # define EPS LDBL_EPSILON
00348 # define R_MIN LDBL_MIN
00349 # define R_MAX LDBL_MAX
00350 # define R_MIN_EXP LDBL_MIN_EXP
00351 # define R_MAX_EXP LDBL_MAX_EXP
00352 # define R_MIN_10_EXP LDBL_MIN_10_EXP
00353 # define R_MAX_10_EXP LDBL_MAX_10_EXP
00354 # define R_DIG LDBL_DIG
00355 # define CREAL creall
00356 # define CIMAG cimagl
00357 #elif defined(NFFT_SINGLE)
00358 # define CEIL ceilf
00359 # define FLOOR floorf
00360 # define ROUND roundf
00361 # define FABS fabsf
00362 # define FMAX fmaxf
00363 # define FMIN fminf
00364 # define FREXP frexpf
00365 # define LDEXP ldexpf
00366 # define EPS FLT_EPSILON
00367 # define R_MIN FLT_MIN
00368 # define R_MAX FLT_MAX
00369 # define R_MIN_EXP FLT_MIN_EXP
00370 # define R_MAX_EXP FLT_MAX_EXP
00371 # define R_MIN_10_EXP FLT_MIN_10_EXP
00372 # define R_MAX_10_EXP FLT_MAX_10_EXP
00373 # define R_DIG FLT_DIG
00374 # define CREAL crealf
00375 # define CIMAG cimagf
00376 #else
00377 # define CEIL ceil
00378 # define FLOOR floor
00379 # define ROUND round
00380 # define FABS fabs
00381 # define FMAX fmax
00382 # define FMIN fmin
00383 # define FREXP frexp
00384 # define LDEXP ldexp
00385 # define EPS DBL_EPSILON
00386 # define R_MIN DBL_MIN
00387 # define R_MAX DBL_MAX
00388 # define R_MIN_EXP DBL_MIN_EXP
00389 # define R_MAX_EXP DBL_MAX_EXP
00390 # define R_MAX_10_EXP DBL_MAX_10_EXP
00391 # define R_MIN_10_EXP DBL_MIN_10_EXP
00392 # define R_DIG DBL_DIG
00393 # define CREAL creal
00394 # define CIMAG cimag
00395 #endif
00396
00397 #if HAVE_DECL_DRAND48 == 0
00398 extern double drand48(void);
00399 #endif
00400 #if HAVE_DECL_SRAND48 == 0
00401 extern void srand48(long int);
00402 #endif
00403 #define RAND ((R)drand48())
00404 #define R_RADIX FLT_RADIX
00405 #define II _Complex_I
00406
00407
00408 #if defined(NFFT_LDOUBLE)
00409 # define FE "LE"
00410 #elif defined(NFFT_SINGLE)
00411 # define FE "E"
00412 #else
00413 # define FE "lE"
00414 #endif
00415
00416 #define TRUE 1
00417 #define FALSE 0
00418
00419 #define KPI K(3.1415926535897932384626433832795028841971693993751)
00420 #define K2PI K(6.2831853071795864769252867665590057683943387987502)
00421
00423 #define UNUSED(x) (void)x
00424
00425 #endif