00001
00022 #ifndef DGVFOURIERTRANSFORM_H
00023 #define DGVFOURIERTRANSFORM_H
00024
00025 #include <complex>
00026 #include <fftw3.h>
00027
00028 #include "DGVAliases.h"
00029
00048 template<typename type, int rank>
00049 class DGV_TEMPLATEDLL DGVFourierTransform
00050 {
00051 public:
00052
00058 DGVFourierTransform();
00062 virtual ~DGVFourierTransform();
00064
00068 void FFT(Array<complex<type>,rank> &field, Array<complex<type>,rank> &fftOfField);
00072 void iFFT(Array<complex<type>,rank> &fftOfField, Array<complex<type>,rank> &field);
00073
00077 void FFT_Real2Complex(Array<type,rank> &field, Array<complex<type>,rank> &fftOfField);
00081 void FFT_Complex2Real(Array<complex<type>,rank> &fftOfField, Array<type,rank> &field);
00082
00086 void normalise(Array<complex<type>,rank> &field);
00090 void normalise(Array<type,rank> &field);
00094 Array<type,rank> shiftDC(Array<type,rank> &field);
00098 Array<complex<type>,rank> shiftDC(Array<complex<type>,rank> &field);
00102 inline void reset()
00103 { mode = FFTW_ESTIMATE; fftw_cleanup(); }
00104
00105 protected:
00106 unsigned mode;
00107 fftw_plan plan;
00108
00112 void execute();
00113 };
00114
00115 template<typename type, int rank>
00116 DGVFourierTransform<type,rank>::DGVFourierTransform()
00117 {
00118 mode = FFTW_ESTIMATE;
00119 }
00120
00121 template<typename type, int rank>
00122 DGVFourierTransform<type,rank>::~DGVFourierTransform()
00123 {
00124 fftw_cleanup();
00125 }
00126
00127 template<typename type, int rank>
00128 void DGVFourierTransform<type,rank>::FFT(Array<complex<type>,rank> &field, Array<complex<type>,rank> &fftOfField)
00129 {
00130 if(rank < 4 && rank > 0)
00131 {
00132 const int extent[3] = {field.rows(),field.cols(),field.depth()};
00134
00135
00136 plan = fftw_plan_dft(rank,extent,reinterpret_cast<fftw_complex*>(field.data()),reinterpret_cast<fftw_complex*>(fftOfField.data()),FFTW_FORWARD,mode);
00137
00138 execute();
00139 }
00140 }
00141
00142 template<typename type, int rank>
00143 void DGVFourierTransform<type,rank>::iFFT(Array<complex<type>,rank> &fftOfField, Array<complex<type>,rank> &field)
00144 {
00145 if(rank < 4 && rank > 0)
00146 {
00147 const int extent[3] = {fftOfField.rows(),fftOfField.cols(),fftOfField.depth()};
00149
00150
00151 plan = fftw_plan_dft(rank,extent,reinterpret_cast<fftw_complex*>(fftOfField.data()),reinterpret_cast<fftw_complex*>(field.data()),FFTW_BACKWARD,mode);
00152
00153 execute();
00154 }
00155 }
00156
00157 template<typename type, int rank>
00158 void DGVFourierTransform<type,rank>::FFT_Real2Complex(Array<type,rank> &field, Array<complex<type>,rank> &fftOfField)
00159 {
00160 if(rank < 4 && rank > 0)
00161 {
00162 const int extent[3] = {field.rows(),field.cols(),field.depth()};
00164
00165
00166 plan = fftw_plan_dft_r2c(rank,extent,field.data(),reinterpret_cast<fftw_complex*>(fftOfField.data()),mode);
00167
00168 execute();
00169 }
00170 }
00171
00172 template<typename type, int rank>
00173 void DGVFourierTransform<type,rank>::FFT_Complex2Real(Array<complex<type>,rank> &fftOfField, Array<type,rank> &field)
00174 {
00175 if(rank < 4 && rank > 0)
00176 {
00177 const int extent[3] = {fftOfField.rows(),fftOfField.cols(),fftOfField.depth()};
00179
00180
00181 plan = fftw_plan_dft_c2r(rank,extent,reinterpret_cast<fftw_complex*>(fftOfField.data()),field.data(),mode);
00182
00183 execute();
00184 }
00185 }
00186
00187 template<typename type, int rank>
00188 void DGVFourierTransform<type,rank>::normalise(Array<complex<type>,rank> &field)
00189 {
00190 type factor;
00191
00192 if(rank == 1)
00193 factor = static_cast<type>(field.size());
00194 else if(rank == 2)
00195 factor = static_cast<type>(field.rows()*field.cols());
00196 else
00197 factor = static_cast<type>(field.rows()*field.cols()*field.depth());
00198
00199
00200 field /= factor;
00201 }
00202
00203 template<typename type, int rank>
00204 void DGVFourierTransform<type,rank>::normalise(Array<type,rank> &field)
00205 {
00206 type factor;
00207
00208 if(rank == 1)
00209 factor = static_cast<type>(field.size());
00210 else if(rank == 2)
00211 factor = static_cast<type>(field.rows()*field.cols());
00212 else
00213 factor = static_cast<type>(field.rows()*field.cols()*field.depth());
00214
00215
00216 field /= factor;
00217 }
00218
00219 template<typename type, int rank>
00220 Array<type,rank> DGVFourierTransform<type,rank>::shiftDC(Array<type,rank> &field)
00221 {
00222 int rowsHalf = static_cast<int>(field.rows()/2.0), colsHalf = static_cast<int>(field.cols()/2.0);
00223 Array<type,rank> tmpField(field.shape());
00224
00225 for(int j = 0; j < field.rows(); j ++)
00226 for(int k = 0; k < field.cols(); k ++)
00227 tmpField( (j+rowsHalf)%field.rows() , (k+colsHalf)%field.cols() ) = field(j,k);
00228
00229 return tmpField;
00230 }
00231
00232 template<typename type, int rank>
00233 Array<complex<type>,rank> DGVFourierTransform<type,rank>::shiftDC(Array<complex<type>,rank> &field)
00234 {
00235 int rowsHalf = static_cast<int>(field.rows()/2.0), colsHalf = static_cast<int>(field.cols()/2.0);
00236 Array<complex<type>,rank> tmpField(field.shape());
00237
00238 for(int j = 0; j < field.rows(); j ++)
00239 for(int k = 0; k < field.cols(); k ++)
00240 tmpField( (j+rowsHalf)%field.rows() , (k+colsHalf)%field.cols() ) = field(j,k);
00241
00242 return tmpField;
00243 }
00244
00245
00246 template<typename type, int rank>
00247 void DGVFourierTransform<type,rank>::execute()
00248 {
00249 if(plan != NULL)
00250 {
00251 fftw_execute(plan);
00252 fftw_destroy_plan(plan);
00253 }
00254 }
00255
00256 #endif // DGVFOURIERTRANSFORM_H