BesselTransformer.cc

Go to the documentation of this file.

#include "BesselTransformer.hh"
#include "QRealComplexFFT.hh"
#include "QVectorView.hh"
#include <math.h>
 
 
BesselTransformer::BesselTransformer(double cutFreq, double duration, double samplFreq, double tauRC, bool zeroPad)
{
    fDuration = duration;
    fSize = (size_t) roundl(duration*samplFreq);
    fZeroPad = zeroPad;
 
    fSamplingFrequency = samplFreq;
    fCutFreq = cutFreq * duration;
    if(tauRC > 0.)  { 
        fRCFreq = 1./(tauRC*2*M_PI) * duration;
    } else fRCFreq = -1.;
 
    double norm = 2.703395061;
    fNormCutFreq = fCutFreq / norm;
  
    QVector besselSestoRe(fSize), besselSestoIm(fSize);
    QVector rcRe(fSize), rcIm(fSize);
    QVector transFunctionRe(fSize), transFunctionIm(fSize);
 
    for(size_t x = 0; x< fSize/2 + 1; x++) { 
        double xx = double(x) / fNormCutFreq;
        double valBesselRe = 10395. - 4725.*pow(xx,2) + 210.*pow(xx,4) - pow(xx,6);
        double valBesselIm =xx*( 10395.-1260.*pow(xx,2)+21.*pow(xx,4));
        besselSestoRe(x) = 10395. * valBesselRe/(pow(valBesselRe,2)+pow(valBesselIm,2));
        besselSestoIm(x) = -10395. * valBesselIm/(pow(valBesselRe,2)+pow(valBesselIm,2));
        if(fRCFreq > 0. )  {
            double xxrc = double(x)/fRCFreq;
            rcRe(x)  = 1./(1.+pow(xxrc,2));
            rcIm(x) = -xxrc/(1.+pow(xxrc,2));
        }
    }
 
    for(size_t y = 0; y<(fSize/2-1); y++) {
        double xx = - (fSize/2-1 - y);
        xx/=fNormCutFreq;
        double valBesselRe= 10395. - 4725.*pow(xx,2) + 210.*pow(xx,4) - pow(xx,6);
        double valBesselIm=xx*( 10395.-1260.*pow(xx,2)+21.*pow(xx,4));
        besselSestoRe(y+fSize/2+1) =  10395.*valBesselRe/(pow(valBesselRe,2)+pow(valBesselIm,2));
        besselSestoIm(y+fSize/2+1) = -10395.*valBesselIm/(pow(valBesselRe,2)+pow(valBesselIm,2));
        if(fRCFreq > 0.)  {
            double xxrc = double(- (fSize/2-1 - y))/fRCFreq;
            rcRe(y+fSize/2+1) = 1./(1.+pow(xxrc,2));
            rcIm(y+fSize/2+1) = -xxrc/(1.+pow(xxrc,2));
        }
    }
 
    for(size_t x = 0; x < fSize; x++) {
        if(fRCFreq > 0.) {
            transFunctionRe[x] = besselSestoRe[x]*rcRe[x] - besselSestoIm[x]*rcIm[x];
            transFunctionIm[x] = besselSestoRe[x]*rcIm[x] + besselSestoIm[x]*rcRe[x];
        } else {
            transFunctionIm[x] = besselSestoIm[x];
            transFunctionRe[x] = besselSestoRe[x];
        }
        //        transFunctionIm[x] = rcIm[x];
        //        transFunctionRe[x] = rcRe[x];
    }
 
    fTransferFunction.SetArray(transFunctionRe.GetConstArray(),transFunctionIm.GetConstArray(),fSize); 
    QRealComplexFFT transformer(fSize);
    transformer.TransformFromFreq(fTransferFunction,fDeltaResponse);
 
    // pad the second half with zeros to have good fft convolution
    if(fZeroPad) {
        size_t newsize = 2;
        while(newsize < fSize) {
            newsize *=2;
        }
        newsize = fSize*2;
        QVector deltaResponsez(newsize);
        QVectorView deltaResponsez_first(deltaResponsez,0,fSize);
        QVectorView deltaResponsez_second(deltaResponsez,fSize,newsize-fSize);
        deltaResponsez_first.GetVector() = fDeltaResponse;
        deltaResponsez_second.GetVector().Initialize(0.);
        fDeltaResponse = deltaResponsez;
        QRealComplexFFT transformer2(newsize);
        transformer2.TransformToFreq(fDeltaResponse,fTransferFunction);
    }
 
    // analyitical formula
    fAN = 3; 
    fAp = new double[fAN];
    fAq = new double[fAN];
    fAtheta = new double[fAN];
    fAphi = new double[fAN];
 
    fAtheta[0] = 4.248359395863364; 
    fAphi[0]   = 0.8675096732313656;
    fAp[0]     = 10.959228792517152;
    fAq[0]     = 39.425157113160175;
    fAtheta[1] = 3.735708356325815;
    fAphi[1]   = 2.6262723114471256;
    fAp[1]     = -14.126767991748853;
    fAq[1]     = -12.701164165560234;
    fAtheta[2] = 2.5159322478108215;
    fAphi[2]   = 4.492672953653942;
    fAp[2]     = 3.1675391992317308;
    fAq[2]     = 0.2024589840416812;
}
 
 
double BesselTransformer::GetDeltaResponse(double t) const 
{
    double result = 0.;
    double omz = fNormCutFreq*2*M_PI;
    for(size_t k = 0; k < fAN; k++) {
        double tmp = 0.;
        tmp += fAp[k] * cos(fAphi[k]*omz*t);
        tmp += fAq[k] * sin(fAphi[k]*omz*t);
        tmp *= exp(-fAtheta[k]*omz*t);
        result += tmp;
    }
    result *= 2*omz;
    return result;
}
 
 
 
