00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifndef IPPL_FFT_FFTPACK_FFT_H
00012 #define IPPL_FFT_FFTPACK_FFT_H
00013
00014
00015 #include "Utility/PAssert.h"
00016 #include "Utility/IpplInfo.h"
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #if defined(IPPL_T3E)
00028
00029 #define cffti_ CFFTI
00030 #define cfftf_ CFFTF
00031 #define cfftb_ CFFTB
00032 #define rffti_ RFFTI
00033 #define rfftf_ RFFTF
00034 #define rfftb_ RFFTB
00035 #define sinti_ SINTI
00036 #define sint_ SINT
00037 #define fcffti_ FCFFTI
00038 #define fcfftf_ FCFFTF
00039 #define fcfftb_ FCFFTB
00040 #define frffti_ FRFFTI
00041 #define frfftf_ FRFFTF
00042 #define frfftb_ FRFFTB
00043 #define fsinti_ FSINTI
00044 #define fsint_ FSINT
00045
00046 #endif
00047
00048
00049 #if defined(IPPL_SP2)
00050
00051 #define cffti_ cffti
00052 #define cfftf_ cfftf
00053 #define cfftb_ cfftb
00054 #define rffti_ rffti
00055 #define rfftf_ rfftf
00056 #define rfftb_ rfftb
00057 #define sinti_ sinti
00058 #define sint_ sint
00059 #define fcffti_ fcffti
00060 #define fcfftf_ fcfftf
00061 #define fcfftb_ fcfftb
00062 #define frffti_ frffti
00063 #define frfftf_ frfftf
00064 #define frfftb_ frfftb
00065 #define fsinti_ fsinti
00066 #define fsint_ fsint
00067
00068 #endif
00069
00070
00071 extern "C" {
00072
00073 void cffti_(int& n, double& wsave);
00074 void cfftf_(int& n, double& r, double& wsave);
00075 void cfftb_(int& n, double& r, double& wsave);
00076
00077 void rffti_(int& n, double& wsave);
00078 void rfftf_(int& n, double& r, double& wsave);
00079 void rfftb_(int& n, double& r, double& wsave);
00080
00081 void sinti_(int& n, double& wsave);
00082 void sint_(int& n, double& r, double& wsave);
00083
00084 void fcffti_(int& n, float& wsave);
00085 void fcfftf_(int& n, float& r, float& wsave);
00086 void fcfftb_(int& n, float& r, float& wsave);
00087
00088 void frffti_(int& n, float& wsave);
00089 void frfftf_(int& n, float& r, float& wsave);
00090 void frfftb_(int& n, float& r, float& wsave);
00091
00092 void fsinti_(int& n, float& wsave);
00093 void fsint_(int& n, float& r, float& wsave);
00094 }
00095
00096
00097
00098
00099 template <class T>
00100 class FFTPACK_wrap {};
00101
00102
00103 template <>
00104 class FFTPACK_wrap<float> {
00105
00106 public:
00107
00108
00109
00110 static void ccffti(int n, float* wsave) { fcffti_(n, *wsave); }
00111 static void rcffti(int n, float* wsave) { frffti_(n, *wsave); }
00112 static void rrffti(int n, float* wsave) { fsinti_(n, *wsave); }
00113
00114 static void ccfftf(int n, float* r, float* wsave) { fcfftf_(n, *r, *wsave); }
00115 static void ccfftb(int n, float* r, float* wsave) { fcfftb_(n, *r, *wsave); }
00116
00117 static void rcfftf(int n, float* r, float* wsave) { frfftf_(n, *r, *wsave); }
00118 static void rcfftb(int n, float* r, float* wsave) { frfftb_(n, *r, *wsave); }
00119
00120 static void rrfft(int n, float* r, float* wsave) { fsint_(n, *r, *wsave); }
00121
00122 };
00123
00124
00125 template <>
00126 class FFTPACK_wrap<double> {
00127
00128 public:
00129
00130
00131
00132 static void ccffti(int n, double* wsave) { cffti_(n, *wsave); }
00133 static void rcffti(int n, double* wsave) { rffti_(n, *wsave); }
00134 static void rrffti(int n, double* wsave) { sinti_(n, *wsave); }
00135
00136 static void ccfftf(int n, double* r, double* wsave) {cfftf_(n, *r, *wsave);}
00137 static void ccfftb(int n, double* r, double* wsave) {cfftb_(n, *r, *wsave);}
00138
00139 static void rcfftf(int n, double* r, double* wsave) {rfftf_(n, *r, *wsave);}
00140 static void rcfftb(int n, double* r, double* wsave) {rfftb_(n, *r, *wsave);}
00141
00142 static void rrfft(int n, double* r, double* wsave) { sint_(n, *r, *wsave); }
00143
00144 };
00145
00146
00147
00148 template <class T>
00149 class FFTPACK {
00150
00151 public:
00152
00153
00154 #ifdef IPPL_HAS_TEMPLATED_COMPLEX
00155 typedef complex<T> Complex_t;
00156 #else
00157 typedef complex Complex_t;
00158 #endif
00159
00160
00161 FFTPACK(void) {}
00162
00163
00164 ~FFTPACK(void);
00165
00166
00167
00168
00169 void setup(unsigned numTransformDims, const int* transformTypes,
00170 const int* axisLengths);
00171
00172
00173 void callFFT(unsigned transformDim, int direction, Complex_t* data);
00174
00175
00176 void callFFT(unsigned transformDim, int direction, T* data);
00177
00178 private:
00179
00180 unsigned numTransformDims_m;
00181 int* transformType_m;
00182 int* axisLength_m;
00183 T** trig_m;
00184
00185 };
00186
00187
00188
00189
00190
00191 template <class T>
00192 inline
00193 FFTPACK<T>::~FFTPACK(void) {
00194
00195 for (unsigned d=0; d<numTransformDims_m; ++d)
00196 delete [] trig_m[d];
00197 delete [] trig_m;
00198 delete [] transformType_m;
00199 delete [] axisLength_m;
00200 }
00201
00202
00203 template <class T>
00204 inline void
00205 FFTPACK<T>::setup(unsigned numTransformDims, const int* transformTypes,
00206 const int* axisLengths) {
00207
00208
00209 numTransformDims_m = numTransformDims;
00210 transformType_m = new int[numTransformDims_m];
00211 axisLength_m = new int[numTransformDims_m];
00212 unsigned d;
00213 for (d=0; d<numTransformDims_m; ++d) {
00214 transformType_m[d] = transformTypes[d];
00215 axisLength_m[d] = axisLengths[d];
00216 }
00217
00218
00219 trig_m = new T*[numTransformDims_m];
00220 for (d=0; d<numTransformDims_m; ++d) {
00221 switch (transformType_m[d]) {
00222 case 0:
00223 trig_m[d] = new T[4 * axisLength_m[d] + 15];
00224 FFTPACK_wrap<T>::ccffti(axisLength_m[d], trig_m[d]);
00225 break;
00226 case 1:
00227 trig_m[d] = new T[2 * axisLength_m[d] + 15];
00228 FFTPACK_wrap<T>::rcffti(axisLength_m[d], trig_m[d]);
00229 break;
00230 case 2:
00231 trig_m[d] = new T[static_cast<int>(2.5 * axisLength_m[d] + 0.5) + 15];
00232 FFTPACK_wrap<T>::rrffti(axisLength_m[d], trig_m[d]);
00233 break;
00234 default:
00235 ERRORMSG("Unknown transform type requested!!" << endl);
00236 break;
00237 }
00238 }
00239
00240 return;
00241 }
00242
00243
00244 template <class T>
00245 inline void
00246 FFTPACK<T>::callFFT(unsigned transformDim, int direction,
00247 FFTPACK<T>::Complex_t* data) {
00248
00249
00250 PAssert(transformDim<numTransformDims_m);
00251 PAssert(direction==+1 || direction==-1);
00252
00253
00254 T* rdata = reinterpret_cast<T*>(data);
00255
00256
00257 switch (transformType_m[transformDim]) {
00258 case 0:
00259 if (direction == +1) {
00260
00261 FFTPACK_wrap<T>::ccfftf(axisLength_m[transformDim], rdata,
00262 trig_m[transformDim]);
00263 }
00264 else {
00265
00266 FFTPACK_wrap<T>::ccfftb(axisLength_m[transformDim], rdata,
00267 trig_m[transformDim]);
00268 }
00269 break;
00270 case 1:
00271 if (direction == +1) {
00272
00273 FFTPACK_wrap<T>::rcfftf(axisLength_m[transformDim], rdata,
00274 trig_m[transformDim]);
00275
00276 int clen = axisLength_m[transformDim]/2 + 1;
00277 data[clen-1] = Complex_t(imag(data[clen-2]),0.0);
00278 for (int i = clen-2; i > 0; --i)
00279 data[i] = Complex_t(imag(data[i-1]),real(data[i]));
00280 data[0] = Complex_t(real(data[0]),0.0);
00281 }
00282 else {
00283
00284 int clen = axisLength_m[transformDim]/2 + 1;
00285 data[0] = Complex_t(real(data[0]),real(data[1]));
00286 for (int i = 1; i < clen-1; ++i)
00287 data[i] = Complex_t(imag(data[i]),real(data[i+1]));
00288
00289 FFTPACK_wrap<T>::rcfftb(axisLength_m[transformDim], rdata,
00290 trig_m[transformDim]);
00291 }
00292 break;
00293 case 2:
00294 ERRORMSG("Input for real-to-real FFT should be real!!" << endl);
00295 break;
00296 default:
00297 ERRORMSG("Unknown transform type requested!!" << endl);
00298 break;
00299 }
00300
00301 return;
00302 }
00303
00304
00305 template <class T>
00306 inline void
00307 FFTPACK<T>::callFFT(unsigned transformDim, int direction, T* data) {
00308
00309
00310 PAssert(transformDim<numTransformDims_m);
00311 PAssert(direction==+1 || direction==-1);
00312
00313
00314 switch (transformType_m[transformDim]) {
00315 case 0:
00316 ERRORMSG("Input for complex-to-complex FFT should be complex!!" << endl);
00317 break;
00318 case 1:
00319 ERRORMSG("real-to-complex FFT uses complex input!!" << endl);
00320 break;
00321 case 2:
00322
00323 FFTPACK_wrap<T>::rrfft(axisLength_m[transformDim], data,
00324 trig_m[transformDim]);
00325 break;
00326 default:
00327 ERRORMSG("Unknown transform type requested!!" << endl);
00328 break;
00329 }
00330
00331 return;
00332 }
00333
00334
00335 #endif // IPPL_FFT_FFTPACK_FFT_H
00336
00337
00338
00339
00340
00341
00342