Index: trunk/MagicSoft/Mars/manalysis/MFFT.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFFT.cc (revision 2807) +++ (revision ) @@ -1,752 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Markus Gaug 01/2004 -! -! Copyright: MAGIC Software Development, 2001-2004 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// Fast Fourier Transforms // -// // -// (Most of the code is adapted from Numerical Recipies in C++, 2nd ed., // -// pp. 509-563) // -// // -// Usage: // -// // -// 1) Functions RealFunctionFFT: (FOURIER TRANSFORM) // -// * Take as argument arrays of real numbers, // -// in some cases the dimension of the array has to be given separately// -// * Return a COMPLEX array with the following meaning: // -// array[0]: The value of F(0) (has only real component) // -// array[1]: The value of F(N/2) (has only real component) // -// array[2i]: The real part of F(i) // -// array[2i+1]: The imaginary part of F(i) // -// * Note that F(N-i)* = F(i), therefore only the positive frequency // -// half is stored. // -// * The dimension MUST be an integer power of 2, // -// otherwise, the array will be shortened!! // -// // -// 2) Functions RealFunctionIFFT: (INVERSER FOURIER TRANSFORM) // -// * Take as argument a COMPLEX array // -// of Fourier-transformed REAL numbers // -// with the following meaning: // -// array[0]: The value of F(0) (has only real component) // -// array[1]: The value of F(N/2) (has only real component) // -// array[2i]: The real part of F(i) // -// array[2i+1]: The imaginary part of F(i) // -// * Returns the original complex array of dimension 2N-1 // -// // -// 3) Functions PowerSpectrumDensity: // -// * Return a histogram with the spectral density, i.e. // -// P(k) = 1/(N*N) * |F(k)|*|F(k)| // -// * The histogram is ranged between 0 and 1./(2*binwidth) // -// * The number of bins equals N/2+1 // -// * Note that histograms with unequal binwidth can not yet be treated! // -// * If the PSD does NOT CONVERGE to 0 at the maximum bin, // -// you HAVE TO sample your data finer! // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MFFT.h" - -#include "TMath.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MFFT); - -using namespace std; - -// --------------------------------------------------------------------------- -// -// Default Constructor -// Initializes random number generator and default variables -// -MFFT::MFFT() : fDim(0) -{ -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MFFT::~MFFT() -{ -} - -void MFFT::TransformF(const Int_t isign) -{ - - UInt_t n,mmax,m,j,istep,i; - Float_t wtemp,wr,wpr,wpi,wi,theta; - Float_t tempr,tempi; - - Int_t nn = fDim/2; - n = nn << 1; - - // - // The bit-reversal section of the routine: - // Exchange the two complex numbers - // - j=1; - for (i=1;i i) { - Swap(fDataF[j-1],fDataF[i-1]); - Swap(fDataF[j],fDataF[i]); - } - m=nn; - while (m >= 2 && j > m) { - j -= m; - m >>= 1; - } - j += m; - } - // - // Here begins the Danielson-Lanczos section of the routine - // - mmax=2; - while (n > mmax) { // Outer loop executed log_2(nn) times - - istep = mmax << 1; - // - // Initialize the trigonometric recurrence: - // - theta = isign*(6.28318530717959/mmax); - - wtemp = TMath::Sin(0.5*theta); - wpr = -2.0*wtemp*wtemp; - wpi = TMath::Sin(theta); - - wr=1.0; - wi=0.0; - - for (m=1; m i) { - Swap(fDataD[j-1],fDataD[i-1]); - Swap(fDataD[j],fDataD[i]); - } - m=nn; - while (m >= 2 && j > m) { - j -= m; - m >>= 1; - } - j += m; - } - // - // Here begins the Danielson-Lanczos section of the routine - // - mmax=2; - while (n > mmax) { // Outer loop executed log_2(nn) times - - istep = mmax << 1; - // - // Initialize the trigonometric recurrence: - // - theta = isign*(6.28318530717959/mmax); - - wtemp = TMath::Sin(0.5*theta); - wpr = -2.0*wtemp*wtemp; - wpi = TMath::Sin(theta); - - wr=1.0; - wi=0.0; - - for (m=1; m>1); - - if(isign==1) // forward transform - { - c2 = -0.5; - TransformF(1); - } - else // set up backward transform - { - c2 = 0.5; - theta = -theta; - } - - wtemp = TMath::Sin(0.5*theta); - wpr = -2.0*wtemp*wtemp; - wpi = TMath::Sin(theta); - - wr = 1.0 + wpr; - wi = wpi; - - for(i=1;i<(fDim>>2);i++) // case i=0 done separately below - { - - i2 = 1 + (i1 = i+i); - i4 = 1 + (i3 = fDim-i1); - - // - // The two separate transforms are separated out of the data - // - h1r = c1*(fDataF[i1]+fDataF[i3]); - h1i = c1*(fDataF[i2]-fDataF[i4]); - h2r = -c2*(fDataF[i2]+fDataF[i4]); - h2i = c2*(fDataF[i1]-fDataF[i3]); - - // - // Here, they are recombined to from the true transform - // of the orginal real data - // - fDataF[i1] = h1r + wr*h2r - wi*h2i; - fDataF[i2] = h1i + wr*h2i + wi*h2r; - fDataF[i3] = h1r - wr*h2r + wi*h2i; - fDataF[i4] = -h1i + wr*h2i + wi*h2r; - - // - // The recurrence - // - wr = (wtemp=wr)*wpr - wi*wpi + wr; - wi = wi*wpr + wtemp*wpi + wi; - } - - // - // Squeeze the first and last data together to get them all - // within the original array - // - if(isign==1) - { - fDataF[0] = (h1r=fDataF[0]) + fDataF[1]; - fDataF[1] = h1r - fDataF[1]; - } - else - { - - fDataF[0] = c1*((h1r=fDataF[0]) + fDataF[1]); - fDataF[1] = c1*(h1r-fDataF[1]); - - // - // The inverse transform for the case isign = -1 - // - TransformF(-1); - - // - // normalize correctly (not done in original NR's) - // - for(i=1;i<=fDim;i++) - fDataF[i] *= (2./fDim); - } -} -void MFFT::RealFTD(const Int_t isign) -{ - - Int_t i,i1,i2,i3,i4; - Float_t c1=0.5,c2,h1r,h1i,h2r,h2i; - Double_t wr,wi,wpr,wpi,wtemp,theta; - - // - // Initialize the recurrence - // - theta=3.141592653589793/(Double_t) (fDim>>1); - - if(isign==1) // forward transform - { - c2 = -0.5; - TransformD(1); - } - else // set up backward transform - { - c2 = 0.5; - theta = -theta; - } - - wtemp = TMath::Sin(0.5*theta); - wpr = -2.0*wtemp*wtemp; - wpi = TMath::Sin(theta); - - wr = 1.0 + wpr; - wi = wpi; - - for(i=1;i<(fDim>>2);i++) // case i=0 done separately below - { - - i2 = 1 + (i1 = i+i); - i4 = 1 + (i3 = fDim-i1); - - // - // The two separate transforms are separated out of the data - // - h1r = c1*(fDataD[i1]+fDataD[i3]); - h1i = c1*(fDataD[i2]-fDataD[i4]); - h2r = -c2*(fDataD[i2]+fDataD[i4]); - h2i = c2*(fDataD[i1]-fDataD[i3]); - - // - // Here, they are recombined to from the true transform - // of the orginal real data - // - fDataD[i1] = h1r + wr*h2r - wi*h2i; - fDataD[i2] = h1i + wr*h2i + wi*h2r; - fDataD[i3] = h1r - wr*h2r + wi*h2i; - fDataD[i4] = -h1i + wr*h2i + wi*h2r; - - // - // The recurrence - // - wr = (wtemp=wr)*wpr - wi*wpi + wr; - wi = wi*wpr + wtemp*wpi + wi; - } - - // - // Squeeze the first and last data together to get them all - // within the original array - // - if(isign==1) - { - fDataD[0] = (h1r=fDataD[0]) + fDataD[1]; - fDataD[1] = h1r - fDataD[1]; - } - else - { - - fDataD[0] = c1*((h1r=fDataD[0]) + fDataD[1]); - fDataD[1] = c1*(h1r-fDataD[1]); - - // - // The inverse transform for the case isign = -1 - // - TransformD(-1); - - // - // normalize correctly (not done in original NR's) - // - for(i=1;i<=fDim;i++) - fDataD[i] *= (2./fDim); - } -} - - -// -// Fast Fourier Transform for float arrays -// -Float_t* MFFT::RealFunctionFFT(const Int_t n, const Float_t *data) -{ - - fDim = n; - CheckDim(n); - - fDataF.Set(fDim); - // - // Clone the array - // - for (Int_t i=0;iGetSize(); - CheckDim(fDim); - - fDataF.Set(fDim); - // - // Clone the array - // - for (Int_t i=0;iAt(i); - - RealFTF(1); - - return new TArrayF(fDim,fDataF.GetArray()); - -} - -// -// Inverse Fast Fourier Transform for TArrayF's -// -TArrayF* MFFT::RealFunctionIFFT(const TArrayF *data) -{ - - fDim = data->GetSize(); - CheckDim(fDim); - - fDataF.Set(fDim); - // - // Clone the array - // - for (Int_t i=0;iAt(i); - - RealFTF(-1); - - return new TArrayF(fDim,fDataF.GetArray()); -} - - -// -// Fast Fourier Transform for TArrayD's -// -TArrayD* MFFT::RealFunctionFFT(const TArrayD *data) -{ - - fDim = data->GetSize(); - CheckDim(fDim); - - fDataD.Set(fDim); - // - // Clone the array - // - for (Int_t i=0;iAt(i); - - RealFTD(1); - - return new TArrayD(fDim,fDataD.GetArray()); - -} - -// -// Inverse Fast Fourier Transform for TArrayD's -// -TArrayD* MFFT::RealFunctionIFFT(const TArrayD *data) -{ - - fDim = data->GetSize(); - CheckDim(fDim); - - fDataD.Set(fDim); - // - // Clone the array - // - for (Int_t i=0;iAt(i); - - RealFTD(-1); - - return new TArrayD(fDim,fDataD.GetArray()); -} - - -// -// Power Spectrum Density Calculation -// -TH1D* MFFT::PowerSpectrumDensity(const TH1D *hist) -{ - - TH1D *newhist = (TH1D*)CheckHist(hist,1); - - fDataD.Set(fDim); - // - // Copy the hist into an array - // - for (Int_t i=0;iGetBinContent(i); - - RealFTD(1); - - Int_t dim2 = fDim*fDim; - Double_t c02; - Double_t ck2; - Double_t cn2; - // - // Fill the new histogram: - // - // 1) P(0) = 1/(N*N) |C(0)|*|C(0)| - // (stored in fData{0]) - // - c02 = fDataD[0]*fDataD[0]; - newhist->Fill(c02/dim2); - // - // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)| + |C(N-k)|*|C(N-k)|) - // - for (Int_t k=2;kFill(ck2/dim2); - } - // - // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|) - // (stored in fData[1]) - // - cn2 = (fDataD[1]*fDataD[1]); - newhist->Fill(cn2/dim2); - - return newhist; -} - - -// -// Power Spectrum Density calculation -// -TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist) -{ - - TH1F *newhist = (TH1F*)CheckHist(hist,0); - - fDataF.Set(fDim); - // - // Copy the hist into an array - // - for (Int_t i=0;iGetBinContent(i); - - RealFTF(1); - - Int_t dim2 = fDim*fDim; - Float_t c02; - Float_t ck2; - Float_t cn2; - // - // Fill the new histogram: - // - // 1) P(0) = 1/(N*N) |C(0)|*|C(0)| - // - c02 = (fDataF[0]*fDataF[0]); - newhist->Fill(c02/dim2); - // - // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|)) - // - for (Int_t k=2;kFill(ck2/dim2); - } - // - // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|) - // - cn2 = (fDataF[1]*fDataF[1]); - newhist->Fill(cn2/dim2); - - return newhist; -} - - -void MFFT::CheckDim(Int_t a) -{ - - // If even number, return 0 - if (a==2) return; - - // If odd number, return the closest power of 2 - if (a & 1) - { - Int_t b = 1; - while (b < fDim/2+1) - b <<= 1; - - fDim = b; - *fLog << warn << "Dimension of Data is not a multiple of 2, will take only first " - << fDim << " entries! " << endl; - return; - } - - CheckDim(a>>1); -} - -TH1* MFFT::CheckHist(const TH1 *hist, const Int_t flag) -{ - - TString name = hist->GetName(); - name += " PSD"; - TString title = hist->GetTitle(); - title += " - Power Spectrum Density"; - - // number of entries - fDim = hist->GetNbinsX(); - CheckDim(fDim); - - // Step width - Double_t delta = hist->GetBinWidth(1); - - // Nyquist frequency - Axis_t fcrit = 1./(2.*delta); - Axis_t low = 0.; - Axis_t up = fcrit; - - switch (flag) - { - case 0: - return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up); - break; - case 1: - return new TH1D(name.Data(),title.Data(),fDim/2+1,low,up); - break; - default: - return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up); - break; - } -} Index: trunk/MagicSoft/Mars/manalysis/MFFT.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFFT.h (revision 2807) +++ (revision ) @@ -1,66 +1,0 @@ -#ifndef MARS_MFFT -#define MARS_MFFT - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -#ifndef ROOT_TArrayF -#include "TArrayF.h" -#endif - -#ifndef ROOT_TArrayD -#include "TArrayD.h" -#endif - -#ifndef ROOT_TH1F -#include "TH1F.h" -#endif - -#ifndef ROOT_TH1D -#include "TH1D.h" -#endif - -class MFFT : public MParContainer -{ -private: - - void Swap(Float_t &a, Float_t &b) { Float_t c = a; a = b; b = c; } - void Swap(Double_t &a, Double_t &b) { Double_t c = a; a = b; b = c; } - - void TransformF(const Int_t isign); - void TransformD(const Int_t isign); - void RealFTF(const Int_t isign); - void RealFTD(const Int_t isign); - - void CheckDim(Int_t a); - TH1 *CheckHist(const TH1 *hist, const Int_t flag); - - Int_t fDim; - TArrayF fDataF; - TArrayD fDataD; - -public: - - MFFT(); - ~MFFT(); - - TArrayF* RealFunctionFFT( const TArrayF *data); - TArrayF* RealFunctionIFFT(const TArrayF *data); - - TArrayD* RealFunctionFFT( const TArrayD *data); - TArrayD* RealFunctionIFFT(const TArrayD *data); - - Float_t* RealFunctionFFT( const Int_t n, const Float_t *data); - Float_t* RealFunctionIFFT(const Int_t n, const Float_t *data); - - Double_t* RealFunctionFFT( const Int_t n, const Double_t *data); - Double_t* RealFunctionIFFT(const Int_t n, const Double_t *data); - - TH1F* PowerSpectrumDensity(const TH1F *hist); - TH1D* PowerSpectrumDensity(const TH1D *hist); - - ClassDef(MFFT,0) // Class to perform a Fast Fourier Transform -}; - -#endif