double BesselTransformer::GetExpResponse(double t_0, double tau, double t) const 
{
    if(t <= t_0) return 0;
    double omz = fNormCutFreq*2*M_PI;
    double zs = t - t_0;
    double result = 0.;
    for(size_t k = 0; k < fAN; k++) {
        double c = fAphi[k] * omz;
        double &p = fAp[k];
        double &q = fAq[k];
        double alpha = -fAtheta[k] * omz;
        double alphatau = alpha - 1./tau;
        double sina = sin(c*zs)*(p*c + q*alphatau);
        double cosa = cos(c*zs)*(p*alphatau - q*c);
        double sx =  exp(zs*alphatau)*(sina + cosa);
        double dx = q*c - p*alphatau;
        result +=  (sx +  dx) / (c*c + alphatau*alphatau);
    }
    result *= 2*omz * exp((t-t_0)/tau);
 
    return result;
}
 
QError BesselTransformer::ConvolveTD(const QVector& pulse, int t_min, QVector& signal) 
{
 
    if(pulse.Size() != fSize) return QERR_SIZE_NOT_MATCH;
    signal.Resize(fSize);
    // MV FIXME
    size_t t_g_max = fSize - 1;; 
    size_t t_g_min = 0;
    for(unsigned int t = 0; t < fSize; t++) {
        int  t_max = fSize;
        int  t_inf,t_sup;
 
        if(t< t_min + t_g_min) {continue;}
        else if (t >= t_g_min + t_min && t< t_min + t_g_max) {
            t_inf = t_min;
            t_sup = t - t_g_min;
        } else if (t>= t_min + t_g_max && t <= t_g_min + t_max) {
            t_inf = t - t_g_max;
            t_sup = t - t_g_min;
 
        } else if (t > t_g_min + t_max && t <= t_max + t_g_max) {
            t_inf = t - t_g_max;
            t_sup = t_max;
        } else continue;
        double fc = 0.;
        for(int m = t_inf; m <= t_sup; m++) {
            fc += pulse[m]*fDeltaResponse[t-m];
        }   
        signal[t] = fc;
    }
    return QERR_SUCCESS;
}
 
