QRealComplexFFTGSL.cc

Go to the documentation of this file.

/*
* 
* Class QRealComplexFFTGSL
*
* Computes the FFT of the double array data[] The result is stored in
* the array result, which must be twice the size N of data.  The
* storage arrangement of the transform is the following:
* the first N elements of result are the real elememants of the FFT, while 
* the remaining ones are the imaginary ones
* 
* This algorithm uses GSL libraries to perform the FFT.
*/
 
#include "QRealComplexFFTGSL.hh"
#include <stdio.h>
#include <iostream>
#include "QVectorC.hh"
#include <gsl/gsl_fft_halfcomplex.h>
 
ClassImp(Diana::QRealComplexFFTGSL);
 
Q_BEGIN_NAMESPACE
 
using namespace std;
 
QRealComplexFFTGSL::QRealComplexFFTGSL(const int size):QFFT(size)
 {
    fWaveTable = gsl_fft_real_wavetable_alloc(size);
    fWork = gsl_fft_real_workspace_alloc(size);
    fHCWaveTable =  gsl_fft_halfcomplex_wavetable_alloc(size);  
    fWindow = NULL;
    fWindowType = WT_None;
    IsPowerOf2();
}
 
QRealComplexFFTGSL::QRealComplexFFTGSL() : QFFT()
{
    fWaveTable = 0;
    fWork = 0;
    fHCWaveTable = 0;
    fData = 0;
    fWindow = NULL;
    fWindowType = WT_None;
    
}
 
QRealComplexFFTGSL::~QRealComplexFFTGSL(){
  if(fWaveTable){
    gsl_fft_real_wavetable_free(fWaveTable);
  }
  if(fHCWaveTable){
    gsl_fft_halfcomplex_wavetable_free(fHCWaveTable);
  }
  if(fWork){
    gsl_fft_real_workspace_free(fWork);
  }
  if(fData){
    delete [] fData;
    fData = 0;
  }
  if(fWindow)
    delete [] fWindow;
  fWindow = 0;
}
 
int QRealComplexFFTGSL::TransformToFreq(const QVector& origData, QVectorC& FT,
                                        bool compress){ 
  Resize(origData.Size());
  SetForward(true);
  SetNormalized(false);
  QVector res(2*origData.Size());
  int r=Transform(origData,res);
  if(compress){
    FT.Resize(fSize/2+1);
    for(size_t i=0;i<=fSize/2;i++){
      FT[i].SetRe(res[i]);
      FT[i].SetIm(res[i+fSize]);
    }
  }
  else{
    FT.Resize(fSize);
    for(size_t i=0;i<fSize;i++){
      FT[i].SetRe(res[i]);
      FT[i].SetIm(res[i+fSize]);
    }
  }
  return r;
}
 
int QRealComplexFFTGSL::TransformFromFreq(const QVectorC& origFT,QVector& data,
                                          bool compress){ 
  if(compress)
    Resize(2*(origFT.Size()-1));
  else
    Resize(origFT.Size());
  SetForward(false);
  SetNormalized(true);
  QVectorC fullFT(fSize);
  int ret=0;
  if(compress){
    for(size_t i=0;i<fSize/2;i++){
      fullFT[i]=origFT[i];
      fullFT[fSize/2+i]=origFT[fSize/2-i].Conj();
    }
    ret = Transform(fullFT.SingleVector(),data);
  }
  else{
    ret = Transform(origFT.SingleVector(),data);
  }
  return ret;
}
 
