#include "TFFTComplexReal.h"
#include "fftw3.h"
#include "TComplex.h"
ClassImp(TFFTComplexReal)
TFFTComplexReal::TFFTComplexReal()
{
fIn = 0;
fOut = 0;
fPlan = 0;
fN = 0;
}
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;
}
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<fNdim; i++){
fN[i] = n[i];
fTotalSize*=n[i];
}
Int_t sizein = Int_t(Double_t(fTotalSize)*(n[ndim-1]/2+1)/n[ndim-1]);
if (!inPlace){
fIn = fftw_malloc(sizeof(fftw_complex)*sizein);
fOut = fftw_malloc(sizeof(Double_t)*fTotalSize);
} else {
fIn = fftw_malloc(sizeof(fftw_complex)*sizein);
fOut = 0;
}
}
TFFTComplexReal::~TFFTComplexReal()
{
fftw_destroy_plan((fftw_plan)fPlan);
fPlan = 0;
fftw_free((fftw_complex*)fIn);
if (fOut)
fftw_free(fOut);
fIn = 0;
fOut = 0;
if (fN)
delete [] fN;
fN = 0;
}
void TFFTComplexReal::Init( Option_t *flags, Int_t ,const Int_t* )
{
fFlags = flags;
if (fOut)
fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN,(fftw_complex*)fIn,(Double_t*)fOut, MapFlag(flags));
else
fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN, (fftw_complex*)fIn, (Double_t*)fIn, MapFlag(flags));
}
void TFFTComplexReal::Transform()
{
if (fPlan)
fftw_execute((fftw_plan)fPlan);
else {
Error("Transform", "transform was not initialized");
return;
}
}
void TFFTComplexReal::GetPoints(Double_t *data, Bool_t fromInput) const
{
if (fromInput){
Error("GetPoints", "the input array has been destroyed");
return;
}
if (fOut){
for (Int_t i=0; i<fTotalSize; i++)
data[i] = ((Double_t*)fOut)[i];
}
else{
for (Int_t i=0; i<fTotalSize; i++)
data[i] = ((Double_t*)fOut)[i];
}
}
Double_t TFFTComplexReal::GetPointReal(Int_t ipoint, Bool_t fromInput) const
{
if (fOut && !fromInput){
return ((Double_t*)fOut)[ipoint];
}
else if (fromInput){
Error("GetPointReal", "Input array was destroyed");
return 0;
}
else
return ((Double_t*)fIn)[ipoint];
}
Double_t TFFTComplexReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
{
Int_t ireal = ipoint[0];
for (Int_t i=0; i<fNdim-1; i++)
ireal=fN[i+1]*ireal + ipoint[i+1];
if (fOut && !fromInput){
return ((Double_t*)fOut)[ireal];
}
else if (fromInput){
Error("GetPointReal", "Input array was destroyed");
return 0;
}
else
return ((Double_t*)fIn)[ireal];
}
void TFFTComplexReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
{
if (fromInput){
Error("GetPointReal", "Input array was destroyed");
return;
} else {
if (fOut){
re = ((Double_t*)fOut)[ipoint];
im = 0;
} else {
re = ((Double_t*)fIn)[ipoint];
im = 0;
}
}
}
void TFFTComplexReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
{
Int_t ireal = ipoint[0];
for (Int_t i=0; i<fNdim-1; i++)
ireal=fN[i+1]*ireal + ipoint[i+1];
if (fromInput){
Error("GetPointReal", "Input array was destroyed");
return;
} else {
if (fOut){
re = ((Double_t*)fOut)[ireal];
im = 0;
} else {
re = ((Double_t*)fIn)[ireal];
im = 0;
}
}
}
Double_t* TFFTComplexReal::GetPointsReal(Bool_t fromInput) const
{
if (fromInput){
Error("GetPointsReal", "Input array has been destroyed");
return 0;
}
if (fOut)
return (Double_t*)fOut;
else
return (Double_t*)fIn;
}
void TFFTComplexReal::GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput) const
{
if (fromInput){
Error("GetPointsComplex", "Input array has been destroyed");
return;
}
if (fOut){
for (Int_t i=0; i<fTotalSize; i++){
re[i] = ((Double_t*)fOut)[i];
im[i] = 0;
}
}
else {
for (Int_t i=0; i<fTotalSize; i++){
re[i] = ((Double_t*)fIn)[i];
im[i] = 0;
}
}
}
void TFFTComplexReal::GetPointsComplex(Double_t *data, Bool_t fromInput) const
{
if (fromInput){
Error("GetPointsComplex", "Input array has been destroyed");
return;
}
if (fOut){
for (Int_t i=0; i<fTotalSize; i+=2){
data[i] = ((Double_t*)fOut)[i/2];
data[i+1]=0;
}
} else {
for (Int_t i=0; i<fTotalSize; i+=2){
data[i] = ((Double_t*)fIn)[i/2];
data[i+1]=0;
}
}
}
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;
}
}
void TFFTComplexReal::SetPoint(const Int_t *ipoint, Double_t re, Double_t im)
{
Int_t ireal = ipoint[0];
for (Int_t i=0; i<fNdim-2; i++){
ireal=fN[i+1]*ireal + ipoint[i+1];
}
ireal = (fN[fNdim-1]/2+1)*ireal+ipoint[fNdim-1];
Int_t realN = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
if (ireal > realN){
Error("SetPoint", "Illegal index value");
return;
}
((fftw_complex*)fIn)[ireal][0] = re;
((fftw_complex*)fIn)[ireal][1] = 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*)fOut)[2*(fN[0]/2)-ipoint][0] = c.Re();
((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][1] = -c.Im();
}
}
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];
}
}
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<sizein; i++){
((fftw_complex*)fIn)[i][0]=re[i];
((fftw_complex*)fIn)[i][1]=im[i];
}
}
UInt_t TFFTComplexReal::MapFlag(Option_t *flag)
{
TString opt = flag;
opt.ToUpper();
if (opt.Contains("ES"))
return FFTW_ESTIMATE;
if (opt.Contains("M"))
return FFTW_MEASURE;
if (opt.Contains("P"))
return FFTW_PATIENT;
if (opt.Contains("EX"))
return FFTW_EXHAUSTIVE;
return FFTW_ESTIMATE;
}
Last change: Wed Jun 25 08:38:54 2008
Last generated: 2008-06-25 08:38
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.