double BesselTransformer::ConvolvedFuncTD(double *t) 
{
    size_t index = (size_t)roundl(double((*t)*fConvolvedVectorTD.Size())/fDuration);
    double ret = fConvolvedVectorTD[index];
    return ret;
}
 
double BesselTransformer::ConvolvedFuncFD(double *t) 
{
    int index = (size_t)roundl(double((*t)*fConvolvedVectorFD.Size())/fDuration);
//    if(index < 0) return fConvolvedVectorFD[0];
//    if(index >= (int)fSize) return fConvolvedVectorFD[fSize-1];
    
    return fConvolvedVectorFD[index];
}
 
void BesselTransformer::ConvolveTD(double (*pulse)(double*)) 
{
    QVector pulseVec(fSize);
    for(size_t i = 0; i < fSize; i++) {
        double convT = double(i)/fSamplingFrequency;
        pulseVec[i] = (*pulse)(&convT);
    }
    ConvolveTD(pulseVec,fConvolvedVectorTD);
}
 
void BesselTransformer::ConvolveFD(double (*pulse)(double*)) 
{
    QVector pulseVec(fSize);
    for(size_t i = 0; i < fSize; i++) {
        double convT = double(i)/fSamplingFrequency;
        pulseVec[i] = (*pulse)(&convT);
    }
    ConvolveFD(pulseVec,fConvolvedVectorFD);
}
 
#ifndef __CINT__
void BesselTransformer::ConvolveTD(pt2FuncParam func,double* param) 
{
    QVector pulseVec(fSize);
    for(size_t i = 0; i < fSize; i++) {
        double convT = double(i)/fSamplingFrequency;
        pulseVec[i] = (*func)(&convT,param);
    }
    ConvolveTD(pulseVec,fConvolvedVectorTD);
}
 
 
void BesselTransformer::ConvolveFD(pt2FuncParam func,double* param) 
{
    QVector pulseVec(fSize);
   for(size_t i = 0; i < fSize; i++) {
        double convT = double(i)/fSamplingFrequency;
        pulseVec[i] = (*func)(&convT,param);
    }
    ConvolveFD(pulseVec,fConvolvedVectorFD);
}
 
#endif
 
QError BesselTransformer::ConvolveFD(const QVector& pulse, QVector& signal) 
{
    if(pulse.Size() != fSize) return QERR_SIZE_NOT_MATCH;
 
    if(fZeroPad) {
        size_t newsize = 2;
        while(newsize < fSize) {
            newsize *=2;
        }
        newsize = fSize*2;
        QVector pulsez(newsize);
        QVectorView pulsez_first(pulsez, 0, newsize-fSize);
        QVectorView pulsez_second(pulsez, newsize-fSize, fSize);
        pulsez_first.GetVector().Initialize(pulse[0]);
        pulsez_second.GetVector() = pulse;
 
        QRealComplexFFT transformer(newsize);
        QVectorC pulseFD;
        transformer.TransformToFreq(pulsez, pulseFD);
        QVectorC convolved = pulseFD.Mult(fTransferFunction);
        //QVectorC convolved = pulseFD;
        QVector signalz(newsize);
        transformer.TransformFromFreq(convolved, signalz);
 
        QVectorView signalz_second(signalz, newsize-fSize, fSize);
        signal = signalz_second.GetVector();
    } else {
        QRealComplexFFT transformer(fSize);
        QVectorC pulseFD;
        transformer.TransformToFreq(pulse, pulseFD);
        QVectorC convolved = pulseFD.Mult(fTransferFunction);
        transformer.TransformFromFreq(convolved, signal);
    }
    return QERR_SUCCESS;
}