QRealComplexFFTW3.cc

Go to the documentation of this file.

#include "QRealComplexFFTW3.hh"
#include <stdio.h>
#include <iostream>
#include "QVectorC.hh"
 
ClassImp(Diana::QRealComplexFFTW3);
 
Q_BEGIN_NAMESPACE
 
using namespace std;
QRealComplexFFTW3::QRealComplexFFTW3(const size_t size)
  : QFFT(), // <Call the empty ctor, do our own allocations
    fStride(1) // We are only doing 1 dim transforms so stride=1
{
#ifdef _HAVE_FFTW3_
  fSizeLock=false;
  fDirection = kForward;
  fNormalized = false;
  // But leave the plans null for now (they will be created when needed)
  fFFTW_Forward_Plan=NULL;
  fFFTW_Backward_Plan=NULL;
  fWindow = NULL;
  fWindowType = WT_None;
  fSize = 0;
  fAllocatedSize = 0;
  Resize(size);
#else
  throw QError(QERR_UNKNOWN_ERR,__FILE__,__LINE__,
               "Package mathtools not compiled with FFTW3 support!");
#endif
}
QRealComplexFFTW3::QRealComplexFFTW3() 
  : QFFT(),
    fStride(1)
{
#ifdef _HAVE_FFTW3_
  fSizeLock=false;
  fFFTW_Forward_Plan=NULL;
  fFFTW_Backward_Plan=NULL;
  fData = 0;
  fComplexData=0;
  fSize = 0;
  fAllocatedSize = 0;
  fWindow = NULL;
  fWindowType = WT_None;
#else
  throw QError(QERR_UNKNOWN_ERR,__FILE__,__LINE__,
               "Package mathtools not compiled with FFTW3 support!");
#endif  
}
 
QRealComplexFFTW3::~QRealComplexFFTW3(){
  ClearCache();
}
 
int QRealComplexFFTW3::TransformToFreq(const QVector& origData, QVectorC& FT,
                                       bool compress){ 
#ifdef _HAVE_FFTW3_
  if(!Resize(origData.Size()))
    return -1;
  if(compress)
    FT.Resize(fSize/2+1);
  else
    FT.Resize(fSize);
  
  SetForward(true);
 
  // re adjust the window if not appropriate
  SetWindowType(fWindowType);
 
  if(fFFTW_Forward_Plan==NULL)
    if(!CreatePlan(true))
      return -1;
 
  // If we are windowing
  if(fWindow){
    for(size_t i=0;i<fSize;i++)
      fData[i]=origData[i]*fWindow[i];
  }
  else{
    for(size_t i=0;i<fSize;i++)
      fData[i]=origData[i];
  }
    
  // Execute the plan
  fftw_execute(fFFTW_Forward_Plan);
 
  // Now we have done the FFT, we should move the data back if needed
  QVector &Re=FT.Re();
  QVector &Im=FT.Im();
  // Do the lower half
  for(size_t i=0;i<fSize/2+1;i++){
    Re[i]=fComplexData[i][0];
    Im[i]=fComplexData[i][1];
  }
  if(!compress){
    // Now we should fill out the complex conjugates
    for(size_t i=0;i<fSize/2-1;i++){
      Re[fSize/2+i+1]=Re[fSize/2-i-1];
      Im[fSize/2+i+1]=-Im[fSize/2-i-1];
    }
  }
#else
  throw QError(QERR_UNKNOWN_ERR,__FILE__,__LINE__,
               "Package mathtools not compiled with FFTW3 support!");
#endif  
  return 0;
}
 
int QRealComplexFFTW3::TransformFromFreq(const QVectorC& origFT, QVector& data,
                                         bool compress){ 
#ifdef _HAVE_FFTW3_
  int ret;
  if(compress)
    ret = Resize(2*(origFT.Size()-1));
  else
    ret = Resize(origFT.Size());
 
  if(!ret)
    return -1;
 
  data.Resize(fSize);
 
  SetForward(false);
 
  // re adjust the window if not appropriate
  SetWindowType(fWindowType);
 
  // Make sure a plan is in place
  if(fFFTW_Backward_Plan==NULL)
    if(!CreatePlan(false))
      return -1;
 
  // Backward transforms
  // Copy the data from the source
  // Ignore all of the data in the second half of the vector.
  const QVector &Re=origFT.Re();
  const QVector &Im=origFT.Im();
  for(size_t i=0;i<fSize/2+1;i++){
    fComplexData[i][0]=Re[i];
    fComplexData[i][1]=Im[i];
  }
  // Execute the plan
  fftw_execute(fFFTW_Backward_Plan);
  
  if(fWindow){
    for(size_t i=0;i<fSize;i++)
      data[i]=fData[i]/fWindow[i];
  }
  else{
    for(size_t i=0;i<fSize;i++)
      data[i]=fData[i];
  }
  data/=fSize;  
#else
  throw QError(QERR_UNKNOWN_ERR,__FILE__,__LINE__,
               "Package mathtools not compiled with FFTW3 support!");
#endif  
  return 0;
}
 
