QFFT.cc

Go to the documentation of this file.

/*
* 
* Class QFFT
*
*/
 
#include "QFFT.hh"
#include "QVectorView.hh"
 
ClassImp(Diana::QFFT);
 
Q_BEGIN_NAMESPACE
 
QFFT::QFFT(const size_t size):fSize(size) {
    fDirection = kForward;
    fNormalized = false;
    IsPowerOf2();
    fData = new double[size]; 
    fWindow = 0;
    fWindowSize = 0;
    fWindowType = WT_None;
}
 
QFFT::QFFT(){
    fSize = 0;
    fDirection = kForward;
    fNormalized = false;
    fIsPowerOf2 = false;
    fData = NULL;
    fWindow = 0;
    fWindowSize = 0;
    fWindowType = WT_None;
}
 
QFFT::~QFFT(){
    if(fData){
        delete [] fData;
        fData = 0;
    }
    if(fWindow){
        delete [] fWindow;
        fWindow = 0;
    }
} 
bool QFFT::IsPowerOf2(UInt_t Value,UInt_t &Power,UInt_t &Remainder){
  Power=Remainder=0;
  if(Value == 0) { Power=0;Remainder=0;return false; }
  size_t Comp=2*sizeof(size_t);
  while(Comp>0){
    if(Value>(size_t)((1<<Comp)-1))
      {Power+=Comp; Remainder+=(Value&((1<<Comp)-1)); Value = Value>>Comp; }
    Comp/=2;
  }
  if(Remainder) return false;
  return true;
}
bool QFFT::IsPowerOf2(UInt_t Value){
  UInt_t P,R;
  return IsPowerOf2(Value,P,R);
}
void QFFT::IsPowerOf2(){
  fIsPowerOf2=IsPowerOf2(fSize);
}
QFFT::WindowType QFFT::StrToWindowType(const std::string& winName)
{
 
  WindowType wType;
 
  if     (!winName.compare("rectangular"))      wType=WT_Rect;     
  else if(!winName.compare("welch"))            wType=WT_Welch;             
  else if(!winName.compare("hann"))             wType=WT_Hann;             
  else if(!winName.compare("hamming"))          wType=WT_Hamming;            
  else if(!winName.compare("cosinus"))          wType=WT_Cosinus;     
  else if(!winName.compare("blackman-harris"))  wType=WT_BlackmanHarris; 
  else if(!winName.compare("kaiser3"))          wType=WT_Kaiser3; 
  else wType = WT_None;
  
  return wType;
 
}
 
QVector QFFT::GetWindow(WindowType wt, size_t size, size_t zeros, int coherent, int param )
{
 
    size_t winSize = size-zeros;
    QVector window(winSize);
 
    switch(wt) {
        case WT_Rect:
            {
                for(size_t i = 0; i < winSize; i++) {
                    window[i] = 1.;
                }
                break;
            }
        case WT_Welch:
            {
                for(size_t i = 0; i < winSize; i++) {
                    window[i] = 1.-  pow( (i-winSize/2.) / (winSize/2.),2.);
                }
                break;
            }
        case WT_Hann: 
            {
                for(size_t i = 0; i < winSize; i++) {
                    window[i] = 0.5*(1 - gsl_sf_cos(2*M_PI*(i+1)/winSize));
                }
                break;
            }
        case WT_Hamming: 
            {
                for(size_t i = 0; i < winSize; i++) {
                    window[i] = 0.53836 - 0.46164*gsl_sf_cos(2*M_PI*(i+1)/(winSize));
                }
                break;
            }
        case WT_Cosinus: 
            {
                size_t risetime = 30;
                if(param > 0) risetime = param;
                for(size_t i = 0; i < winSize; i++) {
                    if(i< risetime) 
                        window[i] = 0.5 - 0.5*cos(M_PI/(risetime)*double(i+1));
                    else if(i < winSize-risetime) 
                        window[i] = 1.;
                    else 
                        window[i] = 0.5 + 0.5*cos(M_PI/(risetime)*double(i+1-(winSize-risetime)));
                }
                break;
            }
        case WT_BlackmanHarris:
            {
                double a0 = 0.35875;
                double a1 = 0.48829;
                double a2 = 0.14128;
                double a3 = 0.01168;
                for(size_t i = 0; i < winSize; i++) {
                    window[i] = a0 - a1*cos(2*M_PI*i/winSize) + a2*cos(4*M_PI*i/winSize) - a3*cos(6*M_PI*i/winSize);
                }
                break;
            }
        case WT_Kaiser3:
            {
                double arg;
                double alpha = 3.;
                for(size_t i = 0; i < winSize; i++) {
                    arg = alpha*sqrt(1.-pow(2.*i/winSize-1.,2.));
                    window[i] = gsl_sf_bessel_I0(arg)/gsl_sf_bessel_I0(alpha);
                }
                break;
            }
        case WT_None:
        default:
            {
                window.Resize(0);
                break;
            }
 
    }
    double corrFactor = 0.;
    if(coherent == 1) {
        for(size_t i = 0; i < window.Size(); i++) {
            corrFactor += window[i];
        }
        corrFactor /= window.Size();
    }  else if(coherent == 2){
        for(size_t i = 0; i < window.Size(); i++) {
            corrFactor += window[i]*window[i];
        }
        corrFactor /= window.Size();
        corrFactor = sqrt(corrFactor);
    } else {
        corrFactor = 1.;
    }
    for(size_t i = 0; i < window.Size(); i++) {
        window[i] /= corrFactor;
    }
 
 
    QVector retwin;
    if(zeros > 0 && window.Size() > 0) {
        retwin.Resize(size);
        // fill left and right side with zeros
        size_t start = zeros/2;
        size_t stop = start + window.Size();
        for(size_t j = 0; j < start; j++) {
            retwin[j] = 0.;
        }
        for(size_t j = start, i = 0; j < stop ; j++,i++) {
            retwin[j] = window[i];
        }
        for(size_t j = stop; j < size; j++) {
            retwin[j] = 0.;
        }
    } else {
        retwin = window;
    }
    return retwin;
}
 
