QMultiChannelDecorrelator.cc

Go to the documentation of this file.

#include "QMultiChannelDecorrelator.hh"
#include <sstream>
 
using namespace Diana;
 
QMultiChannelDecorrelator::QMultiChannelDecorrelator(const int chan, const QChannelCovariance* cov) 
{
    std::vector<int> channels;
    channels.push_back(chan);
    Constructor(channels,cov,false);
}
 
QMultiChannelDecorrelator::QMultiChannelDecorrelator(const int chan, const QChannelCovariance* cov, bool init)
{
    std::vector<int> channels;
    channels.push_back(chan);
    Constructor(channels,cov,init);
}
 
QMultiChannelDecorrelator::QMultiChannelDecorrelator(const std::vector<int>& channels, const QChannelCovariance* cov, bool init)
{
    Constructor(channels,cov,init);
}
 
 
void QMultiChannelDecorrelator::Constructor(const std::vector<int>& channels, const QChannelCovariance* cov, bool init)
{
    fScale = 1.;
    fTransformer = 0;
    fOriginalResolution = Q_DOUBLE_DEFAULT;
    fFilteredResolution = Q_DOUBLE_DEFAULT;
    if(cov == 0) {
        QError err = QERR_UNKNOWN_ERR;
        std::stringstream msg;
        msg<<"Pointer to QChannelCovariance is NULL";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
    }
    if(channels.empty()) {
        QError err = QERR_UNKNOWN_ERR;
        std::stringstream msg;
        msg<<"no channel given";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
    }
    fCovarianceMatrix = cov;
    fTransformer = new QRealComplexFFT(fCovarianceMatrix->GetNFreq());
    fAvailableSideChannels = fCovarianceMatrix->GetChannels();
    fChannels.clear();
    for(size_t i = 0; i < channels.size() ;i++) {
        const int chan = channels[i];
        fAvailableSideChannels.erase(FindChannel(fAvailableSideChannels,chan));
        fChannels.push_back(chan);
    }
    if(init) {
        std::vector<int> empty;
        SetSideChannels(empty);
    }
}
 
QMultiChannelDecorrelator::~QMultiChannelDecorrelator()
{
    if(fTransformer) delete fTransformer;
}
 
 
void QMultiChannelDecorrelator::SetSideChannels(const std::vector<int>& channels)
{
    fSideChannels.clear();
    fWeights.clear(); 
    fOriginalPowerSpectrum.Clear();
    fDecorrelatedPowerSpectrum.Clear();
 
    std::vector<int> allChannels;
    for(size_t i = 0; i < fChannels.size(); i++) {
        allChannels.push_back(fChannels.at(i));
    }
    for(size_t i = 0;  i < channels.size(); i++) {
        const int sideChan = channels[i];
        FindChannel(fAvailableSideChannels, sideChan);
        allChannels.push_back(sideChan); 
        fSideChannels.push_back(sideChan);
    }
    Init(allChannels);
}
 
void QMultiChannelDecorrelator::Init(const std::vector<int>& allChannels)
{
    const int chan = fChannels.at(0);
    fOriginalPowerSpectrum = fCovarianceMatrix->GetCovariance(chan,chan).GetModulus()*fScale*fScale;
    fOriginalResolution = ComputeResolution(fOriginalPowerSpectrum);
 
    QChannelCovariance reducedCovarianceMatrix = (fCovarianceMatrix->GetSubMatrix(allChannels)).Inverse();
    QVectorC Cii = reducedCovarianceMatrix.GetCovariance(chan,chan);
    QVectorC identity(fOriginalPowerSpectrum.Size());
    identity.Initialize(1,0);
    fDecorrelatedPowerSpectrum = (identity.Div(Cii)).GetModulus()*fScale*fScale;
    fFilteredResolution = ComputeResolution(fDecorrelatedPowerSpectrum);
 
    fWeights.push_back(identity);
    for(size_t i = 1;  i < allChannels.size(); i++) {
        QVectorC W = reducedCovarianceMatrix.GetCovariance(chan,allChannels[i]).Div(Cii);
        W[0].SetIm(0);
        W[0].SetRe(0);
        fWeights.push_back(W);
    }
 
}
 
double QMultiChannelDecorrelator::ComputeResolution(const QVector& ps) const
{
    double resolution = 0;
    for(size_t i = 1; i < ps.Size(); i++){
        resolution += ps[i];
    }
    resolution = sqrt(resolution);
    resolution /= ps.Size();
    return resolution;
}
 
void QMultiChannelDecorrelator::Filter(const Diana::QVector& input, const std::vector<const Diana::QVector*>& sideInputs, QVector& output) const
{
    if(fWeights.size() <= 1 || sideInputs.size() != fWeights.size()-1) {
        output.Clear();
        return;
    }
    QVectorC outFreq(fWeights[0].Size());
    QVectorC tmp(fWeights[0].Size());
 
    fTransformer->TransformToFreq(input,tmp);
    outFreq = tmp.Mult(fWeights[0]);
    for(size_t i = 1; i < fWeights.size(); i++) {
        fTransformer->TransformToFreq(*sideInputs[i-1],tmp);
        outFreq += tmp.Mult(fWeights[i]);
    }
    fTransformer->TransformFromFreq(outFreq,output);
}
 
void QMultiChannelDecorrelator::Filter(const std::vector<const Diana::QVector*>& inputs, const std::vector<const Diana::QVector*>& sideInputs, Diana::QVector& output) const
{
    if(inputs.size() != 1) {
        QError err = QERR_SIZE_NOT_MATCH;
        std::stringstream msg;
        msg<<"inputs.Size() = "<<inputs.size()<<" differs from 1";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
        return;    
    }
    Filter(*inputs[0],sideInputs,output);
}
 
std::vector<int>::iterator QMultiChannelDecorrelator::FindChannel(std::vector<int>& channels, const int chan)
{
    std::vector<int>::iterator iter = find(channels.begin(), channels.end(), chan);
    if(iter == channels.end()) {
        QError err = QERR_OUT_OF_RANGE;
        std::stringstream msg;
        msg<<"Channel "<<chan<<" is not contained in input covariance matrix ";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
    }
    return iter;
}