// @(#)root/base:$Name: $:$Id: TVirtualFFT.cxx,v 1.6 2006/11/30 23:19:47 pcanal Exp $ // Author: Anna Kreshuk 10/04/2006 /************************************************************************* * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. * * All rights reserved. * * * * For the licensing terms see $ROOTSYS/LICENSE. * * For the list of contributors see $ROOTSYS/README/CREDITS. * *************************************************************************/ ////////////////////////////////////////////////////////////////////////// // // TVirtualFFT // // TVirtualFFT is an interface class for Fast Fourier Transforms. // // // // The default FFT library is FFTW. To use it, FFTW3 library should already // be installed, and ROOT should be have fftw3 module enabled, with the directories // of fftw3 include file and library specified (see installation instructions). // Function SetDefaultFFT() allows to change the default library. // // Available transform types: // FFT: // - "C2CFORWARD" - a complex input/output discrete Fourier transform (DFT) // in one or more dimensions, -1 in the exponent // - "C2CBACKWARD"- a complex input/output discrete Fourier transform (DFT) // in one or more dimensions, +1 in the exponent // - "R2C" - a real-input/complex-output discrete Fourier transform (DFT) // in one or more dimensions, // - "C2R" - inverse transforms to "R2C", taking complex input // (storing the non-redundant half of a logically Hermitian array) // to real output // - "R2HC" - a real-input DFT with output in ¡Èhalfcomplex¡É format, // i.e. real and imaginary parts for a transform of size n stored as // r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1 // - "HC2R" - computes the reverse of FFTW_R2HC, above // - "DHT" - computes a discrete Hartley transform // // Sine/cosine transforms: // Different types of transforms are specified by parameter kind of the SineCosine() static // function. 4 different kinds of sine and cosine transforms are available // DCT-I (REDFT00 in FFTW3 notation)- kind=0 // DCT-II (REDFT10 in FFTW3 notation)- kind=1 // DCT-III(REDFT01 in FFTW3 notation)- kind=2 // DCT-IV (REDFT11 in FFTW3 notation)- kind=3 // DST-I (RODFT00 in FFTW3 notation)- kind=4 // DST-II (RODFT10 in FFTW3 notation)- kind=5 // DST-III(RODFT01 in FFTW3 notation)- kind=6 // DST-IV (RODFT11 in FFTW3 notation)- kind=7 // Formulas and detailed descriptions can be found in the chapter // "What FFTW really computes" of the FFTW manual // // NOTE: FFTW computes unnormalized transforms, so doing a transform, followed by its // inverse will give the original array, multiplied by normalization constant // (transform size(N) for FFT, 2*(N-1) for DCT-I, 2*(N+1) for DST-I, 2*N for // other sine/cosine transforms) // // How to use it: // Call to the static function FFT returns a pointer to a fast fourier transform // with requested parameters. Call to the static function SineCosine returns a // pointer to a sine or cosine transform with requested parameters. Example: // { // Int_t N = 10; Double_t *in = new Double_t[N]; // TVirtualFFT *fftr2c = TVirtualFFT::FFT(1, &N, "R2C"); // fftr2c->SetPoints(in); // fftr2c->Transform(); // Double_t re, im; // for (Int_t i=0; i<N; i++) // fftr2c->GetPointComplex(i, re, im); // ... // fftr2c->SetPoints(in2); // ... // fftr2c->SetPoints(in3); // ... // } // Different options are explained in the function comments // // // // // ////////////////////////////////////////////////////////////////////////// #include "TROOT.h" #include "TVirtualFFT.h" #include "TPluginManager.h" #include "TEnv.h" TVirtualFFT *TVirtualFFT::fgFFT = 0; TString TVirtualFFT::fgDefault = ""; ClassImp(TVirtualFFT) //_____________________________________________________________________________ TVirtualFFT::~TVirtualFFT() { //destructor if (this==fgFFT) fgFFT = 0; } //_____________________________________________________________________________ TVirtualFFT* TVirtualFFT::FFT(Int_t ndim, Int_t *n, Option_t *option) { //Returns a pointer to the FFT of requested size and type. //Parameters: // -ndim : number of transform dimensions // -n : sizes of each dimension (an array at least ndim long) // -option : consists of 2 parts - flag option and an option to create a new TVirtualFFT // 1) transform type option: // Available transform types are: // C2CForward, C2CBackward, C2R, R2C, R2HC, HC2R, DHT // see class description for details // 2) flag option: choosing how much time should be spent in planning the transform: // Possible options: // "ES" (from "estimate") - no time in preparing the transform, // but probably sub-optimal performance // "M" (from "measure") - some time spend in finding the optimal way // to do the transform // "P" (from "patient") - more time spend in finding the optimal way // to do the transform // "EX" (from "exhaustive") - the most optimal way is found // This option should be chosen depending on how many transforms of the // same size and type are going to be done. // Planning is only done once, for the first transform of this size and type. // 3) option allowing to choose between the global fgFFT and a new TVirtualFFT object // "" - default, changes and returns the global fgFFT variable // "K" (from "keep")- without touching the global fgFFT, // creates and returns a new TVirtualFFT*. User is then responsible for deleting it. // Examples of valid options: "R2C ES K", "C2CF M", "DHT P K", etc. Int_t inputtype=0, currenttype=0; TString opt = option; opt.ToUpper(); //find the tranform flag Option_t *flag; flag = "ES"; if (opt.Contains("ES")) flag = "ES"; if (opt.Contains("M")) flag = "M"; if (opt.Contains("P")) flag = "P"; if (opt.Contains("EX")) flag = "EX"; Int_t ndiff = 0; if (!opt.Contains("K")) { if (fgFFT){ //if the global transform exists, check if it should be changed if (fgFFT->GetNdim()!=ndim) ndiff++; else { Int_t *ncurrent = fgFFT->GetN(); for (Int_t i=0; i<ndim; i++){ if (n[i]!=ncurrent[i]) ndiff++; } } Option_t *t = fgFFT->GetType(); if (!opt.Contains(t)) { if (opt.Contains("HC") || opt.Contains("DHT")) inputtype = 1; if (strcmp(t,"R2HC")==0 || strcmp(t,"HC2R")==0 || strcmp(t,"DHT")==0) currenttype=1; if (!(inputtype==1 && currenttype==1)) ndiff++; } if (ndiff>0){ delete fgFFT; fgFFT = 0; } } } Int_t sign = 0; if (opt.Contains("C2CB") || opt.Contains("C2R")) sign = 1; if (opt.Contains("C2CF") || opt.Contains("R2C")) sign = -1; TVirtualFFT *fft = 0; if (opt.Contains("K") || !fgFFT) { TPluginHandler *h; TString pluginname; if (fgDefault.Length()==0) fgDefault="fftw"; if (strcmp(fgDefault.Data(),"fftw")==0) { if (opt.Contains("C2C")) pluginname = "fftwc2c"; if (opt.Contains("C2R")) pluginname = "fftwc2r"; if (opt.Contains("R2C")) pluginname = "fftwr2c"; if (opt.Contains("HC") || opt.Contains("DHT")) pluginname = "fftwr2r"; if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) { if (h->LoadPlugin()==-1) { printf("handler not found\n"); return 0; } fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE); Int_t *kind = new Int_t[1]; if (pluginname=="fftwr2r") { if (opt.Contains("R2HC")) kind[0] = 10; if (opt.Contains("HC2R")) kind[0] = 11; if (opt.Contains("DHT")) kind[0] = 12; } fft->Init(flag, sign, kind); if (!opt.Contains("K")) { fgFFT = fft; } delete [] kind; return fft; } else { printf("plugin not found\n"); return 0; } } } else { //if the global transform already exists and just needs to be reinitialised //with different parameters if (fgFFT->GetSign()!=sign || !opt.Contains(fgFFT->GetTransformFlag()) || !opt.Contains(fgFFT->GetType())) { Int_t *kind = new Int_t[1]; if (inputtype==1) { if (opt.Contains("R2HC")) kind[0] = 10; if (opt.Contains("HC2R")) kind[0] = 11; if (opt.Contains("DHT")) kind[0] = 12; } fgFFT->Init(flag, sign, kind); delete [] kind; } } return fgFFT; } //_____________________________________________________________________________ TVirtualFFT* TVirtualFFT::SineCosine(Int_t ndim, Int_t *n, Int_t *r2rkind, Option_t *option) { //Returns a pointer to a sine or cosine transform of requested size and kind // //Parameters: // -ndim : number of transform dimensions // -n : sizes of each dimension (an array at least ndim long) // -r2rkind : transform kind for each dimension // 4 different kinds of sine and cosine transforms are available // DCT-I - kind=0 // DCT-II - kind=1 // DCT-III - kind=2 // DCT-IV - kind=3 // DST-I - kind=4 // DST-II - kind=5 // DST-III - kind=6 // DST-IV - kind=7 // -option : consists of 2 parts - flag option and an option to create a new TVirtualFFT // - flag option: choosing how much time should be spent in planning the transform: // Possible options: // "ES" (from "estimate") - no time in preparing the transform, // but probably sub-optimal performance // "M" (from "measure") - some time spend in finding the optimal way // to do the transform // "P" (from "patient") - more time spend in finding the optimal way // to do the transform // "EX" (from "exhaustive") - the most optimal way is found // This option should be chosen depending on how many transforms of the // same size and type are going to be done. // Planning is only done once, for the first transform of this size and type. // - option allowing to choose between the global fgFFT and a new TVirtualFFT object // "" - default, changes and returns the global fgFFT variable // "K" (from "keep")- without touching the global fgFFT, // creates and returns a new TVirtualFFT*. User is then responsible for deleting it. // Examples of valid options: "ES K", "EX", etc TString opt = option; //find the tranform flag Option_t *flag; flag = "ES"; if (opt.Contains("ES")) flag = "ES"; if (opt.Contains("M")) flag = "M"; if (opt.Contains("P")) flag = "P"; if (opt.Contains("EX")) flag = "EX"; if (!opt.Contains("K")) { if (fgFFT){ Int_t ndiff = 0; if (fgFFT->GetNdim()!=ndim || strcmp(fgFFT->GetType(),"R2R")!=0) ndiff++; else { Int_t *ncurrent = fgFFT->GetN(); for (Int_t i=0; i<ndim; i++) { if (n[i] != ncurrent[i]) ndiff++; } } if (ndiff>0) { delete fgFFT; fgFFT = 0; } } } TVirtualFFT *fft = 0; if (!fgFFT || opt.Contains("K")) { TPluginHandler *h; TString pluginname; if (fgDefault.Length()==0) fgDefault="fftw"; if (strcmp(fgDefault.Data(),"fftw")==0) { pluginname = "fftwr2r"; if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) { if (h->LoadPlugin()==-1){ printf("handler not found\n"); return 0; } fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE); fft->Init(flag, 0, r2rkind); if (!opt.Contains("K")) fgFFT = fft; return fft; } else { printf("handler not found\n"); return 0; } } } //if (fgFFT->GetTransformFlag()!=flag) fgFFT->Init(flag,0, r2rkind); return fgFFT; } //_____________________________________________________________________________ TVirtualFFT* TVirtualFFT::GetCurrentTransform() { // static: return current fgFFT if (fgFFT) return fgFFT; else{ printf("fgFFT is not defined yet\n"); return 0; } } //_____________________________________________________________________________ void TVirtualFFT::SetTransform(TVirtualFFT* fft) { // static: set the current transfrom to parameter fgFFT = fft; } //_____________________________________________________________________________ const char *TVirtualFFT::GetDefaultFFT() { // static: return the name of the default fft return fgDefault.Data(); } //______________________________________________________________________________ void TVirtualFFT::SetDefaultFFT(const char *name) { // static: set name of default fft if (fgDefault == name) return; delete fgFFT; fgFFT = 0; fgDefault = name; }