QMultiChannelDecorrelatorOF.cc

Go to the documentation of this file.

#include "QMultiChannelDecorrelatorOF.hh"
#include <sstream>
#include <algorithm>
using namespace Diana;
using namespace std;
 
QMultiChannelDecorrelatorOF::QMultiChannelDecorrelatorOF(const int chan, const QChannelCovariance* cov, const Diana::QVector& averagePulse) : QMultiChannelDecorrelator(chan,cov,false)
{
    fAveragePulses.push_back(averagePulse);
    std::vector<int> empty;
    SetSideChannels(empty);
}
 
QMultiChannelDecorrelatorOF::QMultiChannelDecorrelatorOF(const std::vector<int>& channels, const QChannelCovariance* cov, const std::vector<Diana::QVector>& averagePulses) : QMultiChannelDecorrelator(channels,cov,false)
{
    if(channels.size() != averagePulses.size()) {
        QError err = QERR_OUT_OF_RANGE;
        std::stringstream msg;
        msg<<"averagePulses.size() differs from channels.size(): "<<channels.size()<<" vs "<<averagePulses.size();
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
    }
    fAveragePulses = averagePulses;
    std::vector<int> empty;
    SetSideChannels(empty);
}
 
const Diana::QVector& QMultiChannelDecorrelatorOF::GetAveragePulse(const int chan) const
{
    std::vector<int>::const_iterator iter = std::find(fChannels.begin(),fChannels.end(),chan);
    if(iter == fChannels.end()) {
        QError err = QERR_OUT_OF_RANGE;
        std::stringstream msg;
        msg<<"Channel "<<chan<<" not present in fChannels";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
    }
    return fAveragePulses[*iter];
}
 
void QMultiChannelDecorrelatorOF::Init(const std::vector<int>& allChannels)
{
    
    fOriginalPowerSpectrum = fCovarianceMatrix->GetCovariance(allChannels[0],allChannels[0]).GetModulus()*fScale*fScale;
 
    fOriginalResolution = ComputeResolution(fOriginalPowerSpectrum);
    QChannelCovariance reducedCovarianceMatrix = (fCovarianceMatrix->GetSubMatrix(allChannels)).Inverse();
    size_t size = reducedCovarianceMatrix.GetNFreq();
 
    vector<QVectorC> APFreq(fAveragePulses.size());
    double max;
    for(size_t i = 0; i < fAveragePulses.size(); i++) {
        QVectorC SFreq;
        QVector normAP = fAveragePulses[i];
        int maxpos = normAP.GetMaxIndex();
        QVector baseline(size);
        size_t baseLength = 3*maxpos/4;
        baseline.Initialize(normAP.Sum(baseLength,0)/baseLength);
        normAP -= baseline;
        if(i == 0) {
            fShift = maxpos;
            max = normAP.GetMax();
        }
        normAP /= normAP.GetMax();
        fTransformer->TransformToFreq(normAP,SFreq);
        SFreq[0].SetRe(0);
        SFreq[0].SetIm(0);
        APFreq[i] = SFreq;
    }
 
        
    for(size_t j = 0; j < allChannels.size(); j++) {
        int chanJ = allChannels[j];
        QVectorC W(size);
        W.Initialize(0,0);
        for(size_t i = 0; i < fAveragePulses.size(); i++) {
            int chanI = fChannels[i];
            W+=APFreq[i].Conj().Mult(reducedCovarianceMatrix.GetCovariance(chanJ,chanI));
        }
        W[0].SetRe(0);
        W[0].SetIm(0);
        fWeights.push_back(W);
    }
    
    fNorm = 1.;
    vector<const QVector*> inputs,sideInputs;
    for(size_t i = 0; i < fAveragePulses.size(); i++) {
        inputs.push_back(&fAveragePulses[i]);
    }
    QVector output;
    Filter(inputs,sideInputs,output);
    fNorm = max*fScale/output.GetMax();
 
    QVectorC deco(size);
    deco.Initialize(0,0);
    for(size_t i = 0; i < fAveragePulses.size(); i++) {
        deco += fWeights[i].Mult(APFreq[i]);
    }
    fDecorrelatedPowerSpectrum = deco.GetModulus()*fNorm*fNorm*fScale*fScale;
    fFilteredResolution = ComputeResolution(fDecorrelatedPowerSpectrum);
}
 
void QMultiChannelDecorrelatorOF::Filter(const std::vector<const Diana::QVector*>& inputs, const std::vector<const Diana::QVector*>& sideInputs, Diana::QVector& output) const
{
    QVectorC tmp, outputFreq(fCovarianceMatrix->GetNFreq());
    outputFreq.Initialize(0,0);
    size_t i,j;
    for(i =0; i < inputs.size(); i++) {
        fTransformer->TransformToFreq(*inputs[i],tmp);
        outputFreq += fWeights[i].Mult(tmp);
    }
    for(j =0; j < sideInputs.size(); j++) {
        fTransformer->TransformToFreq(*sideInputs[j],tmp);
        outputFreq += fWeights[i+j].Mult(tmp);
    }
    fTransformer->TransformFromFreq(outputFreq,output);
    output *= fNorm*fScale;
    output.Shift(fShift);
}