// @(#)root/fft:$Id$ // Author: Anna Kreshuk 07/4/2006 /************************************************************************* * Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. * * All rights reserved. * * * * For the licensing terms see $ROOTSYS/LICENSE. * * For the list of contributors see $ROOTSYS/README/CREDITS. * *************************************************************************/ //////////////////////////////////////////////////////////////////////////////// /// \class TFFTComplexReal /// /// One of the interface classes to the FFTW package, can be used directly /// or via the TVirtualFFT class. Only the basic interface of FFTW is implemented. /// /// Computes the inverse of the real-to-complex transforms (class TFFTRealComplex) /// taking complex input (storing the non-redundant half of a logically Hermitian array) /// to real output (see FFTW manual for more details) /// /// How to use it: /// 1. Create an instance of TFFTComplexReal - this will allocate input and output /// arrays (unless an in-place transform is specified) /// 2. Run the Init() function with the desired flags and settings /// 3. Set the data (via SetPoints(), SetPoint() or SetPointComplex() functions) /// 4. Run the Transform() function /// 5. Get the output (via GetPoints(), GetPoint() or GetPointReal() functions) /// 6. Repeat steps 3)-5) as needed /// /// For a transform of the same size, but with different flags, rerun the Init() /// function and continue with steps 3)-5) /// /// NOTE: /// 1. running Init() function will overwrite the input array! Don't set any data /// before running the Init() function /// 2. FFTW computes unnormalized transform, so doing a transform followed by /// its inverse will lead to the original array scaled by the transform size /// 3. In Complex to Real transform the input array is destroyed. It cannot then /// be retrieved when using the Get's methods. /// //////////////////////////////////////////////////////////////////////////////// #include "TFFTComplexReal.h" #include "fftw3.h" #include "TComplex.h" ClassImp(TFFTComplexReal); //////////////////////////////////////////////////////////////////////////////// ///default TFFTComplexReal::TFFTComplexReal() { fIn = 0; fOut = 0; fPlan = 0; fN = 0; fTotalSize = 0; fNdim = 0; } //////////////////////////////////////////////////////////////////////////////// ///For 1d transforms ///Allocates memory for the input array, and, if inPlace = kFALSE, for the output array TFFTComplexReal::TFFTComplexReal(Int_t n, Bool_t inPlace) { if (!inPlace){ fIn = fftw_malloc(sizeof(fftw_complex)*(n/2+1)); fOut = fftw_malloc(sizeof(Double_t)*n); } else { fIn = fftw_malloc(sizeof(fftw_complex)*(n/2+1)); fOut = 0; } fNdim = 1; fN = new Int_t[1]; fN[0] = n; fPlan = 0; fTotalSize = n; } //////////////////////////////////////////////////////////////////////////////// ///For ndim-dimensional transforms ///Second argument contains sizes of the transform in each dimension TFFTComplexReal::TFFTComplexReal(Int_t ndim, Int_t *n, Bool_t inPlace) { fNdim = ndim; fTotalSize = 1; fN = new Int_t[fNdim]; for (Int_t i=0; i n/2, the according ///point before n/2 is set to (re, -im) void TFFTComplexReal::SetPoint(Int_t ipoint, Double_t re, Double_t im) { if (ipoint <= fN[0]/2){ ((fftw_complex*)fIn)[ipoint][0] = re; ((fftw_complex*)fIn)[ipoint][1] = im; } else { ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][0] = re; ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][1] = -im; } } //////////////////////////////////////////////////////////////////////////////// ///Set the point #ipoint. Since the input is Hermitian, only the first (roughly)half of ///the points have to be set. void TFFTComplexReal::SetPoint(const Int_t *ipoint, Double_t re, Double_t im) { Int_t ireal = ipoint[0]; for (Int_t i=0; i realN){ Error("SetPoint", "Illegal index value"); return; } ((fftw_complex*)fIn)[ireal][0] = re; ((fftw_complex*)fIn)[ireal][1] = im; } //////////////////////////////////////////////////////////////////////////////// ///since the input must be complex-Hermitian, if the ipoint > n/2, the according ///point before n/2 is set to (re, -im) void TFFTComplexReal::SetPointComplex(Int_t ipoint, TComplex &c) { if (ipoint <= fN[0]/2){ ((fftw_complex*)fIn)[ipoint][0] = c.Re(); ((fftw_complex*)fIn)[ipoint][1] = c.Im(); } else { ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][0] = c.Re(); ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][1] = -c.Im(); } } //////////////////////////////////////////////////////////////////////////////// ///set all points. the values are copied. points should be ordered as follows: ///[re_0, im_0, re_1, im_1, ..., re_n, im_n) void TFFTComplexReal::SetPoints(const Double_t *data) { Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]); for (Int_t i=0; i<2*(sizein); i+=2){ ((fftw_complex*)fIn)[i/2][0]=data[i]; ((fftw_complex*)fIn)[i/2][1]=data[i+1]; } } //////////////////////////////////////////////////////////////////////////////// ///Set all points. The values are copied. void TFFTComplexReal::SetPointsComplex(const Double_t *re, const Double_t *im) { Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]); for (Int_t i=0; i