QFir.cc

Go to the documentation of this file.

#include "QFir.hh"
#include <sstream>
#include <cmath>
#include <gsl/gsl_sf_bessel.h>
#include <iostream>
 
QObjectImp(QFir);
 
 
using namespace Diana;
using namespace std;
 
QFir::QFir(const QCFirData& configuration) :
fConfiguration(configuration),
fCoefficients(configuration.fM)
{
   
    if(fConfiguration.fM%2 == 0) {
        QError err = QERR_SIZE_NOT_MATCH;
        std::stringstream msg;
        msg<<"Size of filter must be odd. You supplied ("<<fConfiguration.fM<<")";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        throw err;
    }
    fFilterReduction = 0;
    if(fConfiguration.fMethod == QCFirData::KayserBessel) KayserBesselAlgorithm();
    QVector theta(fConfiguration.fM*2), outputTheta;
    theta.Initialize(0);
    for(int i = fConfiguration.fM; i < 2*fConfiguration.fM; i++) {
        theta[i] = 1;
    }
    Filter(theta,outputTheta);
    double t10=-1,t50=-1,t90=-1;
    for(size_t i = 0; i < outputTheta.Size(); i++) {
        if(outputTheta[i]>0.9 && t90<0) {
            if(i>0) t90=i+(0.9-outputTheta[i])/(outputTheta[i]-outputTheta[i-1]);
            else t90 = i;
        }
        if(outputTheta[i]>0.5 && t50<0) {
            if(i>0) t50=i+(0.5-outputTheta[i])/(outputTheta[i]-outputTheta[i-1]);
            else t50 = i;
        }
        if(outputTheta[i]>0.1 && t10<0) {
            if(i>0) t10=i+(0.1-outputTheta[i])/(outputTheta[i]-outputTheta[i-1]);
            else t10 = i;
        }
    }
    fDelay = t50-fConfiguration.fM;
    fRiseTime = t90-t10;
}
 
void QFir::Clear()
{
    fConfiguration.Clear();
    fCoefficients.clear();
    fFilteredSize = 0;
    fFilterReduction = 0;
    fRiseTime = Q_DOUBLE_DEFAULT;
    fDelay = Q_DOUBLE_DEFAULT;
}
 
 
QError QFir::Filter(const QVector& input, QVector& output) const
{
    size_t M = fCoefficients.size();
    if(input.Size() <= fFilterReduction) {
        QError err = QERR_SIZE_NOT_MATCH;
        std::stringstream msg;
        msg<<"Input vector size ("<<input.Size()<<")"
            <<" is <= the filter reduction ("<<fFilterReduction<<")"; 
        err.SetDescription(__FILE__,__LINE__,msg.str());
        return err;
    }
    size_t N = input.Size()-fFilterReduction;
    output.Resize(N);
    for(size_t i = 0; i < N; i++) {
        output[i] = 0;
        for(size_t j = 0; j < M; j++) {
            output[i] += fCoefficients[j]*input[i+j];
        }
    }
 
    return QERR_SUCCESS;
}
 
void QFir::KayserBesselAlgorithm()
{
    // from http://www.arc.id.au/FilterDesign.html
    int M = fCoefficients.size();
    int Np = (M-1)/2;
 
    // filter
    double period = (2.*M_PI);
    fCoefficients[Np] = period*fConfiguration.fCutoff/M_PI ;
    for(int j=1; j<Np; j++)
        fCoefficients[j+Np] = sin(period*fConfiguration.fCutoff*(j))/(M_PI*(j));
 
    // window
    double Alpha;
    if (fConfiguration.fAttDB<21)
        Alpha = 0;
    else if (fConfiguration.fAttDB>50)
        Alpha = 0.1102*(fConfiguration.fAttDB-8.7);
    else
        Alpha = 0.5842*pow((fConfiguration.fAttDB-21), 0.4)+0.07886*(fConfiguration.fAttDB-21);
    double IAlpha = gsl_sf_bessel_I0 (Alpha);
 
    for(int j=0; j<=Np; j++){
        double arg = Alpha*sqrt(1.-pow(double(j)/Np,2));
        double window = gsl_sf_bessel_I0(arg)/IAlpha;
        //        cout<<"j: "<<j<<" w:"<<window<<" s:"<<fCoefficients[j+Np]<<endl;
        fCoefficients[j+Np] *= window;
    }
 
    // symmetrize
    for(int j = 0; j < Np; j++){
        fCoefficients[j] = fCoefficients[M-j-1];
    }
 
    // correct gain if != 1
    double area = 0;
    for(size_t j = 0; j < fCoefficients.size(); j++) {
        area += fCoefficients[j];
    }
    for(size_t j = 0; j < fCoefficients.size(); j++) {
        fCoefficients[j]/=area;
    }
 
 
    fFilterReduction = M-1;
    /*
       for(int j = 0; j < M; j++) {
       std::cout<<fCoefficients[j]<<" ";
           }
       */
    //    std::cout<<std::endl;
}