Open citibeth opened 8 years ago
Might be somehow relevant - an example Blitz++/HDF5 combination code: https://github.com/igfuw/libmpdataxx/blob/master/libmpdata%2B%2B/output/hdf5.hpp
A relevant comment by "Shakes" (@shakes76) from SF (https://sourceforge.net/p/blitz/patches/2/):
The following is the code for getting FFT's on Blitz Arrays using the mature and highly recognised FFTW C Library.
template <typename type, int size> void FourierTransform<type,size>::FFTW(Array<complex<double>,2> &field, Array<complex<double>,2> &fftOfField) { complex<double> *ptrField, *ptrFFT; fftw_plan plan;
ptrField = field.data(); ptrFFT = fftOfField.data(); plan = fftw_plan_dft_2d(field.rows(),field.cols(),reinterpret_cast<fftw_complex>(ptrField),reinterpret_cast<fftw_complex>(ptrFFT),FFTW_FORWARD,FFTW_ESTIMATE);
fftw_execute(plan);
fftw_destroy_plan(plan); }
> I have tested and its does work. Note that this works because as FFTW points out:
>
> "C++ has its own complex<T> template class, defined in the standard <complex> header file. Reportedly, the C++ standards committee has recently agreed to mandate that the storage format used for this type be binary-compatible with the C99 type, i.e. an array T[2] with consecutive real [0] and imaginary [1] parts. (See report WG21/N1388.) Although not part of the official standard as of this writing, the proposal stated that: “This solution has been tested with all current major implementations of the standard library and shown to
be working.”
>
> To the extent that this is true, if you have a variable complex<double> *x, you can pass it directly to FFTW via reinterpret_cast<fftw_complex*>(x)."
>
> Also data reordering has to be done if DC is wanted in the centre (http://www.fftw.org/faq/section3.html#centerorigin).
> I have done it in the following:
```c++
void
FourierTransform<type,size>::centerDC(Array<complex<double>,2>
&field, Array<complex<double>,2> &result)
{
Array<complex<double>,2> tmp(field.rows(),field.cols());
for(int k = 0; k < field.cols(); k ++)
for(int j = 0; j < field.rows(); j ++)
tmp((j+field.rows()/2)%field.rows(),(k+field.cols()/2)%field.cols())
= field(j,k);
result = tmp;
}
Code Excerpt is taken from the FourierTransform Class in the up and coming Quantum Mechanics and Discrete Geometry Toolkit based on Blitz++ and Qt - qC++ (qcplusplus.sourceforge.net). It currently supports BEC simulations and the Fourier class can be found in the Subversion area.
Keep up the good work guys, more work on the tiny stuff would be appreciated :)
Cheers Shakes
Patrick Guio wrote:
I like your utilities idea. I also have some wrapper functions to read/write HDF4/5 files and use FFTW 1D transforms (in/out of place, real/complex). The Fits file format for multidimensional array might also be of interest to consider.