int QRealComplexFFTW3::Transform(const QVector& origData, QVector& result){
  return 0;
}
bool QRealComplexFFTW3::CreatePlan(bool IsForward){
#ifdef _HAVE_FFTW3_
  if(IsForward){
    fFFTW_Forward_Plan=fftw_plan_dft_r2c_1d(fSize,fData,fComplexData,
                                            FFTW_MEASURE);
    if(fFFTW_Forward_Plan==NULL)
      return false;
  }
  else {
    fFFTW_Backward_Plan=fftw_plan_dft_c2r_1d(fSize,fComplexData,fData,
                                             FFTW_MEASURE);
    if(fFFTW_Backward_Plan==NULL)
      return false;
  }
#endif
  return true;
}
/***********************************************************************
  This code requires FFTW >3.3.4 to work. Currently, by default it
  will return false.
bool QRealComplexFFTW3::CheckPlanCompatible(const QVector& data,
                                            const QVector& result) const {
 
  // Check if we can execute plan without moving data
  bool externMem=true;
 
  double *real=(fDirection==kForward ? data.fData : result.fData);
  double *comp=(fDirection==kForward ? result.fData : data.fData);
 
  // We cant be windowing
  externMem &= !fWindow;
  // The strides must be the same (this is trivially true here).
  externMem &= (data.fStride==fStride && result.fStride==fStride);
  // The alignments must be the same
  externMem &= 
    (fftw_alignment_of(real) == fftw_alignment_of(fData));
  externMem &= 
    (fftw_alignment_of(comp) == fftw_alignment_of(fComplexData));
  
  return externMem;
}
************************************************************************/
bool QRealComplexFFTW3::Resize(size_t size){ 
#ifdef _HAVE_FFTW3_
  if(fSize==size) return true;
  else if(fSizeLock) return false;
  else if(size < fAllocatedSize){
    // Clear the cache, but leave the memory
    ClearCache(false);
    // Set a new size
    fSize = size;
  }
  else {
    // Clear the cache and the allocated memory
    ClearCache(true);
    // Set a new size
    fSize=size;
    // Set the new allocate size
    fAllocatedSize=size;
    // Allocate the memory for the real data
    fData = QVector::ArrayAlloc(size);
    // Allocate the memory for the complex data.  Since this is a real
    // FFT, the resulting array only needs to hold half of the complex
    // values, but the data are not 'packed', and so we require only n+1
    // (n is odd) or n+2 (n is even) double elements. (Or n/2+1 complex
    // elements)
    fComplexData=
      (fftw_complex *)QVector::ArrayAlloc(size+((size % 2) ? 1 : 2));
 
  }
#endif
  return true;
}
bool QRealComplexFFTW3::ClearCache(bool clearMem){
#ifdef _HAVE_FFTW3_
  // Delete the allocated memory but don't delete the FFTW3
  // wisdom. (But dont save it either.)
  if(clearMem && fData){
    QVector::ArrayFree(fData);
    QVector::ArrayFree((double *)fComplexData);
    fData=NULL;
    fComplexData=NULL;
  }
  if(fFFTW_Forward_Plan){
    fftw_destroy_plan(fFFTW_Forward_Plan);
    fFFTW_Forward_Plan=NULL;
  }
  if(fFFTW_Backward_Plan){
    fftw_destroy_plan(fFFTW_Backward_Plan);
    fFFTW_Backward_Plan=NULL;
  }
#endif
  return true;
}
void QRealComplexFFTW3::SetWindowType(WindowType w, int coherent){
#ifdef _HAVE_FFTW3_
  if(w != fWindowType || fWindowSize != fSize){
    QVector tmp  = GetWindow(w,fSize,0,coherent);
    size_t wsize = tmp.Size();
    if(wsize > 0 ) {
      fWindow = tmp.GetArray();
      fWindowSize = fSize;
    } else { 
      fWindow = NULL;
      fWindowSize = 0;
    }
    fWindowType = w;
  }
#endif
}
 
// inside this else statement is code if we dont have the fftw3
// headers
 
 
Q_END_NAMESPACE