QVector QFFT::FFTAntiSym(const QVector& input)
{
  Diana::QVector antiSym(input);
  size_t size=antiSym.Size()/2-1;
  antiSym.Shift(-size);
  return antiSym;
}
QVector QFFT::FFTSym(const QVector& input)
{
  Diana::QVector sym(input);
  size_t size=sym.Size()/2-1;
  sym.Shift(size);
  return sym;
}
 
QVector QFFT::ZeroPad(const QVector& input,size_t n_zeros,int isSym,
                      double zeroVal){
    // Calculate the new size of the array
    size_t totsize = input.Size() + n_zeros;
    // Offset for the input array
    size_t offset=(isSym == kMiddle ? input.Size()/2+1 : 0);
    // New array
    QVector zeroPadded(totsize);
 
    // Put the zeros at the front
    for(size_t i = 0; i < n_zeros; i++) {
      zeroPadded[i] = zeroVal;
    }
    // Put the input array at the end of the new array
    // If we are using kMiddle, the offset will tell this for loop 
    // to start in the middle of the input array
    for(size_t i=n_zeros,j=0;j<input.Size();i++,j++){
      zeroPadded[i]=input[(offset+j) % input.Size()];
    }
    // If we are symmetric, shift half the zeros to the back of the
    // array
    if(isSym==kSymmetric)
      zeroPadded.Shift(-(n_zeros/2));
    // If we are right padding, shift ALL the zeros to the back of the
    // array
    else if(isSym==kRight)
      zeroPadded.Shift(-n_zeros);
    // If we are left or middle padding, then shift the array by
    // offset, for kLeft this is zero, for kMiddle this is
    // input.Size()/2+1 (effectively undoing the offset that was done
    // before)
    else if(isSym==kMiddle)
      zeroPadded.Shift(offset);
    // Return the new array
    return zeroPadded;
}
 
QVectorC QFFT::ZeroPad(const QVectorC& input,size_t n_zeros,int isSym,
                       Diana::QComplex zeroVal){
    // Calculate the new size of the array
    size_t totsize = input.Size() + n_zeros;
    // Offset for the input array
    size_t offset=(isSym == kMiddle ? input.Size()/2+1 : 0);
    // New array
    QVectorC zeroPadded(totsize);
 
    // Put the zeros at the front
    for(size_t i = 0; i < n_zeros; i++) {
      zeroPadded[i] = zeroVal;
    }
    // Put the input array at the end of the new array
    // If we are using kMiddle, the offset will tell this for loop 
    // to start in the middle of the input array
    for(size_t i=n_zeros,j=0;j<input.Size();i++,j++){
      zeroPadded[i]=input[(offset+j) % input.Size()];
    }
    // If we are symmetric, shift half the zeros to the back of the
    // array
    if(isSym==kSymmetric)
      zeroPadded.Shift(-(n_zeros/2));
    // If we are right padding, shift ALL the zeros to the back of the
    // array
    else if(isSym==kRight)
      zeroPadded.Shift(-n_zeros);
    // If we are left or middle padding, then shift the array by
    // offset, for kLeft this is zero, for kMiddle this is
    // input.Size()/2+1 (effectively undoing the offset that was done
    // before)
    else if(isSym==kMiddle)
      zeroPadded.Shift(offset);
    // Return the new array
    return zeroPadded;
}
 
QVector QFFT::CutSides(const QVector& input, size_t ncut, bool isSym)
{
    size_t totsize = input.Size() - ncut;
    QVector theInput;
    if(isSym) {
        theInput = input;
    } else {
        theInput = FFTSym(input);
    }
    QVectorConstView subv(theInput,ncut/2,totsize);
 
    return FFTAntiSym(subv.GetVector());
}
 
Q_END_NAMESPACE