FancyFFTs.h

/* -*- C++ -*-.***************************************************************************************************
 Author: Ben Strutt
 Email: b.strutt.12@ucl.ac.uk

 Description: 
             C++ ROOT friendly namespace to do FFTs faster than Ryan. 
         Will probably be pretty bare bones intially.
         I only really want this for doing Cross Correlations.
*************************************************************************************************************** */

#ifndef FANCYFFTS_H
#define FANCYFFTS_H

#include <iostream>
#include <map>
#include <algorithm>
#include "TObject.h"
#include "TSystem.h"
#include "TMath.h"
#include "TGraph.h"
#include "TThread.h"


// Needs to be in front of including fftw3.h!
#ifdef __CINT__
#ifdef FFTW_64_BIT // Hack for Hawaii install of FFTW. Is there ever a reason to now have this defed? 
typedef struct {char a[16];} __float128; /* 16 chars have the same size as one __float128 */
#endif
#endif 


/* 
   Will use std::complex<double> for i/o as should be bit-to-bit identical to typdef fftw_complex double[2].
   So long as <complex> is included in front of fftw3.h.
*/
#include <complex>
#include <fftw3.h>



using std::complex;

namespace FancyFFTs {

  enum conventionFlag {
    kSum = 0,
    kAverage = 1,
    kTimeIntegral = 2,
    kPowSpecDensity = 3,
    kPowSpecDensity_dBm = 4
  };

  complex<double>* doFFT(int len, double* input, bool copyOutputToNewArray, int threadInd=0);
  complex<double>* doFFT(int len, double* input, complex<double>* output, int threadInd=0);
  complex<double>* doFFT(int len, double* input, complex<double>* output,
                     bool copyOutputToNewArray, int threadInd=0);  
  double* doInvFFT(int len, complex<double>* input, bool copyOutputToNewArray, int threadInd=0);
  double* doInvFFT(int len, complex<double>* input, double* output, int threadInd=0);  
  double* doInvFFT(int len, complex<double>* input, double* output,
              bool copyOutputToNewArray, int threadInd=0);

  double* getPowerSpectrum(int len, double* input, double dt, 
                  conventionFlag normFlag,
                  int threadInd=0);

  double* getPowerSpectrum(int len, double* input, double dt, 
                  conventionFlag normFlag,
                  double* outputPtr, int threadInd=0);

  TGraph* getPowerSpectrumTGraph(int len, double* input, double dt, 
                    conventionFlag normFlag,
                    bool dBScale, int threadInd=0);
  double* getFreqArray(int len, double dt);
  int getNumFreqs(int len);
  int printListOfKeys();
  double* crossCorrelate(int len, double* v1, double* v2, int threadInd=0);
  double* crossCorrelate(int len, complex<double>* fft1, complex<double>* fft2, int threadInd=0);
  double* crossCorrelate(int len, complex<double>* fft1, complex<double>* fft2,
                double* output, int threadInd=0);
  int extendToPowerOfTwo(int len);
  complex<double>* zeroPadFFT(complex<double>* fft, int numSamples, int numSamplesUpsampled);
  complex<double>* zeroPadFFT(complex<double>* fft, complex<double>* output,
                  int numSamples, int numSamplesUpsampled);  


 /* Takes care of checking whether a plan exists or not */
  bool makeNewPlanIfNeeded(int len,int threadInd=0);

  double* getRealArray(std::pair<Int_t, Int_t> key);
  

};


#endif //FANCYFFTS_H