QNoisePowerSpectrumEstimator.cc

Go to the documentation of this file.

#include "QNoisePowerSpectrumEstimator.hh"
#include "QVectorC.hh"
 
using namespace Diana;
 
QNoisePowerSpectrumEstimator::QNoisePowerSpectrumEstimator(const int size,  Diana::QFFT::WindowType w, int gain) : fTransformer(size)
{
    fTransformer.SetWindowType(w, gain);
    fNoisePowerSpectrum.Resize(size);
    Clear();
}
void QNoisePowerSpectrumEstimator::Clear()
{
    const int size = fNoisePowerSpectrum.Size();
    fNoisePowerSpectrum.Clear();
    fNoisePowerSpectrum.Resize(size);
    fNoisePowerSpectrum.Initialize(0);
    fNumberOfPulses = 0;
}
 
QError QNoisePowerSpectrumEstimator::Add(const QVector& pulse)
{
    QError err = QERR_SUCCESS;
    if(pulse.Size() != fNoisePowerSpectrum.Size()) return QERR_SIZE_NOT_MATCH;
 
    double mean = pulse.Sum(pulse.Size())/pulse.Size();
    QVector meanPulse(pulse.Size()); meanPulse.Initialize(mean);
    QVector zeroMeanPulse = pulse;
    zeroMeanPulse -= meanPulse;
 
    QVectorC transformedPulse(pulse.Size());
    int ierr = fTransformer.TransformToFreq(zeroMeanPulse, transformedPulse);
    if(ierr != 0) return QERR_UNKNOWN_ERR;
    QVector realPart = transformedPulse.Re();
    QVector imagPart = transformedPulse.Im();
    for (size_t k = 0; k < fNoisePowerSpectrum.Size(); k++) { 
        double magSquared = realPart[k]*realPart[k] + imagPart[k]*imagPart[k];
        fNoisePowerSpectrum[k] += magSquared;
    }
    fNumberOfPulses++;
    return err;
}
 
QVector QNoisePowerSpectrumEstimator::GetNoisePowerSpectrum() const
{
    if(fNumberOfPulses <= 0) return fNoisePowerSpectrum;
    QVector ret = fNoisePowerSpectrum;
    ret /= fNumberOfPulses;
 
    return ret;
}