int QRealComplexFFTGSL::Transform(const QVector& origData, QVector& result){
  if(fDirection==kForward)
    result.Resize(fSize*2);
  else 
    result.Resize(fSize);
 
  // readjust the window if not appropriate
  SetWindowType(fWindowType);
  
  if(fWindow && (fDirection==kForward)){
    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];
    } 
  }
  
  int rval = -1;
  if(fDirection == kForward && fNormalized == false){
    if (fIsPowerOf2){
      rval = gsl_fft_real_radix2_transform(fData,1,fSize);
      if(rval == GSL_SUCCESS){
        //real elements
        for (size_t i = 0; i <= fSize/2; i++){
          result[i] = fData[i];
        }
        for (size_t i = fSize/2 + 1; i < fSize; i++){
          result[i] = fData[fSize - i];
        }
        //first complex element
        size_t s32=size_t(1.5*fSize);
        result[fSize] = 0;result[s32]=0;
        for (size_t i = fSize + 1; i < s32; i++){
          result[i] = fData[2*fSize - i];
        }
        for (size_t i = s32+1; i < 2*fSize; i++){
          result[i] = -fData[i - fSize];
        }
        return 0;
      }
      //error
      else if(rval != GSL_EDOM)
        return -1;
    }
    //fData size is not power of 2
    rval = gsl_fft_real_transform(fData, 1, fSize, fWaveTable, fWork);
    if(rval == GSL_SUCCESS){
      if (fSize%2){
        printf("QRealComplexFFTGSL::Transform(forward): odd lenght data size"
               " (%zd) not implemented yet\n",fSize);
        return -1;
      } else {
        size_t imOff = fSize;
        result[0] = fData[0]; 
        result[imOff] = 0.;
        for(size_t i = 1; i < fSize/2;  i++) {
          //real elements
          result[fSize-i] = result[i] = fData[i*2-1];
          //imag elements
          result[imOff+fSize-i] = -(result[imOff+i] = fData[i*2]);
        }
        result[fSize/2] = fData[fSize-1];
        result[imOff + fSize/2] = 0.;
        return 0;
      }
    }
    //error
    else if(rval != GSL_EDOM)
      return -1;
    
  }
  
  if(fDirection == kBackward){
    if(fIsPowerOf2) {
      // fData size is power of 2
      fData[0]=origData[0];
      fData[fSize/2]=origData[fSize/2];
      for (size_t i=1;i<fSize/2;i++){
        fData[i]=origData[i];
        fData[fSize-i]=origData[fSize+i];
      }
      int rval = fNormalized ? 
        gsl_fft_halfcomplex_radix2_inverse(fData,1,fSize) :
        gsl_fft_halfcomplex_radix2_backward(fData,1,fSize);
      if(rval == GSL_SUCCESS)result.SetArray(fData,fSize);
      return rval;
    } else {
      if (fSize%2){
        //fData size not multiple of 2
        printf("QRealComplexFFTGSL::Transform(backward): odd lenght data size"
               " (%zd) not implemented yet\n",fSize);
        return -1;
      } else {
        //fData size is multiple of 2
        fData[0] =  origData[0];
        fData[fSize-1] = origData[fSize/2];
        for(size_t i = 1; i < fSize/2; i++) {
          // real
          fData[i*2-1] = origData[i];
          //imag
          fData[i*2] = origData[fSize + i];
        }
        
        int rval = 
          gsl_fft_halfcomplex_inverse(fData, 1, fSize,fHCWaveTable, fWork);
        if(rval == GSL_SUCCESS)result.SetArray(fData,fSize);
        return rval;
      }
      
    }
  }
  return rval;
}
bool QRealComplexFFTGSL::Resize(size_t n){
  if (n != fSize){
    if(fWaveTable)
      gsl_fft_real_wavetable_free(fWaveTable);
    if(fHCWaveTable)
      gsl_fft_halfcomplex_wavetable_free(fHCWaveTable);
    if(fWork)
      gsl_fft_real_workspace_free(fWork);
    if(fData)
      delete [] fData;
    
    fData = new double[n]; 
    fWaveTable = gsl_fft_real_wavetable_alloc(n);
    fWork = gsl_fft_real_workspace_alloc(n);
    fHCWaveTable =  gsl_fft_halfcomplex_wavetable_alloc(n);  
    fSize = n;
    IsPowerOf2();
    if(fWindow){
      SetWindowType(fWindowType);
    }
  }
  return true;
}
 
void QRealComplexFFTGSL::SetWindowType(WindowType w, int coherent){
  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;
  }
}
 
Q_END_NAMESPACE