From 8704b736d63ccc5c5d4452df5c705d48571b8e6c Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Thu, 15 Apr 2021 21:17:37 +0200 Subject: [PATCH 01/48] Modify KTFrequencySpectrumFFTW to use Eigen --- .../Data/Transform/KTFrequencySpectrumFFTW.cc | 101 ++++++++---------- .../Data/Transform/KTFrequencySpectrumFFTW.hh | 67 ++++++------ 2 files changed, 84 insertions(+), 84 deletions(-) diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc index 7aa45b4f7..343acbb66 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc @@ -24,7 +24,7 @@ namespace Katydid KTLOGGER(fslog, "KTFrequencySpectrumFFTW"); KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW() : - KTPhysicalArray< 1, fftw_complex >(), + KTPhysicalArray< 1, std::complex >(), KTFrequencySpectrum(), fIsArrayOrderFlipped(false), fIsSizeEven(true), @@ -38,7 +38,7 @@ namespace Katydid } KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(size_t nBins, double rangeMin, double rangeMax, bool arrayOrderIsFlipped) : - KTPhysicalArray< 1, fftw_complex >(nBins, rangeMin, rangeMax), + KTPhysicalArray< 1, std::complex >(nBins, rangeMin, rangeMax), KTFrequencySpectrum(), fIsArrayOrderFlipped(arrayOrderIsFlipped), fIsSizeEven(nBins%2 == 0), @@ -64,14 +64,11 @@ namespace Katydid KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(std::initializer_list value, size_t nBins, double rangeMin, double rangeMax, bool arrayOrderIsFlipped) : KTFrequencySpectrumFFTW(nBins, rangeMin, rangeMax, arrayOrderIsFlipped) { - for (unsigned index = 0; index < nBins; ++index) - { - std::copy(value.begin(), value.end(), fData[index]); - } + std::copy(value.begin(), value.end(), this->begin()); } KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig) : - KTPhysicalArray< 1, fftw_complex >(orig), + KTPhysicalArray< 1, std::complex >(orig), KTFrequencySpectrum(), fIsArrayOrderFlipped(orig.fIsArrayOrderFlipped), fIsSizeEven(orig.fIsSizeEven), @@ -88,18 +85,20 @@ namespace Katydid { } - KTFrequencySpectrumFFTW& KTFrequencySpectrumFFTW::operator=(const KTFrequencySpectrumFFTW& rhs) - { - KTPhysicalArray< 1, fftw_complex >::operator=(rhs); - fIsArrayOrderFlipped = rhs.fIsArrayOrderFlipped; - fIsSizeEven = rhs.fIsSizeEven; - fLeftOfCenterOffset = rhs.fLeftOfCenterOffset; - fCenterBin = rhs.fCenterBin; - fConstBinAccess = rhs.fConstBinAccess; - fBinAccess = rhs.fBinAccess; - fNTimeBins = rhs.fNTimeBins; - return *this; - } + //Is the copy assignment operator necessary? I think the default one will do the right thing + + //~ KTFrequencySpectrumFFTW& KTFrequencySpectrumFFTW::operator=(const KTFrequencySpectrumFFTW& rhs) + //~ { + //~ KTPhysicalArray< 1, std::complex >::operator=(rhs); + //~ fIsArrayOrderFlipped = rhs.fIsArrayOrderFlipped; + //~ fIsSizeEven = rhs.fIsSizeEven; + //~ fLeftOfCenterOffset = rhs.fLeftOfCenterOffset; + //~ fCenterBin = rhs.fCenterBin; + //~ fConstBinAccess = rhs.fConstBinAccess; + //~ fBinAccess = rhs.fBinAccess; + //~ fNTimeBins = rhs.fNTimeBins; + //~ return *this; + //~ } const KTAxisProperties< 1 >& KTFrequencySpectrumFFTW::GetAxis() const { @@ -113,14 +112,9 @@ namespace Katydid KTFrequencySpectrumFFTW& KTFrequencySpectrumFFTW::CConjugate() { - unsigned nBins = size(); -#pragma omp parallel for - for (unsigned iBin=0; iBin {0.0}; + return *this; } KTFrequencySpectrumFFTW& KTFrequencySpectrumFFTW::Scale(double scale) { - unsigned nBins = size(); -#pragma omp parallel for - for (unsigned iBin=0; iBin Replace KTFrequencySpectrumPolar KTFrequencySpectrumPolar* KTFrequencySpectrumFFTW::CreateFrequencySpectrumPolar() const { unsigned nBins = size(); @@ -169,7 +149,7 @@ namespace Katydid #pragma omp parallel for for (unsigned iBin=0; iBinGetData().setZero(); //probably unnecessary, Eigen array should be zero-initialized int dcBin = FindBin(0.); // default case: dcBin >= 0 && dcBin < size() @@ -214,22 +197,32 @@ namespace Katydid double scaling = 1. / KTPowerSpectrum::GetResistance() / (double)GetNTimeBins(); + // to replace after Eigen replaces all KTPhysicalArrays + //int nPosBins = lastPosFreqBin - firstPosFreqBin; + //int nNegBins = lastNegFreqBin - firstNegFreqBin; + // + //newPS->GetData().segment(firstPosFreqBin, nPosBins) = + // fData.segment(firstPosFreqBin, nPosBins).abs2() * scaling; + //newPS->GetData().segment(firstNegFreqBin, nNegBins) = + // fData.segment(firstNegFreqBin, nNegBins).abs2() * scaling; + + // to replace after Eigen replaces all KTPhysicalArrays double valueImag, valueReal; #pragma omp parallel for private(valueReal, valueImag) for (unsigned iBin = firstPosFreqBin; iBin < lastPosFreqBin; ++iBin) { - valueReal = (*this)(iBin)[0]; - valueImag = (*this)(iBin)[1]; + valueReal = (*this)(iBin).real(); + valueImag = (*this)(iBin).imag(); (*newPS)(iBin) = (valueReal * valueReal + valueImag * valueImag) * scaling; } #pragma omp parallel for private(valueReal, valueImag) for (unsigned iBin = firstNegFreqBin; iBin < lastNegFreqBin; ++iBin) { - valueReal = (*this)(iBin)[0]; - valueImag = (*this)(iBin)[1]; + valueReal = (*this)(iBin).real(); + valueImag = (*this)(iBin).imag(); (*newPS)(iBin) = (valueReal * valueReal + valueImag * valueImag) * scaling; } - + return newPS; } diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh index 42893f285..4092683f9 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh @@ -10,7 +10,7 @@ #define KTFREQUENCYSPECTRUMFFTW_HH_ #include "KTFrequencySpectrum.hh" -#include "KTPhysicalArrayFFTW.hh" +#include "KTPhysicalArrayComplex.hh" #include #include @@ -20,7 +20,7 @@ namespace Katydid class KTFrequencySpectrumPolar; - class KTFrequencySpectrumFFTW : public KTPhysicalArray< 1, fftw_complex >, public KTFrequencySpectrum + class KTFrequencySpectrumFFTW : public KTPhysicalArray< 1, std::complex >, public KTFrequencySpectrum { public: KTFrequencySpectrumFFTW(); @@ -50,8 +50,8 @@ namespace Katydid // replace some of the KTPhysicalArray interface - const fftw_complex& operator()(unsigned i) const; - fftw_complex& operator()(unsigned i); + const std::complex& operator()(unsigned i) const; + std::complex& operator()(unsigned i); virtual double GetReal(unsigned bin) const; virtual double GetImag(unsigned bin) const; @@ -70,22 +70,23 @@ namespace Katydid virtual void SetNTimeBins(unsigned bins); private: - typedef const fftw_complex& (KTFrequencySpectrumFFTW::*ConstOrderedBinAccessFunc)(unsigned) const; - typedef fftw_complex& (KTFrequencySpectrumFFTW::*OrderedBinAccessFunc)(unsigned); + typedef const std::complex& (KTFrequencySpectrumFFTW::*ConstOrderedBinAccessFunc)(unsigned) const; + typedef std::complex& (KTFrequencySpectrumFFTW::*OrderedBinAccessFunc)(unsigned); ConstOrderedBinAccessFunc fConstBinAccess; OrderedBinAccessFunc fBinAccess; - const fftw_complex& ReorderedBinAccess(unsigned i) const; - fftw_complex& ReorderedBinAccess(unsigned i); + const std::complex& ReorderedBinAccess(unsigned i) const; + std::complex& ReorderedBinAccess(unsigned i); - const fftw_complex& AsIsBinAccess(unsigned i) const; - fftw_complex& AsIsBinAccess(unsigned i); + const std::complex& AsIsBinAccess(unsigned i) const; + std::complex& AsIsBinAccess(unsigned i); public: // normal KTFrequencySpectrumPolar functions - virtual KTFrequencySpectrumFFTW& operator=(const KTFrequencySpectrumFFTW& rhs); + //see comment in cc file + //virtual KTFrequencySpectrumFFTW& operator=(const KTFrequencySpectrumFFTW& rhs); /// In-place calculation of the complex conjugate virtual KTFrequencySpectrumFFTW& CConjugate(); @@ -103,7 +104,7 @@ namespace Katydid unsigned fNTimeBins; protected: - mutable const fftw_complex* fPointCache; + mutable const std::complex* fPointCache; }; inline bool KTFrequencySpectrumFFTW::GetIsArrayOrderFlipped() const @@ -131,71 +132,77 @@ namespace Katydid return GetDataLabel(); } - inline const fftw_complex& KTFrequencySpectrumFFTW::operator()(unsigned i) const + inline const std::complex& KTFrequencySpectrumFFTW::operator()(unsigned i) const { return (this->*fConstBinAccess)(i); } - inline fftw_complex& KTFrequencySpectrumFFTW::operator()(unsigned i) + inline std::complex& KTFrequencySpectrumFFTW::operator()(unsigned i) { return (this->*fBinAccess)(i); } - inline const fftw_complex& KTFrequencySpectrumFFTW::ReorderedBinAccess(unsigned i) const + inline const std::complex& KTFrequencySpectrumFFTW::ReorderedBinAccess(unsigned i) const { return (i >= fCenterBin) ? fData[i - fCenterBin] : fData[i + fLeftOfCenterOffset]; } - inline fftw_complex& KTFrequencySpectrumFFTW::ReorderedBinAccess(unsigned i) + inline std::complex& KTFrequencySpectrumFFTW::ReorderedBinAccess(unsigned i) { return (i >= fCenterBin) ? fData[i - fCenterBin] : fData[i + fLeftOfCenterOffset]; } - inline const fftw_complex& KTFrequencySpectrumFFTW::AsIsBinAccess(unsigned i) const + inline const std::complex& KTFrequencySpectrumFFTW::AsIsBinAccess(unsigned i) const { return fData[i]; } - inline fftw_complex& KTFrequencySpectrumFFTW::AsIsBinAccess(unsigned i) + inline std::complex& KTFrequencySpectrumFFTW::AsIsBinAccess(unsigned i) { return fData[i]; } inline double KTFrequencySpectrumFFTW::GetReal(unsigned bin) const { - return (*this)(bin)[0]; + return (*this)(bin).real(); } inline double KTFrequencySpectrumFFTW::GetImag(unsigned bin) const { - return (*this)(bin)[1]; + return (*this)(bin).imag(); } inline void KTFrequencySpectrumFFTW::SetRect(unsigned bin, double real, double imag) { - fPointCache = &(*this)(bin); - (*const_cast< fftw_complex* >(fPointCache))[0] = real; - (*const_cast< fftw_complex* >(fPointCache))[1] = imag; + //~ fPointCache = &(*this)(bin); + //~ (*const_cast< fftw_complex* >(fPointCache))[0] = real; + //~ (*const_cast< fftw_complex* >(fPointCache))[1] = imag; + + (*this)(bin) = std::complex(real, imag); return; } inline double KTFrequencySpectrumFFTW::GetAbs(unsigned bin) const { - fPointCache = &(*this)(bin); - return sqrt((*fPointCache)[0]*(*fPointCache)[0] + (*fPointCache)[1]*(*fPointCache)[1]); + //fPointCache = &(*this)(bin); + //return sqrt((*fPointCache)[0]*(*fPointCache)[0] + (*fPointCache)[1]*(*fPointCache)[1]); + return std::abs((*this)(bin)); } inline double KTFrequencySpectrumFFTW::GetArg(unsigned bin) const { - fPointCache = &(*this)(bin); - return atan2((*fPointCache)[1], (*fPointCache)[0]); + //fPointCache = &(*this)(bin); + //return atan2((*fPointCache)[1], (*fPointCache)[0]); + return std::arg((*this)(bin)); } inline void KTFrequencySpectrumFFTW::SetPolar(unsigned bin, double abs, double arg) { - fPointCache = &(*this)(bin); - (*const_cast< fftw_complex* >(fPointCache))[0] = abs * cos(arg); - (*const_cast< fftw_complex* >(fPointCache))[1] = abs * sin(arg); + //~ fPointCache = &(*this)(bin); + //~ (*const_cast< fftw_complex* >(fPointCache))[0] = abs * cos(arg); + //~ (*const_cast< fftw_complex* >(fPointCache))[1] = abs * sin(arg); + + (*this)(bin) = std::polar(abs, arg); return; } From 59e148709f66edf850b9157e70016f7ac415c27a Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Fri, 16 Apr 2021 15:44:52 +0200 Subject: [PATCH 02/48] Modify KTPhysicalArryComplex Moves all function definitions to the cc file --- Source/Utility/KTPhysicalArrayComplex.cc | 707 +++++++++++++++++++++++ Source/Utility/KTPhysicalArrayComplex.hh | 706 +--------------------- 2 files changed, 722 insertions(+), 691 deletions(-) diff --git a/Source/Utility/KTPhysicalArrayComplex.cc b/Source/Utility/KTPhysicalArrayComplex.cc index 6e55ff166..cc4f74ad4 100644 --- a/Source/Utility/KTPhysicalArrayComplex.cc +++ b/Source/Utility/KTPhysicalArrayComplex.cc @@ -42,6 +42,240 @@ namespace Katydid KTPhysicalArray< 1, value_type >::~KTPhysicalArray() {} + const Eigen::Array< std::complex, Eigen::Dynamic, 1, Eigen::ColMajor >& KTPhysicalArray< 1, std::complex >::GetData() const + { + return fData; + } + + Eigen::Array< std::complex, Eigen::Dynamic, 1, Eigen::ColMajor >& KTPhysicalArray< 1, std::complex >::GetData() + { + return fData; + } + + const std::string& KTPhysicalArray< 1, std::complex >::GetDataLabel() const + { + return fLabel; + } + + void KTPhysicalArray< 1, std::complex >::SetDataLabel(const std::string& label) + { + fLabel = label; + return; + } + + const std::complex& KTPhysicalArray< 1, std::complex >::operator()(unsigned i) const + { +#ifndef NDEBUG + if (i >= size()) + { + std::stringstream msg; + msg << "Out of bounds: " << i << " >= " << size(); + KTERROR(utillog_physarr, msg.str()); + throw std::out_of_range(msg.str()); + } +#endif + return fData[i]; + } + + std::complex& KTPhysicalArray< 1, std::complex >::operator()(unsigned i) + { +#ifndef NDEBUG + if (i >= size()) + { + std::stringstream msg; + msg << "Out of bounds: " << i << " >= " << size(); + KTERROR(utillog_physarr, msg.str()); + throw std::out_of_range(msg.str()); + } +#endif + return fData[i]; + } + + bool KTPhysicalArray< 1, std::complex >::IsCompatibleWith(const KTPhysicalArray< 1, std::complex >& rhs) const + { + //return (this->size() == rhs.size() && this->GetRangeMin() == rhs.GetRangeMin() && this->GetRangeMax() == GetRangeMax()); + return (this->size() == rhs.size()); + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator+=(const KTPhysicalArray< 1, std::complex>& rhs) + { + if (this->IsCompatibleWith(rhs)) + { + fData += rhs.fData; + } + + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator-=(const KTPhysicalArray< 1, std::complex>& rhs) + { + if (this->IsCompatibleWith(rhs)) + { + fData -= rhs.fData; + } + + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator*=(const KTPhysicalArray< 1, std::complex>& rhs) + { + if (this->IsCompatibleWith(rhs)) + { + fData *= rhs.fData; + } + + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator/=(const KTPhysicalArray< 1, std::complex>& rhs) + { + if (this->IsCompatibleWith(rhs)) + { + fData /= rhs.fData; + } + + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator+=(const std::complex& rhs) + { + fData += rhs; + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator-=(const std::complex& rhs) + { + fData -= rhs; + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator*=(const std::complex& rhs) + { + fData *= rhs; + return *this; + } + + KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator/=(const std::complex& rhs) + { + fData /= rhs; + return *this; + } + + //************************* + // Operator implementations + //************************* + + /// Add two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + KTPhysicalArray< 1, std::complex > operator+(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); + + lhs += rhs; + return lhs; + } + + /// Subtracts two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + KTPhysicalArray< 1, std::complex > operator-(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); + + lhs -= rhs; + return lhs; + } + + /// Multiplies two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + KTPhysicalArray< 1, std::complex > operator*(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); + + lhs *= rhs; + return lhs; + } + + /// Divides two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + KTPhysicalArray< 1, std::complex > operator/(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); + + lhs /= rhs; + return lhs; + } + + std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 1, std::complex >& rhs) + { + ostr << rhs.GetData(); + return ostr; + } + + //************************* + // Iterators + //************************* + + KTPhysicalArray< 1, std::complex >::const_iterator KTPhysicalArray< 1, std::complex >::begin() const + { + return fData.data(); + } + + KTPhysicalArray< 1, std::complex >::const_iterator KTPhysicalArray< 1, std::complex >::end() const + { + return fData.data() + size(); + } + + KTPhysicalArray< 1, std::complex >::iterator KTPhysicalArray< 1, std::complex >::begin() + { + return fData.data(); + } + + KTPhysicalArray< 1, std::complex >::iterator KTPhysicalArray< 1, std::complex >::end() + { + return fData.data() + size(); + } + + KTPhysicalArray< 1, std::complex >::const_reverse_iterator KTPhysicalArray< 1, std::complex >::rbegin() const + { + //a reverse iterator points to the first element outside the container + return std::reverse_iterator* >(fData.data() + size()); + } + + KTPhysicalArray< 1, std::complex >::const_reverse_iterator KTPhysicalArray< 1, std::complex >::rend() const + { + return std::reverse_iterator* >(fData.data()); + } + + KTPhysicalArray< 1, std::complex >::reverse_iterator KTPhysicalArray< 1, std::complex >::rbegin() + { + return std::reverse_iterator< std::complex* >(fData.data() + size()); + } + + KTPhysicalArray< 1, std::complex >::reverse_iterator KTPhysicalArray< 1, std::complex >::rend() + { + return std::reverse_iterator< std::complex* >(fData.data()); + } + + //************************* + // Begin, end free functions for range based for-loop + //************************* + + KTPhysicalArray< 1, std::complex >::iterator begin(KTPhysicalArray< 1, std::complex > & array) + { + return array.begin(); + } + + KTPhysicalArray< 1, std::complex >::iterator end(KTPhysicalArray< 1, std::complex > & array) + { + return array.end(); + } + + KTPhysicalArray< 1, std::complex >::const_iterator begin(const KTPhysicalArray< 1, std::complex > & array) + { + return array.begin(); + } + + KTPhysicalArray< 1, std::complex >::const_iterator end(const KTPhysicalArray< 1, std::complex > & array) + { + return array.end(); + } + //******************************* // 2D implementation @@ -87,5 +321,478 @@ namespace Katydid KTPhysicalArray< 2, std::complex >::~KTPhysicalArray() { } + + size_t KTPhysicalArray< 2, std::complex >::cols() const + { + return static_cast(fData.cols()); + } + + size_t KTPhysicalArray< 2, std::complex >::rows() const + { + return static_cast(fData.rows()); + } + + + const Eigen::Array< std::complex, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor >& KTPhysicalArray< 2, std::complex >::GetData() const + { + return fData; + } + + + typename Eigen::Array< std::complex, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor >& KTPhysicalArray< 2, std::complex >::GetData() + { + return fData; + } + + + const std::string& KTPhysicalArray< 2, std::complex >::GetDataLabel() const + { + return fLabel; + } + + + void KTPhysicalArray< 2, std::complex >::SetDataLabel(const std::string& label) + { + fLabel = label; + return; + } + + + const std::complex& KTPhysicalArray< 2, std::complex >::operator()(unsigned i, unsigned j) const + { +#ifdef Katydid_DEBUG + if (i >= size(1)) + { + std::stringstream msg; + msg << "Out of bounds on axis 1: " << i << " >= " << size(1); + KTERROR(utillog_physarr, msg.str()); + throw std::out_of_range(msg.str()); + } + if (j >= size(2)) + { + std::stringstream msg; + msg << "Out of bounds on axis 2: " << j << " >= " << size(2); + KTERROR(utillog_physarr, msg.str()); + throw std::out_of_range(msg.str()); + } +#endif + return fData(i,j); + } + + + std::complex& KTPhysicalArray< 2, std::complex >::operator()(unsigned i, unsigned j) + { +#ifdef Katydid_DEBUG + if (i >= size(1)) + { + std::stringstream msg; + msg << "Out of bounds on axis 1: " << i << " >= " << size(1); + KTERROR(utillog_physarr, msg.str()); + throw std::out_of_range(msg.str()); + } + if (j >= size(2)) + { + std::stringstream msg; + msg << "Out of bounds on axis 2: " << j << " >= " << size(2); + KTERROR(utillog_physarr, msg.str()); + throw std::out_of_range(msg.str()); + } +#endif + return fData(i,j); + } + + + bool KTPhysicalArray< 2, std::complex >::IsCompatibleWith(const KTPhysicalArray< 2, std::complex >& rhs) const + { + //return (this->size(1) == rhs.size(1) && this->GetRangeMin(1) == rhs.GetRangeMin(1) && this->GetRangeMax(1) == GetRangeMax(1) && + // this->size(2) == rhs.size(2) && this->GetRangeMin(2) == rhs.GetRangeMin(2) && this->GetRangeMax(2) == GetRangeMax(2)); + return (this->size(1) == rhs.size(1) && + this->size(2) == rhs.size(2)); + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator+=(const KTPhysicalArray< 2, std::complex>& rhs) + { + if ( this->IsCompatibleWith(rhs) ) + { + fData += rhs.fData; + } + + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator-=(const KTPhysicalArray< 2, std::complex>& rhs) + { + if ( this->IsCompatibleWith(rhs) ) + { + fData -= rhs.fData; + } + + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator*=(const KTPhysicalArray< 2, std::complex>& rhs) + { + if ( this->IsCompatibleWith(rhs) ) + { + fData *= rhs.fData; + } + + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator/=(const KTPhysicalArray< 2, std::complex>& rhs) + { + if ( this->IsCompatibleWith(rhs) ) + { + fData /= rhs.fData; + } + + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator+=(const std::complex& rhs) + { + fData += rhs; + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator-=(const std::complex& rhs) + { + fData -= rhs; + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator*=(const std::complex& rhs) + { + fData *= rhs; + return *this; + } + + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator/=(const std::complex& rhs) + { + fData /= rhs; + return *this; + } + + + KTPhysicalArray< 2, std::complex >::const_row_iterator KTPhysicalArray< 2, std::complex >::begin2() const + { + return const_row_iterator{ fData.data(), 1 }; + } + + + KTPhysicalArray< 2, std::complex >::const_col_iterator KTPhysicalArray< 2, std::complex >::begin1() const + { + return const_col_iterator{ fData.data(), static_cast(size(1)) }; + } + + + KTPhysicalArray< 2, std::complex >::const_row_iterator KTPhysicalArray< 2, std::complex >::end2() const + { + return const_row_iterator{ fData.data() + size(1), 1 }; + } + + + KTPhysicalArray< 2, std::complex >::const_col_iterator KTPhysicalArray< 2, std::complex >::end1() const + { + return const_col_iterator{ fData.data() + size(1)*size(2), static_cast(size(1)) }; + } + + + KTPhysicalArray< 2, std::complex >::row_iterator KTPhysicalArray< 2, std::complex >::begin2() + { + return row_iterator{ fData.data(), 1 }; + } + + + KTPhysicalArray< 2, std::complex >::col_iterator KTPhysicalArray< 2, std::complex >::begin1() + { + return col_iterator{ fData.data(), static_cast(size(1)) }; + } + + + KTPhysicalArray< 2, std::complex >::row_iterator KTPhysicalArray< 2, std::complex >::end2() + { + return row_iterator{ fData.data() + size(1), 1 }; + } + + + KTPhysicalArray< 2, std::complex >::col_iterator KTPhysicalArray< 2, std::complex >::end1() + { + return col_iterator{ fData.data() + size(1)*size(2), static_cast(size(1)) }; + } + + //************************* + // Reverse iterators + //************************* + + KTPhysicalArray< 2, std::complex >::const_reverse_row_iterator KTPhysicalArray< 2, std::complex >::rbegin2() const + { + return const_reverse_row_iterator{ fData.data() + size(1), -1 }; + } + + + KTPhysicalArray< 2, std::complex >::const_reverse_col_iterator KTPhysicalArray< 2, std::complex >::rbegin1() const + { + return const_reverse_col_iterator{ fData.data() + size(1)*size(2), -static_cast(size(1)) }; + } + + + KTPhysicalArray< 2, std::complex >::const_reverse_row_iterator KTPhysicalArray< 2, std::complex >::rend2() const + { + return const_reverse_row_iterator{ fData.data(), -1 }; + } + + + KTPhysicalArray< 2, std::complex >::const_reverse_col_iterator KTPhysicalArray< 2, std::complex >::rend1() const + { + return const_reverse_col_iterator{ fData.data(), -static_cast(size(1)) }; + } + + + KTPhysicalArray< 2, std::complex >::reverse_row_iterator KTPhysicalArray< 2, std::complex >::rbegin2() + { + return reverse_row_iterator{ fData.data() + size(1), -1 }; + } + + + KTPhysicalArray< 2, std::complex >::reverse_col_iterator KTPhysicalArray< 2, std::complex >::rbegin1() + { + return reverse_col_iterator{ fData.data() + size(1)*size(2), -static_cast(size(1)) }; + } + + + KTPhysicalArray< 2, std::complex >::reverse_row_iterator KTPhysicalArray< 2, std::complex >::rend2() + { + return reverse_row_iterator{ fData.data(), -1 }; + } + + + KTPhysicalArray< 2, std::complex >::reverse_col_iterator KTPhysicalArray< 2, std::complex >::rend1() + { + return reverse_col_iterator{ fData.data(), -static_cast(size(1)) }; + } + + + std::complex KTPhysicalArray< 2, std::complex >::GetMaximumBin(unsigned& maxXBin, unsigned& maxYBin) const + { + //return fData.maxCoeff(&maxXBin, &maxYBin); // can be used for real values + + //functor for finding the maximum + struct MaxVisitor + { + + unsigned row; + unsigned col; + value_type max; + + void init( const value_type& value, Eigen::Index i, Eigen::Index j ) + { + max=value; + row=i; + col=j; + }; + + void operator()( const value_type& value, Eigen::Index i, Eigen::Index j ) + { + if(std::abs(value)>std::abs(max)) + { + max=value; + row=i; + col=j; + } + }; + }; + + MaxVisitor mv; + fData.visit(mv); + maxXBin = mv.row; + maxYBin = mv.col; + + return mv.max; + } + + + std::complex KTPhysicalArray< 2, std::complex >::GetMinimumBin(unsigned& minXBin, unsigned& minYBin) const + { + //return fData.minCoeff(&minXBin, &minYBin); // can be used for real values + + //functor for finding the minimum + struct MinVisitor + { + + unsigned row; + unsigned col; + value_type min; + + void init( const value_type& value, Eigen::Index i, Eigen::Index j ) + { + min=value; + row=i; + col=j; + }; + + void operator()( const value_type& value, Eigen::Index i, Eigen::Index j ) + { + if(std::abs(value), std::complex > KTPhysicalArray< 2, std::complex >::GetMinMaxBin(unsigned& minXBin, unsigned& minYBin, unsigned& maxXBin, unsigned& maxYBin) + { + //functor for finding minimum and maximum at the same time + struct MinMaxVisitor + { + + unsigned rowMax; + unsigned colMax; + unsigned rowMin; + unsigned colMin; + value_type max; + value_type min; + + void init( const value_type& value, Eigen::Index i, Eigen::Index j ) + { + max=value; + min=value; + rowMax=i; + colMax=j; + rowMin=i; + colMin=j; + }; + + void operator()( const value_type& value, Eigen::Index i, Eigen::Index j ) + { + if(std::abs(value)>std::abs(max)) + { + max=value; + rowMax=i; + colMax=j; + } + if(std::abs(value) > operator+(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); + + lhs += rhs; + return lhs; + } + + /// Subtracts two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + + KTPhysicalArray< 2, std::complex > operator-(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); + + lhs -= rhs; + return lhs; + } + + /// Multiplies two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + + KTPhysicalArray< 2, std::complex > operator*(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); + + lhs *= rhs; + return lhs; + } + + /// Divides two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + + KTPhysicalArray< 2, std::complex > operator/(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); + + lhs /= rhs; + return lhs; + } + + std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 2, std::complex >& rhs) + { + ostr << rhs.GetData(); + return ostr; + } + + + //************************* + // Matrix and Vector operations + //************************* + + KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator%=(const KTPhysicalArray< 2, std::complex >& rhs) + { + + fData = fData.matrix()*rhs.fData.matrix(); + + SetRangeMin(2, rhs.GetRangeMin(2)); + SetRangeMax(2, rhs.GetRangeMax(2)); + + return *this; + } + + KTPhysicalArray< 2, std::complex > operator%(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) + { + + lhs %= rhs; + return lhs; + } + + //matrix-vector-multiplication + KTPhysicalArray< 1, std::complex > operator%(const KTPhysicalArray< 2, std::complex >& lhs, KTPhysicalArray< 1, std::complex > rhs) + { + + rhs.fData = lhs.fData.matrix()*rhs.fData.matrix(); + + rhs.SetRangeMin(lhs.GetRangeMin(1)); + rhs.SetRangeMax(lhs.GetRangeMax(1)); + + return rhs; + } } /* namespace Katydid */ diff --git a/Source/Utility/KTPhysicalArrayComplex.hh b/Source/Utility/KTPhysicalArrayComplex.hh index ff10fc076..dca0a0e07 100644 --- a/Source/Utility/KTPhysicalArrayComplex.hh +++ b/Source/Utility/KTPhysicalArrayComplex.hh @@ -137,241 +137,16 @@ namespace Katydid }; - - inline const Eigen::Array< std::complex, Eigen::Dynamic, 1, Eigen::ColMajor >& KTPhysicalArray< 1, std::complex >::GetData() const - { - return fData; - } - - inline Eigen::Array< std::complex, Eigen::Dynamic, 1, Eigen::ColMajor >& KTPhysicalArray< 1, std::complex >::GetData() - { - return fData; - } - - inline const std::string& KTPhysicalArray< 1, std::complex >::GetDataLabel() const - { - return fLabel; - } - - inline void KTPhysicalArray< 1, std::complex >::SetDataLabel(const std::string& label) - { - fLabel = label; - return; - } - - inline const std::complex& KTPhysicalArray< 1, std::complex >::operator()(unsigned i) const - { -#ifndef NDEBUG - if (i >= size()) - { - std::stringstream msg; - msg << "Out of bounds: " << i << " >= " << size(); - KTERROR(utillog_physarr, msg.str()); - throw std::out_of_range(msg.str()); - } -#endif - return fData[i]; - } - - inline std::complex& KTPhysicalArray< 1, std::complex >::operator()(unsigned i) - { -#ifndef NDEBUG - if (i >= size()) - { - std::stringstream msg; - msg << "Out of bounds: " << i << " >= " << size(); - KTERROR(utillog_physarr, msg.str()); - throw std::out_of_range(msg.str()); - } -#endif - return fData[i]; - } - - inline bool KTPhysicalArray< 1, std::complex >::IsCompatibleWith(const KTPhysicalArray< 1, std::complex >& rhs) const - { - //return (this->size() == rhs.size() && this->GetRangeMin() == rhs.GetRangeMin() && this->GetRangeMax() == GetRangeMax()); - return (this->size() == rhs.size()); - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator+=(const KTPhysicalArray< 1, std::complex>& rhs) - { - if (this->IsCompatibleWith(rhs)) - { - fData += rhs.fData; - } - - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator-=(const KTPhysicalArray< 1, std::complex>& rhs) - { - if (this->IsCompatibleWith(rhs)) - { - fData -= rhs.fData; - } - - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator*=(const KTPhysicalArray< 1, std::complex>& rhs) - { - if (this->IsCompatibleWith(rhs)) - { - fData *= rhs.fData; - } - - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator/=(const KTPhysicalArray< 1, std::complex>& rhs) - { - if (this->IsCompatibleWith(rhs)) - { - fData /= rhs.fData; - } - - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator+=(const std::complex& rhs) - { - fData += rhs; - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator-=(const std::complex& rhs) - { - fData -= rhs; - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator*=(const std::complex& rhs) - { - fData *= rhs; - return *this; - } - - inline KTPhysicalArray< 1, std::complex >& KTPhysicalArray< 1, std::complex >::operator/=(const std::complex& rhs) - { - fData /= rhs; - return *this; - } - - //************************* - // Operator implementations - //************************* - - /// Add two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - inline KTPhysicalArray< 1, std::complex > operator+(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs += rhs; - return lhs; - } - - /// Subtracts two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - inline KTPhysicalArray< 1, std::complex > operator-(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs -= rhs; - return lhs; - } - - /// Multiplies two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - inline KTPhysicalArray< 1, std::complex > operator*(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs *= rhs; - return lhs; - } - - /// Divides two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - inline KTPhysicalArray< 1, std::complex > operator/(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs /= rhs; - return lhs; - } - - std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 1, std::complex >& rhs) - { - ostr << rhs.GetData(); - return ostr; - } - - //************************* - // Iterators - //************************* - - inline KTPhysicalArray< 1, std::complex >::const_iterator KTPhysicalArray< 1, std::complex >::begin() const - { - return fData.data(); - } - - inline KTPhysicalArray< 1, std::complex >::const_iterator KTPhysicalArray< 1, std::complex >::end() const - { - return fData.data() + size(); - } - - inline KTPhysicalArray< 1, std::complex >::iterator KTPhysicalArray< 1, std::complex >::begin() - { - return fData.data(); - } - - inline KTPhysicalArray< 1, std::complex >::iterator KTPhysicalArray< 1, std::complex >::end() - { - return fData.data() + size(); - } - - inline KTPhysicalArray< 1, std::complex >::const_reverse_iterator KTPhysicalArray< 1, std::complex >::rbegin() const - { - //a reverse iterator points to the first element outside the container - return std::reverse_iterator* >(fData.data() + size()); - } - - inline KTPhysicalArray< 1, std::complex >::const_reverse_iterator KTPhysicalArray< 1, std::complex >::rend() const - { - return std::reverse_iterator* >(fData.data()); - } - - inline KTPhysicalArray< 1, std::complex >::reverse_iterator KTPhysicalArray< 1, std::complex >::rbegin() - { - return std::reverse_iterator< std::complex* >(fData.data() + size()); - } - - inline KTPhysicalArray< 1, std::complex >::reverse_iterator KTPhysicalArray< 1, std::complex >::rend() - { - return std::reverse_iterator< std::complex* >(fData.data()); - } - //************************* - // Begin, end free functions for range based for-loop + // Free function operators //************************* - KTPhysicalArray< 1, std::complex >::iterator begin(KTPhysicalArray< 1, std::complex > & array) - { - return array.begin(); - } + KTPhysicalArray< 1, std::complex > operator+(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); + KTPhysicalArray< 1, std::complex > operator-(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); + KTPhysicalArray< 1, std::complex > operator*(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); + KTPhysicalArray< 1, std::complex > operator/(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); + std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 1, std::complex >& rhs); - KTPhysicalArray< 1, std::complex >::iterator end(KTPhysicalArray< 1, std::complex > & array) - { - return array.end(); - } - - KTPhysicalArray< 1, std::complex >::const_iterator begin(const KTPhysicalArray< 1, std::complex > & array) - { - return array.begin(); - } - - KTPhysicalArray< 1, std::complex >::const_iterator end(const KTPhysicalArray< 1, std::complex > & array) - { - return array.end(); - } - //************************* // 2-D array implementation //************************* @@ -469,475 +244,24 @@ namespace Katydid }; - size_t KTPhysicalArray< 2, std::complex >::cols() const - { - return static_cast(fData.cols()); - } - - size_t KTPhysicalArray< 2, std::complex >::rows() const - { - return static_cast(fData.rows()); - } - - - inline const Eigen::Array< std::complex, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor >& KTPhysicalArray< 2, std::complex >::GetData() const - { - return fData; - } - - - inline typename Eigen::Array< std::complex, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor >& KTPhysicalArray< 2, std::complex >::GetData() - { - return fData; - } - - - inline const std::string& KTPhysicalArray< 2, std::complex >::GetDataLabel() const - { - return fLabel; - } - - - inline void KTPhysicalArray< 2, std::complex >::SetDataLabel(const std::string& label) - { - fLabel = label; - return; - } - - - inline const std::complex& KTPhysicalArray< 2, std::complex >::operator()(unsigned i, unsigned j) const - { -#ifdef Katydid_DEBUG - if (i >= size(1)) - { - std::stringstream msg; - msg << "Out of bounds on axis 1: " << i << " >= " << size(1); - KTERROR(utillog_physarr, msg.str()); - throw std::out_of_range(msg.str()); - } - if (j >= size(2)) - { - std::stringstream msg; - msg << "Out of bounds on axis 2: " << j << " >= " << size(2); - KTERROR(utillog_physarr, msg.str()); - throw std::out_of_range(msg.str()); - } -#endif - return fData(i,j); - } - - - inline std::complex& KTPhysicalArray< 2, std::complex >::operator()(unsigned i, unsigned j) - { -#ifdef Katydid_DEBUG - if (i >= size(1)) - { - std::stringstream msg; - msg << "Out of bounds on axis 1: " << i << " >= " << size(1); - KTERROR(utillog_physarr, msg.str()); - throw std::out_of_range(msg.str()); - } - if (j >= size(2)) - { - std::stringstream msg; - msg << "Out of bounds on axis 2: " << j << " >= " << size(2); - KTERROR(utillog_physarr, msg.str()); - throw std::out_of_range(msg.str()); - } -#endif - return fData(i,j); - } - - - inline bool KTPhysicalArray< 2, std::complex >::IsCompatibleWith(const KTPhysicalArray< 2, std::complex >& rhs) const - { - //return (this->size(1) == rhs.size(1) && this->GetRangeMin(1) == rhs.GetRangeMin(1) && this->GetRangeMax(1) == GetRangeMax(1) && - // this->size(2) == rhs.size(2) && this->GetRangeMin(2) == rhs.GetRangeMin(2) && this->GetRangeMax(2) == GetRangeMax(2)); - return (this->size(1) == rhs.size(1) && - this->size(2) == rhs.size(2)); - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator+=(const KTPhysicalArray< 2, std::complex>& rhs) - { - if ( this->IsCompatibleWith(rhs) ) - { - fData += rhs.fData; - } - - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator-=(const KTPhysicalArray< 2, std::complex>& rhs) - { - if ( this->IsCompatibleWith(rhs) ) - { - fData -= rhs.fData; - } - - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator*=(const KTPhysicalArray< 2, std::complex>& rhs) - { - if ( this->IsCompatibleWith(rhs) ) - { - fData *= rhs.fData; - } - - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator/=(const KTPhysicalArray< 2, std::complex>& rhs) - { - if ( this->IsCompatibleWith(rhs) ) - { - fData /= rhs.fData; - } - - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator+=(const std::complex& rhs) - { - fData += rhs; - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator-=(const std::complex& rhs) - { - fData -= rhs; - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator*=(const std::complex& rhs) - { - fData *= rhs; - return *this; - } - - - KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator/=(const std::complex& rhs) - { - fData /= rhs; - return *this; - } - - - inline KTPhysicalArray< 2, std::complex >::const_row_iterator KTPhysicalArray< 2, std::complex >::begin2() const - { - return const_row_iterator{ fData.data(), 1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::const_col_iterator KTPhysicalArray< 2, std::complex >::begin1() const - { - return const_col_iterator{ fData.data(), static_cast(size(1)) }; - } - - - inline KTPhysicalArray< 2, std::complex >::const_row_iterator KTPhysicalArray< 2, std::complex >::end2() const - { - return const_row_iterator{ fData.data() + size(1), 1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::const_col_iterator KTPhysicalArray< 2, std::complex >::end1() const - { - return const_col_iterator{ fData.data() + size(1)*size(2), static_cast(size(1)) }; - } - - - inline KTPhysicalArray< 2, std::complex >::row_iterator KTPhysicalArray< 2, std::complex >::begin2() - { - return row_iterator{ fData.data(), 1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::col_iterator KTPhysicalArray< 2, std::complex >::begin1() - { - return col_iterator{ fData.data(), static_cast(size(1)) }; - } - - - inline KTPhysicalArray< 2, std::complex >::row_iterator KTPhysicalArray< 2, std::complex >::end2() - { - return row_iterator{ fData.data() + size(1), 1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::col_iterator KTPhysicalArray< 2, std::complex >::end1() - { - return col_iterator{ fData.data() + size(1)*size(2), static_cast(size(1)) }; - } - //************************* - // Reverse iterators + // Free function operators 2D //************************* - - inline KTPhysicalArray< 2, std::complex >::const_reverse_row_iterator KTPhysicalArray< 2, std::complex >::rbegin2() const - { - return const_reverse_row_iterator{ fData.data() + size(1), -1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::const_reverse_col_iterator KTPhysicalArray< 2, std::complex >::rbegin1() const - { - return const_reverse_col_iterator{ fData.data() + size(1)*size(2), -static_cast(size(1)) }; - } - - - inline KTPhysicalArray< 2, std::complex >::const_reverse_row_iterator KTPhysicalArray< 2, std::complex >::rend2() const - { - return const_reverse_row_iterator{ fData.data(), -1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::const_reverse_col_iterator KTPhysicalArray< 2, std::complex >::rend1() const - { - return const_reverse_col_iterator{ fData.data(), -static_cast(size(1)) }; - } - - - inline KTPhysicalArray< 2, std::complex >::reverse_row_iterator KTPhysicalArray< 2, std::complex >::rbegin2() - { - return reverse_row_iterator{ fData.data() + size(1), -1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::reverse_col_iterator KTPhysicalArray< 2, std::complex >::rbegin1() - { - return reverse_col_iterator{ fData.data() + size(1)*size(2), -static_cast(size(1)) }; - } - - - inline KTPhysicalArray< 2, std::complex >::reverse_row_iterator KTPhysicalArray< 2, std::complex >::rend2() - { - return reverse_row_iterator{ fData.data(), -1 }; - } - - - inline KTPhysicalArray< 2, std::complex >::reverse_col_iterator KTPhysicalArray< 2, std::complex >::rend1() - { - return reverse_col_iterator{ fData.data(), -static_cast(size(1)) }; - } - - - std::complex KTPhysicalArray< 2, std::complex >::GetMaximumBin(unsigned& maxXBin, unsigned& maxYBin) const - { - //return fData.maxCoeff(&maxXBin, &maxYBin); // can be used for real values - - //functor for finding the maximum - struct MaxVisitor - { - unsigned row; - unsigned col; - value_type max; - - void init( const value_type& value, Eigen::Index i, Eigen::Index j ) - { - max=value; - row=i; - col=j; - }; - - void operator()( const value_type& value, Eigen::Index i, Eigen::Index j ) - { - if(std::abs(value)>std::abs(max)) - { - max=value; - row=i; - col=j; - } - }; - }; - - MaxVisitor mv; - fData.visit(mv); - maxXBin = mv.row; - maxYBin = mv.col; - - return mv.max; - } - - - std::complex KTPhysicalArray< 2, std::complex >::GetMinimumBin(unsigned& minXBin, unsigned& minYBin) const - { - //return fData.minCoeff(&minXBin, &minYBin); // can be used for real values - - //functor for finding the minimum - struct MinVisitor - { + KTPhysicalArray< 2, std::complex > operator+(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); + KTPhysicalArray< 2, std::complex > operator-(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); + KTPhysicalArray< 2, std::complex > operator*(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); + KTPhysicalArray< 2, std::complex > operator/(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); + std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 2, std::complex >& rhs); - unsigned row; - unsigned col; - value_type min; - - void init( const value_type& value, Eigen::Index i, Eigen::Index j ) - { - min=value; - row=i; - col=j; - }; - - void operator()( const value_type& value, Eigen::Index i, Eigen::Index j ) - { - if(std::abs(value), std::complex > KTPhysicalArray< 2, std::complex >::GetMinMaxBin(unsigned& minXBin, unsigned& minYBin, unsigned& maxXBin, unsigned& maxYBin) - { - //functor for finding minimum and maximum at the same time - struct MinMaxVisitor - { - - unsigned rowMax; - unsigned colMax; - unsigned rowMin; - unsigned colMin; - value_type max; - value_type min; - - void init( const value_type& value, Eigen::Index i, Eigen::Index j ) - { - max=value; - min=value; - rowMax=i; - colMax=j; - rowMin=i; - colMin=j; - }; - - void operator()( const value_type& value, Eigen::Index i, Eigen::Index j ) - { - if(std::abs(value)>std::abs(max)) - { - max=value; - rowMax=i; - colMax=j; - } - if(std::abs(value) > operator+(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs += rhs; - return lhs; - } - - /// Subtracts two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator-(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs -= rhs; - return lhs; - } - /// Multiplies two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator*(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs *= rhs; - return lhs; - } - - /// Divides two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator/(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs /= rhs; - return lhs; - } - - std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 2, std::complex >& rhs) - { - ostr << rhs.GetData(); - return ostr; - } - - //************************* - // Matrix and Vector operations + // Free function operators for matrix and vector operations //************************* - inline KTPhysicalArray< 2, std::complex >& KTPhysicalArray< 2, std::complex >::operator%=(const KTPhysicalArray< 2, std::complex >& rhs) - { - - fData = fData.matrix()*rhs.fData.matrix(); - - SetRangeMin(2, rhs.GetRangeMin(2)); - SetRangeMax(2, rhs.GetRangeMax(2)); - - return *this; - } + KTPhysicalArray< 2, std::complex > operator%(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); + KTPhysicalArray< 1, std::complex > operator%(const KTPhysicalArray< 2, std::complex >& lhs, KTPhysicalArray< 1, std::complex > rhs); - inline KTPhysicalArray< 2, std::complex > operator%(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - - lhs %= rhs; - return lhs; - } - //matrix-vector-multiplication - inline KTPhysicalArray< 1, std::complex > operator%(const KTPhysicalArray< 2, std::complex >& lhs, KTPhysicalArray< 1, std::complex > rhs) - { - - rhs.fData = lhs.fData.matrix()*rhs.fData.matrix(); - - rhs.SetRangeMin(lhs.GetRangeMin(1)); - rhs.SetRangeMax(lhs.GetRangeMax(1)); - - return rhs; - } - - } /* namespace Katydid */ #endif /* KTPHYSICALARRAYCOMPLEX_HH_ */ From ec9a277ca242a986aa02085e5180e28a07f62b38 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Fri, 16 Apr 2021 16:31:50 +0200 Subject: [PATCH 03/48] Modify KTMultiFSDataFFTW to use Eigen --- Source/Data/Transform/KTMultiFSDataFFTW.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Data/Transform/KTMultiFSDataFFTW.cc b/Source/Data/Transform/KTMultiFSDataFFTW.cc index eca38f508..247d33624 100644 --- a/Source/Data/Transform/KTMultiFSDataFFTW.cc +++ b/Source/Data/Transform/KTMultiFSDataFFTW.cc @@ -83,7 +83,7 @@ namespace Katydid if (fs == NULL) continue; for (int iBinY=1; iBinY<=hist->GetNbinsY(); ++iBinY) { - hist->SetBinContent(iBinX, iBinY, sqrt((*fs)(iBinY-1)[0] * (*fs)(iBinY-1)[0] + (*fs)(iBinY-1)[1] * (*fs)(iBinY-1)[1])); + hist->SetBinContent(iBinX, iBinY, std::abs((*fs)(iBinY-1))); } } @@ -110,7 +110,7 @@ namespace Katydid if (fs == NULL) continue; for (int iBinY=1; iBinY<=hist->GetNbinsY(); ++iBinY) { - hist->SetBinContent(iBinX, iBinY, atan2((*fs)(iBinY-1)[1], (*fs)(iBinY-1)[0])); + hist->SetBinContent(iBinX, iBinY, std::arg((*fs)(iBinY-1))); } } @@ -138,7 +138,7 @@ namespace Katydid if (fs == NULL) continue; for (int iBinY=1; iBinY<=hist->GetNbinsY(); ++iBinY) { - value = (*fs)(iBinY-1)[0] * (*fs)(iBinY-1)[0] + (*fs)(iBinY-1)[1] * (*fs)(iBinY-1)[1]; + value = std::norm((*fs)(iBinY-1)); hist->SetBinContent(iBinX, iBinY, value); } } From b5efc90af1966fdad1909ad488d5023e15bd461a Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Fri, 16 Apr 2021 16:32:39 +0200 Subject: [PATCH 04/48] Modify KT2ROOT.cc Modifies the functions accepting KTFrequencySpectrumFFTW types --- Source/IO/Conversions/KT2ROOT.cc | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Source/IO/Conversions/KT2ROOT.cc b/Source/IO/Conversions/KT2ROOT.cc index d40dbb543..8beaa847d 100644 --- a/Source/IO/Conversions/KT2ROOT.cc +++ b/Source/IO/Conversions/KT2ROOT.cc @@ -392,7 +392,7 @@ namespace Katydid TH1D* hist = new TH1D(name.c_str(), "Frequency Spectrum: Magnitude", (int) nBins, fs->GetRangeMin(), fs->GetRangeMax()); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, std::sqrt((*fs)(iBin)[0] * (*fs)(iBin)[0] + (*fs)(iBin)[1] * (*fs)(iBin)[1])); + hist->SetBinContent((int) iBin + 1, fs->GetAbs(iBin)); } hist->SetXTitle("Frequency (Hz)"); hist->SetYTitle("Voltage (V)"); @@ -405,7 +405,7 @@ namespace Katydid TH1D* hist = new TH1D(name.c_str(), "Frequency Spectrum: Phase", (int) nBins, fs->GetRangeMin(), fs->GetRangeMax()); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, std::atan2((*fs)(iBin)[1], (*fs)(iBin)[0])); + hist->SetBinContent((int) iBin + 1, fs->GetArg(iBin)); } hist->SetXTitle("Frequency (Hz)"); hist->SetYTitle("Phase"); @@ -421,7 +421,7 @@ namespace Katydid for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, scaling * ((*fs)(iBin)[0] * (*fs)(iBin)[0] + (*fs)(iBin)[1] * (*fs)(iBin)[1])); + hist->SetBinContent((int) iBin + 1, scaling * fs->GetAbs(iBin)); } hist->SetXTitle("Frequency (Hz)"); @@ -438,7 +438,7 @@ namespace Katydid // skip the DC bin; start at iBin = 1 for (unsigned iBin = 1; iBin < nBins; ++iBin) { - value = std::sqrt((*fs)(iBin)[0] * (*fs)(iBin)[0] + (*fs)(iBin)[1] * (*fs)(iBin)[1]); + value = fs->GetAbs(iBin); if (value < tMinMag) tMinMag = value; if (value > tMaxMag) tMaxMag = value; } @@ -446,7 +446,7 @@ namespace Katydid TH1D* hist = new TH1D(name.c_str(), "Magnitude Distribution", 100, tMinMag * 0.95, tMaxMag * 1.05); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - value = sqrt((*fs)(iBin)[0] * (*fs)(iBin)[0] + (*fs)(iBin)[1] * (*fs)(iBin)[1]); + value = fs->GetAbs(iBin); hist->Fill(value); } hist->SetXTitle("Voltage (V)"); @@ -463,7 +463,7 @@ namespace Katydid // skip the DC bin; start at iBin = 1 for (unsigned iBin = 1; iBin < nBins; ++iBin) { - value = ((*fs)(iBin)[0] * (*fs)(iBin)[0] + (*fs)(iBin)[1] * (*fs)(iBin)[1]) * scaling; + value = fs->GetAbs(iBin) * scaling; if (value < tMinMag) tMinMag = value; if (value > tMaxMag) tMaxMag = value; } @@ -471,7 +471,7 @@ namespace Katydid TH1D* hist = new TH1D(name.c_str(), "Power Distribution", 100, tMinMag * 0.95, tMaxMag * 1.05); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - value = (*fs)(iBin)[0] * (*fs)(iBin)[0] + (*fs)(iBin)[1] * (*fs)(iBin)[1]; + value = fs->GetAbs(iBin); hist->Fill(value * scaling); } hist->SetXTitle("Power (W)"); From f5d3ce73611a35ac3284d4010438da7bebf18f70 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Fri, 16 Apr 2021 22:52:18 +0200 Subject: [PATCH 05/48] Add functionality to KTPhysicalArrayComplex Modifications of the arithmetic operators to enable full compatibility for child classes --- Source/Utility/KTPhysicalArrayComplex.cc | 76 ------------------ Source/Utility/KTPhysicalArrayComplex.hh | 98 ++++++++++++++++++++++-- 2 files changed, 90 insertions(+), 84 deletions(-) diff --git a/Source/Utility/KTPhysicalArrayComplex.cc b/Source/Utility/KTPhysicalArrayComplex.cc index cc4f74ad4..b74271a46 100644 --- a/Source/Utility/KTPhysicalArrayComplex.cc +++ b/Source/Utility/KTPhysicalArrayComplex.cc @@ -165,42 +165,6 @@ namespace Katydid // Operator implementations //************************* - /// Add two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - KTPhysicalArray< 1, std::complex > operator+(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs += rhs; - return lhs; - } - - /// Subtracts two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - KTPhysicalArray< 1, std::complex > operator-(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs -= rhs; - return lhs; - } - - /// Multiplies two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - KTPhysicalArray< 1, std::complex > operator*(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs *= rhs; - return lhs; - } - - /// Divides two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - KTPhysicalArray< 1, std::complex > operator/(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 1, std::complex >(); - - lhs /= rhs; - return lhs; - } - std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 1, std::complex >& rhs) { ostr << rhs.GetData(); @@ -713,46 +677,6 @@ namespace Katydid //************************* // Operators //************************* - - /// Add two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator+(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs += rhs; - return lhs; - } - - /// Subtracts two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator-(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs -= rhs; - return lhs; - } - - /// Multiplies two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator*(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs *= rhs; - return lhs; - } - - /// Divides two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. - - KTPhysicalArray< 2, std::complex > operator/(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs) - { - if (! lhs.IsCompatibleWith(rhs)) return KTPhysicalArray< 2, std::complex >(); - - lhs /= rhs; - return lhs; - } std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 2, std::complex >& rhs) { diff --git a/Source/Utility/KTPhysicalArrayComplex.hh b/Source/Utility/KTPhysicalArrayComplex.hh index dca0a0e07..37e8cce77 100644 --- a/Source/Utility/KTPhysicalArrayComplex.hh +++ b/Source/Utility/KTPhysicalArrayComplex.hh @@ -16,6 +16,7 @@ #include #include #include +#include #include @@ -141,12 +142,52 @@ namespace Katydid // Free function operators //************************* - KTPhysicalArray< 1, std::complex > operator+(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); - KTPhysicalArray< 1, std::complex > operator-(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); - KTPhysicalArray< 1, std::complex > operator*(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); - KTPhysicalArray< 1, std::complex > operator/(KTPhysicalArray< 1, std::complex > lhs, const KTPhysicalArray< 1, std::complex >& rhs); std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 1, std::complex >& rhs); + // this template enables usability of the arithmetic operators with all derived classes + template + using IsKTPhysicalArray1D = typename std::enable_if >, T>::value, T>::type; + + /// Add two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + template < typename T> + IsKTPhysicalArray1D operator+(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T();//KTPhysicalArray< 1, std::complex >(); + + lhs.KTPhysicalArray<1, std::complex >::operator+=(rhs); + return lhs; + } + + /// Subtracts two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + template < typename T> + IsKTPhysicalArray1D operator-(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<1, std::complex >::operator-=(rhs); + return lhs; + } + + /// Multiplies two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + template < typename T> + IsKTPhysicalArray1D operator*(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<1, std::complex >::operator*=(rhs); + return lhs; + } + + /// Divides two 1-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + template < typename T> + IsKTPhysicalArray1D operator/(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<1, std::complex >::operator/=(rhs); + return lhs; + } + //************************* // 2-D array implementation //************************* @@ -248,12 +289,53 @@ namespace Katydid // Free function operators 2D //************************* - KTPhysicalArray< 2, std::complex > operator+(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); - KTPhysicalArray< 2, std::complex > operator-(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); - KTPhysicalArray< 2, std::complex > operator*(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); - KTPhysicalArray< 2, std::complex > operator/(KTPhysicalArray< 2, std::complex > lhs, const KTPhysicalArray< 2, std::complex >& rhs); std::ostream& operator<< (std::ostream& ostr, const KTPhysicalArray< 2, std::complex >& rhs); + template + using IsKTPhysicalArray2D = typename std::enable_if >, T>::value, T>::type; + + /// Add two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + + template < typename T> + IsKTPhysicalArray2D operator+(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<2, std::complex >::operator+=(rhs); + return lhs; + } + + /// Subtracts two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + + template < typename T> + IsKTPhysicalArray2D operator-(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<2, std::complex >::operator-=(rhs); + return lhs; + } + + /// Multiplies two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + + template < typename T> + IsKTPhysicalArray2D operator*(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<2, std::complex >::operator*=(rhs); + return lhs; + } + + /// Divides two 2-D KTPhysicalArrays; requires lhs.size() == rhs.size(); axis range set to that of lhs. + template < typename T> + IsKTPhysicalArray2D operator/(T lhs, const T& rhs) + { + if (! lhs.IsCompatibleWith(rhs)) return T(); + + lhs.KTPhysicalArray<2, std::complex >::operator/=(rhs); + return lhs; + } //************************* // Free function operators for matrix and vector operations From 19a216b74bf615b9785b75e183d2b6ad832d3faf Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Mon, 19 Apr 2021 20:06:11 +0200 Subject: [PATCH 06/48] Add new function to KTFrequencySpectrumFFTW The function returns the norm of a bin, i.e. the square of the absolute value without any unnecessary calculations of squareroots. --- Source/Data/Transform/KTFrequencySpectrumFFTW.hh | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh index 4092683f9..26fb2e3fc 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh @@ -60,6 +60,7 @@ namespace Katydid virtual double GetAbs(unsigned bin) const; virtual double GetArg(unsigned bin) const; + virtual double GetNorm(unsigned bin) const; virtual void SetPolar(unsigned bin, double abs, double arg); @@ -188,6 +189,13 @@ namespace Katydid //return sqrt((*fPointCache)[0]*(*fPointCache)[0] + (*fPointCache)[1]*(*fPointCache)[1]); return std::abs((*this)(bin)); } + + inline double KTFrequencySpectrumFFTW::GetNorm(unsigned bin) const + { + //fPointCache = &(*this)(bin); + //return sqrt((*fPointCache)[0]*(*fPointCache)[0] + (*fPointCache)[1]*(*fPointCache)[1]); + return std::norm((*this)(bin)); + } inline double KTFrequencySpectrumFFTW::GetArg(unsigned bin) const { From 5940f666d32b09eb43c78c6faf9e038ee8e76db4 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Mon, 19 Apr 2021 21:04:13 +0200 Subject: [PATCH 07/48] Add comment to KTFrequencySpectrumFFTW.hh --- Source/Data/Transform/KTFrequencySpectrumFFTW.hh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh index 26fb2e3fc..668e62e65 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh @@ -179,6 +179,9 @@ namespace Katydid //~ (*const_cast< fftw_complex* >(fPointCache))[0] = real; //~ (*const_cast< fftw_complex* >(fPointCache))[1] = imag; + //having the point cache for this really is not important + // for >=O2 the assembly code should be exactly the same + // can check it in https://godbolt.org/ (*this)(bin) = std::complex(real, imag); return; } From 31a10d18fce0ec5ccbbb2aeac4458389e34f2276 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 19:44:06 +0200 Subject: [PATCH 08/48] Modify FFT processors The FFT processors now take the modified KTFrequencySpectrumFFTW. The call to the FFTW routine is now done with the data pointer of the Eigen object which is cast to be a fftw_complex* --- Source/Transform/KTForwardFFTW.cc | 6 ++--- Source/Transform/KTFractionalFFT.cc | 34 +++++++++++++++++------------ Source/Transform/KTReverseFFTW.cc | 14 ++++++++++-- 3 files changed, 35 insertions(+), 19 deletions(-) diff --git a/Source/Transform/KTForwardFFTW.cc b/Source/Transform/KTForwardFFTW.cc index 35eac0f0e..a18cdcf7a 100644 --- a/Source/Transform/KTForwardFFTW.cc +++ b/Source/Transform/KTForwardFFTW.cc @@ -473,7 +473,7 @@ namespace Katydid void KTForwardFFTW::DoTransform(const KTTimeSeriesReal* tsIn, KTFrequencySpectrumFFTW* fsOut) const { std::copy(tsIn->begin(), tsIn->end(), fRInputArray); - fftw_execute_dft_r2c(fForwardPlan, fRInputArray, fsOut->GetData()); + fftw_execute_dft_r2c(fForwardPlan, fRInputArray, reinterpret_cast(fsOut->GetData().data())); (*fsOut) *= sqrt(2. / (double)fTimeSize); return; } @@ -510,7 +510,7 @@ namespace Katydid fCInputArray[iBin][0] = tsIn->GetData()[iBin]; fCInputArray[iBin][1] = 0; } - fftw_execute_dft(fForwardPlan, fCInputArray, fsOut->GetData()); + fftw_execute_dft(fForwardPlan, fCInputArray, reinterpret_cast(fsOut->GetData().data())); (*fsOut) *= sqrt(1. / (double)fTimeSize); return; } @@ -542,7 +542,7 @@ namespace Katydid void KTForwardFFTW::DoTransform(const KTTimeSeriesFFTW* tsIn, KTFrequencySpectrumFFTW* fsOut) const { - fftw_execute_dft(fForwardPlan, tsIn->GetData(), fsOut->GetData()); + fftw_execute_dft(fForwardPlan, tsIn->GetData(), reinterpret_cast(fsOut->GetData().data())); (*fsOut) *= sqrt(1. / (double) fTimeSize); return; } diff --git a/Source/Transform/KTFractionalFFT.cc b/Source/Transform/KTFractionalFFT.cc index e464d743a..4d930ab4f 100644 --- a/Source/Transform/KTFractionalFFT.cc +++ b/Source/Transform/KTFractionalFFT.cc @@ -184,19 +184,26 @@ namespace Katydid fs->SetNTimeBins( slice.GetSliceSize() ); // Second chirp transform - for( unsigned iBin = 0; iBin < fs->GetNFrequencyBins(); ++iBin ) - { - // Obtain norm and phase from components - norm = std::hypot( (*fs)(iBin)[0], (*fs)(iBin)[1] ); - phase = atan2( (*fs)(iBin)[1], (*fs)(iBin)[0] ); - - // Shift phase - phase -= q2 * (double)iBin * (double)iBin; + + //~ for( unsigned iBin = 0; iBin < fs->GetNFrequencyBins(); ++iBin ) + //~ { + //~ // Obtain norm and phase from components + //~ norm = std::hypot( (*fs)(iBin)[0], (*fs)(iBin)[1] ); + //~ phase = atan2( (*fs)(iBin)[1], (*fs)(iBin)[0] ); + + //~ // Shift phase + //~ phase -= q2 * (double)iBin * (double)iBin; + + //~ // Assign components from norm and new phase + //~ (*fs)(iBin)[0] = norm * cos( phase ); + //~ (*fs)(iBin)[1] = norm * sin( phase ); + //~ } + + // Shift phase + auto bins = Eigen::ArrayXd::LinSpaced(fs->GetNFrequencyBins(), 0, + fs->GetNFrequencyBins()-1); - // Assign components from norm and new phase - (*fs)(iBin)[0] = norm * cos( phase ); - (*fs)(iBin)[1] = norm * sin( phase ); - } + fs->GetData() *= exp(-std::complex{0.0, 1.0}*q2*bins*bins); // Reverse FFT KTTimeSeriesFFTW* newTS = fReverseFFT.TransformToComplex( fs ); @@ -218,8 +225,7 @@ namespace Katydid // Assign components from norm and new phase (*newTS)(iBin)[0] = norm * cos( phase ); (*newTS)(iBin)[1] = norm * sin( phase ); - (*newFS)(iBin)[0] = norm * cos( phase ); - (*newFS)(iBin)[1] = norm * sin( phase ); + (*newFS)(iBin) = std::polar(norm, phase); } newTSData.SetTimeSeries( newTS, iComponent ); // newTS now owned by newTSData diff --git a/Source/Transform/KTReverseFFTW.cc b/Source/Transform/KTReverseFFTW.cc index 1576078ec..68c49b655 100644 --- a/Source/Transform/KTReverseFFTW.cc +++ b/Source/Transform/KTReverseFFTW.cc @@ -278,7 +278,12 @@ namespace Katydid void KTReverseFFTW::DoTransform(const KTFrequencySpectrumFFTW* fsIn, KTTimeSeriesReal* tsOut) const { - fftw_execute_dft_c2r(fReversePlan, fsIn->GetData(), fROutputArray); + fftw_complex *data = + const_cast( + reinterpret_cast( + fsIn->GetData().data())); + + fftw_execute_dft_c2r(fReversePlan, data, fROutputArray); std::copy(fROutputArray, fROutputArray+fTimeSize, tsOut->begin()); (*tsOut) *= sqrt(1. / double(fTimeSize)); return; @@ -360,7 +365,12 @@ namespace Katydid void KTReverseFFTW::DoTransform(const KTFrequencySpectrumFFTW* fsIn, KTTimeSeriesFFTW* tsOut) const { - fftw_execute_dft(fReversePlan, fsIn->GetData(), tsOut->GetData()); + + fftw_complex *data = + const_cast( + reinterpret_cast( + fsIn->GetData().data())); + fftw_execute_dft(fReversePlan, data, tsOut->GetData()); (*tsOut) *= sqrt(1. / double(fTimeSize)); return; } From 35181ffe181ee884a7e9526129c5813ead8e8146 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 19:48:46 +0200 Subject: [PATCH 09/48] Modify files depending on KTFrequencySpectrumFFTW Mostly modifies a few lines per file that did applied operations specifig to fftw_complex --- .../KTFrequencyCandidateIdentifier.cc | 2 +- .../Validation/TestChannelAggregator.cc | 6 +- .../Validation/TestConvolution1D.cc | 6 +- .../Executables/Validation/TestCorrelator.cc | 77 ++++++++++--------- .../Executables/Validation/TestForwardFFTW.cc | 8 +- .../Executables/Validation/TestReverseFFTW.cc | 5 +- .../Validation/TestSpectrogramStriper.cc | 4 +- .../Validation/TestSpectrumDiscriminator.cc | 11 ++- .../Executables/Validation/TestWignerVille.cc | 2 +- .../KTAmplitudeDistributor.cc | 4 +- Source/SpectrumAnalysis/KTAxialAggregator.cc | 3 +- .../SpectrumAnalysis/KTChannelAggregator.cc | 3 +- Source/SpectrumAnalysis/KTConvolution.cc | 22 +++--- Source/SpectrumAnalysis/KTDataAccumulator.cc | 16 ++-- .../SpectrumAnalysis/KTGainNormalization.cc | 15 ++-- .../KTGainVariationProcessor.cc | 3 +- .../SpectrumAnalysis/KTSpectrogramStriper.cc | 4 +- .../KTSpectrumDiscriminator.cc | 4 +- .../KTVariableSpectrumDiscriminator.cc | 6 +- 19 files changed, 100 insertions(+), 101 deletions(-) diff --git a/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc b/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc index 564510275..d6028b731 100644 --- a/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc +++ b/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc @@ -150,7 +150,7 @@ namespace Katydid double peakValue = 0.; for (unsigned iBin=firstBin; iBin <= lastBin; iBin++) { - value = sqrt((*freqSpec)(iBin)[0] * (*freqSpec)(iBin)[0] + (*freqSpec)(iBin)[1] * (*freqSpec)(iBin)[1]); + value = freqSpec->GetAbs(iBin); // sqrt((*freqSpec)(iBin)[0] * (*freqSpec)(iBin)[0] + (*freqSpec)(iBin)[1] * (*freqSpec)(iBin)[1]); weightedMean += freqSpec->GetBinCenter(iBin) * value; integral += value; if (value > peakValue) peakValue = value; diff --git a/Source/Executables/Validation/TestChannelAggregator.cc b/Source/Executables/Validation/TestChannelAggregator.cc index 4e5cf0bfc..7bb42bd45 100644 --- a/Source/Executables/Validation/TestChannelAggregator.cc +++ b/Source/Executables/Validation/TestChannelAggregator.cc @@ -79,13 +79,11 @@ int main() //Arbitrarily chose the middle bin for the signal frequency if(iFreqBin==(int)(nFreqBins/2)) { - (*fftwSpectrum)(iFreqBin)[0]= cos(phaseShift); - (*fftwSpectrum)(iFreqBin)[1]= -sin(phaseShift); + fftwSpectrum->SetRect(iFreqBin, cos(phaseShift), -sin(phaseShift)); } // Since non-signal frequencies are irrelavant, manually setthing them to 0 else{ - (*fftwSpectrum)(iFreqBin)[0]= 0.0; - (*fftwSpectrum)(iFreqBin)[1]= 0.0; + fftwSpectrum->SetRect(iFreqBin, 0.0, 0.0); } } newFreqData.SetSpectrum(fftwSpectrum, iComponent); diff --git a/Source/Executables/Validation/TestConvolution1D.cc b/Source/Executables/Validation/TestConvolution1D.cc index 12b980b47..fe54be4e5 100644 --- a/Source/Executables/Validation/TestConvolution1D.cc +++ b/Source/Executables/Validation/TestConvolution1D.cc @@ -43,15 +43,13 @@ int main() if( powerSpect->GetBinCenter( iBin ) >= pulseStart && powerSpect->GetBinCenter( iBin ) <= pulseEnd ) { (*powerSpect)(iBin) = 1.0; - (*fftwSpect)(iBin)[0] = 1.0; - (*fftwSpect)(iBin)[1] = 0.0; + fftwSpect->SetRect(iBin, 1.0, 0.0); polarSpect->SetRect( iBin, 1.0, 0.0 ); } else { (*powerSpect)(iBin) = 0.0; - (*fftwSpect)(iBin)[0] = 0.0; - (*fftwSpect)(iBin)[1] = 0.0; + fftwSpect->SetRect( iBin, 0.0, 0.0); polarSpect->SetRect( iBin, 0.0, 0.0 ); } diff --git a/Source/Executables/Validation/TestCorrelator.cc b/Source/Executables/Validation/TestCorrelator.cc index 9535c019b..99bd662e9 100644 --- a/Source/Executables/Validation/TestCorrelator.cc +++ b/Source/Executables/Validation/TestCorrelator.cc @@ -55,46 +55,47 @@ int main() dataInput->SetNComponents(2); KTFrequencySpectrumFFTW* spectrum0 = new KTFrequencySpectrumFFTW(19, 0, 20); - (*spectrum0)(0)[0] = 0.; (*spectrum0)(0)[1] = 0.; - (*spectrum0)(1)[0] = 0.; (*spectrum0)(1)[1] = 0.; - (*spectrum0)(2)[0] = 1.; (*spectrum0)(2)[1] = 0.; - (*spectrum0)(3)[0] = 2.; (*spectrum0)(3)[1] = 0.; - (*spectrum0)(4)[0] = 3.; (*spectrum0)(4)[1] = 0.; - (*spectrum0)(5)[0] = 4.; (*spectrum0)(5)[1] = 0.; - (*spectrum0)(6)[0] = 5.; (*spectrum0)(6)[1] = 0.; - (*spectrum0)(7)[0] = 6.; (*spectrum0)(7)[1] = 0.; - (*spectrum0)(8)[0] = 0.; (*spectrum0)(8)[1] = 1.; - (*spectrum0)(9)[0] = 0.; (*spectrum0)(9)[1] = 10.; - (*spectrum0)(10)[0] = 0.; (*spectrum0)(10)[1] = 1.; - (*spectrum0)(11)[0] = 6.; (*spectrum0)(11)[1] = 0.; - (*spectrum0)(12)[0] = 5.; (*spectrum0)(12)[1] = 0.; - (*spectrum0)(13)[0] = 4.; (*spectrum0)(13)[1] = 0.; - (*spectrum0)(14)[0] = 3.; (*spectrum0)(14)[1] = 0.; - (*spectrum0)(15)[0] = 2.; (*spectrum0)(15)[1] = 0.; - (*spectrum0)(16)[0] = 1.; (*spectrum0)(16)[1] = 0.; - (*spectrum0)(17)[0] = 0.; (*spectrum0)(17)[1] = 0.; - (*spectrum0)(18)[0] = 0.; (*spectrum0)(18)[1] = 0.; + spectrum0->SetRect(0, 0., 0.); + spectrum0->SetRect(1, 0., 0.); + spectrum0->SetRect(2, 1., 0.); + spectrum0->SetRect(3, 2., 0.); + spectrum0->SetRect(4, 3., 0.); + spectrum0->SetRect(5, 4., 0.); + spectrum0->SetRect(6, 5., 0.); + spectrum0->SetRect(7, 6., 0.); + spectrum0->SetRect(8, 0., 1.); + spectrum0->SetRect(9, 0., 10.); + spectrum0->SetRect(10, 0., 1.); + spectrum0->SetRect(11, 6., 0.); + spectrum0->SetRect(12, 5., 0.); + spectrum0->SetRect(13, 4., 0.); + spectrum0->SetRect(14, 3., 0.); + spectrum0->SetRect(15, 2., 0.); + spectrum0->SetRect(16, 1., 0.); + spectrum0->SetRect(17, 0., 0.); + spectrum0->SetRect(18, 0., 0.); KTFrequencySpectrumFFTW* spectrum1 = new KTFrequencySpectrumFFTW(19, 0, 20); - (*spectrum1)(0)[0] = 0.; (*spectrum1)(0)[1] = 0.; - (*spectrum1)(1)[0] = 0.; (*spectrum1)(1)[1] = 0.; - (*spectrum1)(2)[0] = 0.; (*spectrum1)(2)[1] = 0.; - (*spectrum1)(3)[0] = 5.; (*spectrum1)(3)[1] = 0.; - (*spectrum1)(4)[0] = 5.; (*spectrum1)(4)[1] = 0.; - (*spectrum1)(5)[0] = 5.; (*spectrum1)(5)[1] = 0.; - (*spectrum1)(6)[0] = 5.; (*spectrum1)(6)[1] = 0.; - (*spectrum1)(7)[0] = 0.; (*spectrum1)(7)[1] = 0.; - (*spectrum1)(8)[0] = 2.; (*spectrum1)(8)[1] = 2.; - (*spectrum1)(9)[0] = 2.; (*spectrum1)(9)[1] = 2.; - (*spectrum1)(10)[0] = 2.; (*spectrum1)(10)[1] = 2.; - (*spectrum1)(11)[0] = 0.; (*spectrum1)(11)[1] = 0.; - (*spectrum1)(12)[0] = 5.; (*spectrum1)(12)[1] = 0.; - (*spectrum1)(13)[0] = 5.; (*spectrum1)(13)[1] = 0.; - (*spectrum1)(14)[0] = 5.; (*spectrum1)(14)[1] = 0.; - (*spectrum1)(15)[0] = 5.; (*spectrum1)(15)[1] = 0.; - (*spectrum1)(16)[0] = 0.; (*spectrum1)(16)[1] = 0.; - (*spectrum1)(17)[0] = 0.; (*spectrum1)(17)[1] = 0.; - (*spectrum1)(18)[0] = 0.; (*spectrum1)(18)[1] = 0.; + spectrum1->SetRect(0, 0., 0.); + spectrum1->SetRect(1, 0., 0.); + spectrum1->SetRect(2, 0., 0.); + spectrum1->SetRect(3, 5., 0.); + spectrum1->SetRect(4, 5., 0.); + spectrum1->SetRect(5, 5., 0.); + spectrum1->SetRect(6, 5., 0.); + spectrum1->SetRect(7, 0., 0.); + spectrum1->SetRect(8, 2., 2.); + spectrum1->SetRect(9, 2., 2.); + spectrum1->SetRect(10, 2., 2.); + spectrum1->SetRect(11, 0., 0.); + spectrum1->SetRect(12, 5., 0.); + spectrum1->SetRect(13, 5., 0.); + spectrum1->SetRect(14, 5., 0.); + spectrum1->SetRect(15, 5., 0.); + spectrum1->SetRect(16, 0., 0.); + spectrum1->SetRect(17, 0., 0.); + spectrum1->SetRect(18, 0., 0.); + /**/ dataInput->SetSpectrum(spectrum0, 0); diff --git a/Source/Executables/Validation/TestForwardFFTW.cc b/Source/Executables/Validation/TestForwardFFTW.cc index 21dcb21a2..5fc0d87d2 100644 --- a/Source/Executables/Validation/TestForwardFFTW.cc +++ b/Source/Executables/Validation/TestForwardFFTW.cc @@ -95,7 +95,7 @@ int main() for (unsigned iBin = 0; iBin < nFreqBins; iBin++) { - value = (*frequencySpectrum)(iBin)[0]*(*frequencySpectrum)(iBin)[0] + (*frequencySpectrum)(iBin)[1]*(*frequencySpectrum)(iBin)[1]; + value = frequencySpectrum->GetNorm(iBin); if (value > maxValue) { maxValue = value; @@ -140,7 +140,7 @@ int main() double fsSum = 0.; // units: volts^2 for (unsigned iBin=0; iBinGetNorm(iBin); } KTINFO(vallog, "sum(freqSpectrum[i]^2) = " << fsSum << " V^2"); @@ -169,7 +169,7 @@ int main() for (unsigned iBin = 0; iBin < nFreqBins; iBin++) { - value = (*frequencySpectrum2)(iBin)[0]*(*frequencySpectrum2)(iBin)[0] + (*frequencySpectrum2)(iBin)[1]*(*frequencySpectrum2)(iBin)[1]; + value = frequencySpectrum2->GetNorm(iBin); if (value > maxValue) { maxValue = value; @@ -214,7 +214,7 @@ int main() fsSum = 0.; // units: volts^2 for (unsigned iBin=0; iBinGetNorm(iBin); } KTINFO(vallog, "sum(freqSpectrum2[i]^2) = " << fsSum << " V^2"); diff --git a/Source/Executables/Validation/TestReverseFFTW.cc b/Source/Executables/Validation/TestReverseFFTW.cc index 5f2106704..fc9e98d70 100644 --- a/Source/Executables/Validation/TestReverseFFTW.cc +++ b/Source/Executables/Validation/TestReverseFFTW.cc @@ -51,8 +51,7 @@ int main() // Fill with a sinusoid for (unsigned iBin=0; iBinGetBinCenter(iBin) * mult); - (*fsFFTW)(iBin)[1] = 0.0; + fsFFTW->SetRect( iBin, cos(fsFFTW->GetBinCenter(iBin) * mult), 0.0); } // Create and prepare the FFT @@ -122,7 +121,7 @@ int main() double fsSum = 0.; // units: volts^2 for (unsigned iBin=0; iBinGetNorm(iBin); } KTINFO(vallog, "sum(freqSpectrum[i]^2) = " << fsSum << " V^2"); diff --git a/Source/Executables/Validation/TestSpectrogramStriper.cc b/Source/Executables/Validation/TestSpectrogramStriper.cc index aef19f005..519cd42e9 100644 --- a/Source/Executables/Validation/TestSpectrogramStriper.cc +++ b/Source/Executables/Validation/TestSpectrogramStriper.cc @@ -118,12 +118,12 @@ int main() for (unsigned iSlice = 0; iSlice < slicesPerAcq; ++iSlice) { KTDEBUG(testlog, "New slice"); - (*spectrum)(peakPos)[0] = 10.; + spectrum->SetRect(peakPos, 10., spectrum->GetImag(peakPos)); striper.AddData(header, spectrumData); // update data for the next slice - (*spectrum)(peakPos)[0] = 0.; + spectrum->SetRect(peakPos, 0.0, spectrum->GetImag(peakPos)); peakPos++; if (peakPos == nFreqBins) peakPos = 0; diff --git a/Source/Executables/Validation/TestSpectrumDiscriminator.cc b/Source/Executables/Validation/TestSpectrumDiscriminator.cc index 6f86c16d2..889fc2d8a 100644 --- a/Source/Executables/Validation/TestSpectrumDiscriminator.cc +++ b/Source/Executables/Validation/TestSpectrumDiscriminator.cc @@ -115,13 +115,11 @@ int main() { #ifdef ROOT_FOUND (*spectrum)(iBin).set_polar(rand.Gaus(meanValue, noiseSigma), 0.); - (*spectrumFFTW)(iBin)[0] = rand.Gaus(meanValue, noiseSigma); - (*spectrumFFTW)(iBin)[1] = 0.; + spectrumFFTW->SetRect(iBin, rand.Gaus(meanValue, noiseSigma), 0.); (*powerSpectrum)(iBin) = rand.Gaus(meanValue, noiseSigma); #else (*spectrum)(iBin).set_polar(meanValue, 0.); - (*spectrumFFTW)(iBin)[0] = meanValue; - (*spectrumFFTW)(iBin)[1] = 0.; + spectrumFFTW->SetRect(iBin, meanValue, 0.); (*powerSpectrum)(iBin) = meanValue; #endif } @@ -137,8 +135,9 @@ int main() double multiplier = meanPeakMult; #endif (*spectrum)(iBin).set_polar((*spectrum)(iBin).abs() * multiplier, 0.); - (*spectrumFFTW)(iBin)[0] =(*spectrumFFTW)(iBin)[0] * multiplier; - (*spectrumFFTW)(iBin)[1] =(*spectrumFFTW)(iBin)[1] * multiplier; + spectrumFFTW->SetRect(iBin, + spectrumFFTW->GetReal(iBin)*multiplier, + spectrumFFTW->GetImag(iBin)*multiplier); (*powerSpectrum)(iBin) = (*powerSpectrum)(iBin) * multiplier; //KTINFO(testlog, "Adding peak at bin " << iBin << "; new value: " << (*spectrum)(iBin).abs()); //KTINFO(testlog, "Adding peak at bin " << iBin << "; new value: " << std::sqrt((*spectrumFFTW)(iBin)[0] * (*spectrumFFTW)(iBin)[0] + (*spectrumFFTW)(iBin)[1] * (*spectrumFFTW)(iBin)[1])); diff --git a/Source/Executables/Validation/TestWignerVille.cc b/Source/Executables/Validation/TestWignerVille.cc index 1e5f6d082..f1d3316f6 100644 --- a/Source/Executables/Validation/TestWignerVille.cc +++ b/Source/Executables/Validation/TestWignerVille.cc @@ -137,7 +137,7 @@ int main() KTFrequencySpectrumFFTW* spectrum = spectra(iX); for (unsigned iY=0; iYGetAbs(iY); histOut->SetBinContent(iX+1, iY+1, value); } } diff --git a/Source/SpectrumAnalysis/KTAmplitudeDistributor.cc b/Source/SpectrumAnalysis/KTAmplitudeDistributor.cc index 9d310ab18..0a299e99f 100644 --- a/Source/SpectrumAnalysis/KTAmplitudeDistributor.cc +++ b/Source/SpectrumAnalysis/KTAmplitudeDistributor.cc @@ -331,7 +331,7 @@ namespace Katydid for (unsigned iBin = fMinBin; iBin <= fMaxBin; iBin++) { - fBuffer[fNSlicesProcessed][component][iBin] = sqrt((*spectrum)(iBin)[0]*(*spectrum)(iBin)[0] + (*spectrum)(iBin)[1]*(*spectrum)(iBin)[1]); + fBuffer[fNSlicesProcessed][component][iBin] = spectrum->GetAbs(iBin); } fNBuffered++; KTDEBUG(adlog, "Buffer now contains " << fNBuffered << " distributions; " << fNSlicesProcessed + 1 << " slices processed"); @@ -342,7 +342,7 @@ namespace Katydid { for (unsigned iBin = fMinBin; iBin <= fMaxBin; iBin++) { - fDistributions->AddToDist(iBin, sqrt((*spectrum)(iBin)[0]*(*spectrum)(iBin)[0] + (*spectrum)(iBin)[1]*(*spectrum)(iBin)[1]), component); + fDistributions->AddToDist(iBin, spectrum->GetAbs(iBin), component); } return true; } diff --git a/Source/SpectrumAnalysis/KTAxialAggregator.cc b/Source/SpectrumAnalysis/KTAxialAggregator.cc index 705a3aa21..9567cd0c4 100644 --- a/Source/SpectrumAnalysis/KTAxialAggregator.cc +++ b/Source/SpectrumAnalysis/KTAxialAggregator.cc @@ -64,8 +64,7 @@ namespace Katydid double imagVal = freqSpectrum->GetImag(iFreqBin); double summedRealVal = realVal + newFreqSpectrum->GetReal(iFreqBin); double summedImagVal = imagVal + newFreqSpectrum->GetImag(iFreqBin); - (*newFreqSpectrum)(iFreqBin)[0] = summedRealVal; - (*newFreqSpectrum)(iFreqBin)[1] = summedImagVal; + newFreqSpectrum->SetRect(iFreqBin, summedRealVal, summedImagVal); // This doesn't seem to work //(*newFreqSpectrum)+=*fftwData.GetSpectrumFFTW(iComponent+nComponents*iRing); } diff --git a/Source/SpectrumAnalysis/KTChannelAggregator.cc b/Source/SpectrumAnalysis/KTChannelAggregator.cc index 52e90aa2e..943cb3597 100644 --- a/Source/SpectrumAnalysis/KTChannelAggregator.cc +++ b/Source/SpectrumAnalysis/KTChannelAggregator.cc @@ -274,8 +274,7 @@ namespace Katydid ApplyPhaseShift(realVal, imagVal, phaseShift); double summedRealVal = realVal + newFreqSpectrum->GetReal(iFreqBin); double summedImagVal = imagVal + newFreqSpectrum->GetImag(iFreqBin); - (*newFreqSpectrum)(iFreqBin)[0] = summedRealVal; - (*newFreqSpectrum)(iFreqBin)[1] = summedImagVal; + (*newFreqSpectrum)(iFreqBin) = std::complex{ summedRealVal, summedImagVal }; } // End of loop over freq bins } // End of loop over all comps newFreqSpectrum->SetNTimeBins(nTimeBins); diff --git a/Source/SpectrumAnalysis/KTConvolution.cc b/Source/SpectrumAnalysis/KTConvolution.cc index 39e77881b..e46b0baf5 100644 --- a/Source/SpectrumAnalysis/KTConvolution.cc +++ b/Source/SpectrumAnalysis/KTConvolution.cc @@ -439,16 +439,14 @@ namespace Katydid void KTConvolution1D::ConjugateAndReverse( KTFrequencySpectrumFFTW& spectrum ) { - fftw_complex temp; int nBins = spectrum.GetNFrequencyBins(); for( int iBin = 0; iBin < (nBins+1)/2; ++iBin ) { - temp[0] = spectrum(iBin)[0]; - temp[1] = spectrum(iBin)[1]; - spectrum(iBin)[0] = spectrum(nBins - iBin - 1)[0]; - spectrum(iBin)[1] = -1. * spectrum(nBins - iBin - 1)[1]; - spectrum(nBins - iBin - 1)[0] = temp[0]; - spectrum(nBins - iBin - 1)[1] = -1. * temp[1]; + auto tmpReal = spectrum.GetReal(iBin); + auto tmpImag = spectrum.GetImag(iBin); + spectrum.SetRect(iBin, spectrum.GetReal(nBins - iBin - 1), + -1. * spectrum.GetImag(nBins - iBin - 1)); + spectrum.SetRect(nBins - iBin -1, tmpReal, tmpImag); } return; @@ -487,8 +485,8 @@ namespace Katydid } else { - fInputArrayComplex[nBin][0] = (*initialSpectrum)(position)[0]; - fInputArrayComplex[nBin][1] = (*initialSpectrum)(position)[1]; + fInputArrayComplex[nBin][0] = initialSpectrum->GetReal(position); + fInputArrayComplex[nBin][1] = initialSpectrum->GetImag(position); } return; @@ -519,8 +517,10 @@ namespace Katydid void KTConvolution1D::SetOutputArray( int position, int nBin, KTFrequencySpectrumFFTW& transformedFSFFTW, double norm ) { - transformedFSFFTW(position)[0] = fOutputArrayComplex[nBin][0] / (double)norm; - transformedFSFFTW(position)[1] = fOutputArrayComplex[nBin][1] / (double)norm; + + transformedFSFFTW.SetRect(position, + fOutputArrayComplex[nBin][0] / (double)norm, + fOutputArrayComplex[nBin][1] / (double)norm); return; } diff --git a/Source/SpectrumAnalysis/KTDataAccumulator.cc b/Source/SpectrumAnalysis/KTDataAccumulator.cc index 18c480f8c..5f0ec951f 100644 --- a/Source/SpectrumAnalysis/KTDataAccumulator.cc +++ b/Source/SpectrumAnalysis/KTDataAccumulator.cc @@ -440,11 +440,15 @@ namespace Katydid KTFrequencySpectrumVariance* varSpect = devData.GetSpectrum(iComponent); for (unsigned iBin = 0; iBin < arraySize; ++iBin) { - (*avSpect)(iBin)[0] = (*avSpect)(iBin)[0] * remainingFrac + (*newSpect)(iBin)[0] * fAveragingFrac; - (*avSpect)(iBin)[1] = (*avSpect)(iBin)[1] * remainingFrac + (*newSpect)(iBin)[1] * fAveragingFrac; - - (*varSpect)(iBin) = (*varSpect)(iBin) * remainingFrac + ((*newSpect)(iBin)[0] * (*newSpect)(iBin)[0] + (*newSpect)(iBin)[1] * (*newSpect)(iBin)[1]) * fAveragingFrac; + // (*avSpect)(iBin)[0] = (*avSpect)(iBin)[0] * remainingFrac + (*newSpect)(iBin)[0] * fAveragingFrac; + // (*avSpect)(iBin)[1] = (*avSpect)(iBin)[1] * remainingFrac + (*newSpect)(iBin)[1] * fAveragingFrac; + + (*varSpect)(iBin) = (*varSpect)(iBin) * remainingFrac + + newSpect->GetNorm(iBin) * fAveragingFrac; } + + (*avSpect) = avSpect->Scale(remainingFrac) + newSpect->Scale(fAveragingFrac); + } return true; @@ -587,7 +591,7 @@ namespace Katydid unsigned nBins = varSpect->GetNFrequencyBins(); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - (*varSpect)(iBin) = (*varSpect)(iBin) - (*avSpect)(iBin)[0] * (*avSpect)(iBin)[0] - (*avSpect)(iBin)[1] * (*avSpect)(iBin)[1]; + (*varSpect)(iBin) = (*varSpect)(iBin) - avSpect->GetNorm(iBin); } } return true; @@ -650,7 +654,7 @@ namespace Katydid unsigned nBins = varSpect->GetNFrequencyBins(); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - (*varSpect)(iBin) = (*varSpect)(iBin) - (*avSpect)(iBin)[0] * (*avSpect)(iBin)[0] - (*avSpect)(iBin)[1] * (*avSpect)(iBin)[1]; + (*varSpect)(iBin) = (*varSpect)(iBin) - avSpect->GetNorm(iBin); } } return true; diff --git a/Source/SpectrumAnalysis/KTGainNormalization.cc b/Source/SpectrumAnalysis/KTGainNormalization.cc index 6587b3de2..04fe2836b 100644 --- a/Source/SpectrumAnalysis/KTGainNormalization.cc +++ b/Source/SpectrumAnalysis/KTGainNormalization.cc @@ -284,14 +284,16 @@ namespace Katydid #pragma omp for private(iBin) for (iBin=0; iBin < fMinBin; ++iBin) { - (*newSpectrum)(iBin)[0] = (*frequencySpectrum)(iBin)[0]; - (*newSpectrum)(iBin)[1] = (*frequencySpectrum)(iBin)[1]; + newSpectrum->SetRect(iBin, + frequencySpectrum->GetReal(iBin), + frequencySpectrum->GetImag(iBin)); } #pragma omp for private(iBin) for (iBin=fMaxBin+1; iBin < nSpectrumBins; ++iBin) { - (*newSpectrum)(iBin)[0] = (*frequencySpectrum)(iBin)[0]; - (*newSpectrum)(iBin)[1] = (*frequencySpectrum)(iBin)[1]; + newSpectrum->SetRect(iBin, + frequencySpectrum->GetReal(iBin), + frequencySpectrum->GetImag(iBin)); } // Then scale the bins within the scaling range @@ -299,10 +301,9 @@ namespace Katydid #pragma omp for private(iBin) for (iBin=fMinBin; iBin < fMaxBin+1; ++iBin) { - value.set_rect((*frequencySpectrum)(iBin)[0], (*frequencySpectrum)(iBin)[1]); + value.set_rect(frequencySpectrum->GetReal(iBin), frequencySpectrum->GetImag(iBin)); value.set_polar(normalizedMean + (value.abs() - (*splineImp)(iBin - fMinBin)) * sqrt(normalizedVariance / (*varSplineImp)(iBin - fMinBin)), value.arg()); - (*newSpectrum)(iBin)[0] = real(value); - (*newSpectrum)(iBin)[1] = imag(value); + newSpectrum->SetRect(iBin, real(value), imag(value)); } } diff --git a/Source/SpectrumAnalysis/KTGainVariationProcessor.cc b/Source/SpectrumAnalysis/KTGainVariationProcessor.cc index dbc267ddf..759195ce1 100644 --- a/Source/SpectrumAnalysis/KTGainVariationProcessor.cc +++ b/Source/SpectrumAnalysis/KTGainVariationProcessor.cc @@ -462,7 +462,8 @@ namespace Katydid double mean = 0.; for (unsigned iBin=fitPointStartBin; iBinGetAbs(iBin); } mean /= (double)nBinsPerFitPoint; yVals[iFitPoint] = mean; diff --git a/Source/SpectrumAnalysis/KTSpectrogramStriper.cc b/Source/SpectrumAnalysis/KTSpectrogramStriper.cc index 5dce598d9..2e36e96fb 100644 --- a/Source/SpectrumAnalysis/KTSpectrogramStriper.cc +++ b/Source/SpectrumAnalysis/KTSpectrogramStriper.cc @@ -110,12 +110,12 @@ namespace Katydid return true; } + // Why this function? The copy constructor of KTFrequencySpectrumFFTW should be responsible for this void KTSpectrogramStriper::CopySpectrum(const KTFrequencySpectrumFFTW* source, KTFrequencySpectrumFFTW* dest, unsigned arraySize) { for (unsigned iBin = 0; iBin < arraySize; ++iBin) { - (*dest)(iBin)[0] = (*source)(iBin)[0]; - (*dest)(iBin)[1] = (*source)(iBin)[1]; + dest->SetRect(iBin, source->GetReal(iBin), source->GetImag(iBin)); } } diff --git a/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc b/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc index 3f80dc461..90e4300ed 100644 --- a/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc +++ b/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc @@ -223,7 +223,7 @@ namespace Katydid #pragma omp parallel for reduction(+:mean) for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { - magnitude[iBin] = sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); + magnitude[iBin] = spectrum->GetAbs(iBin); //sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); mean += magnitude[iBin]; variance += magnitude[iBin] * magnitude[iBin]; } @@ -472,7 +472,7 @@ namespace Katydid neighborhoodAmplitude = 0; for (unsigned jBin = iBin-fNeighborhoodRadius; jBin<= iBin+fNeighborhoodRadius; ++jBin) { - neighborhoodAmplitude += sqrt((*spectrum)(jBin)[0] * (*spectrum)(jBin)[0] + (*spectrum)(jBin)[1] * (*spectrum)(jBin)[1]); + neighborhoodAmplitude += spectrum->GetAbs(jBin); //sqrt((*spectrum)(jBin)[0] * (*spectrum)(jBin)[0] + (*spectrum)(jBin)[1] * (*spectrum)(jBin)[1]); } } void KTSpectrumDiscriminator::SumAdjacentBinAmplitude(const KTFrequencySpectrumPolar* spectrum, double& neighborhoodAmplitude, const unsigned& iBin) diff --git a/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc b/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc index 8e54838f3..c60a15821 100644 --- a/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc +++ b/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc @@ -536,7 +536,7 @@ namespace Katydid #pragma omp parallel for private(value) for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { - double value = sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); + double value = spectrum->GetAbs(iBin); //sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); double threshold = thresholdMult * (*splineImp)(iBin - fMinBin); double mean = (*splineImp)(iBin - fMinBin); double variance = (*varSplineImp)(iBin - fMinBin); @@ -572,7 +572,7 @@ namespace Katydid double mean = (*splineImp)(iBin - fMinBin); double variance = (*varSplineImp)(iBin - fMinBin); double threshold = mean + fSigmaThreshold * sqrt( variance ); - double value = sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); + double value = spectrum->GetAbs(iBin); // sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); if (value >= threshold) { double neighborhoodAmplitude = 0.; @@ -714,7 +714,7 @@ namespace Katydid neighborhoodAmplitude = 0; for (unsigned jBin = iBin-fNeighborhoodRadius; jBin <= iBin+fNeighborhoodRadius; ++jBin) { - neighborhoodAmplitude += sqrt((*spectrum)(jBin)[0] * (*spectrum)(jBin)[0] + (*spectrum)(jBin)[1] * (*spectrum)(jBin)[1]); + neighborhoodAmplitude += spectrum->GetAbs(jBin); //sqrt((*spectrum)(jBin)[0] * (*spectrum)(jBin)[0] + (*spectrum)(jBin)[1] * (*spectrum)(jBin)[1]); } } void KTVariableSpectrumDiscriminator::SumAdjacentBinAmplitude(const KTFrequencySpectrumPolar* spectrum, double& neighborhoodAmplitude, const unsigned& iBin) From 568cfb5b68e2a7c6f7998239be523246ad1650ff Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 20:37:37 +0200 Subject: [PATCH 10/48] Fix in KTFrequencySpectrumFFTW --- Source/Data/Transform/KTFrequencySpectrumFFTW.cc | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc index 343acbb66..eeed5e985 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc @@ -61,10 +61,15 @@ namespace Katydid //KTINFO(fslog, "neg freq offset: " << fLeftOfCenterOffset); } + // I have doubts about this function KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(std::initializer_list value, size_t nBins, double rangeMin, double rangeMax, bool arrayOrderIsFlipped) : KTFrequencySpectrumFFTW(nBins, rangeMin, rangeMax, arrayOrderIsFlipped) { - std::copy(value.begin(), value.end(), this->begin()); + for (unsigned index = 0; index < nBins; ++index) + { + std::copy(value.begin(), value.end(), fData[index]); + } + //std::copy(value.begin(), value.end(), this->begin()); } KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig) : From 2bd273b7b1bc38ad2a16453f705c612e5f8f0456 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 20:46:17 +0200 Subject: [PATCH 11/48] Modify KTFrequenySpectrumFFTW Removes the copy constructor to replace it with the default one --- .../Data/Transform/KTFrequencySpectrumFFTW.cc | 37 ++++++++++--------- .../Data/Transform/KTFrequencySpectrumFFTW.hh | 5 ++- 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc index eeed5e985..657f9fa9b 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc @@ -72,25 +72,26 @@ namespace Katydid //std::copy(value.begin(), value.end(), this->begin()); } - KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig) : - KTPhysicalArray< 1, std::complex >(orig), - KTFrequencySpectrum(), - fIsArrayOrderFlipped(orig.fIsArrayOrderFlipped), - fIsSizeEven(orig.fIsSizeEven), - fLeftOfCenterOffset(orig.fLeftOfCenterOffset), - fCenterBin(orig.fCenterBin), - fConstBinAccess(orig.fConstBinAccess), - fBinAccess(orig.fBinAccess), - fNTimeBins(orig.fNTimeBins), - fPointCache() - { - } - - KTFrequencySpectrumFFTW::~KTFrequencySpectrumFFTW() - { - } + //copy constructor, destructor and copy assignment operator + //shouldn't be necessary. I think the default ones will do the right thing + + //~ KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig) : + //~ KTPhysicalArray< 1, std::complex >(orig), + //~ KTFrequencySpectrum(), + //~ fIsArrayOrderFlipped(orig.fIsArrayOrderFlipped), + //~ fIsSizeEven(orig.fIsSizeEven), + //~ fLeftOfCenterOffset(orig.fLeftOfCenterOffset), + //~ fCenterBin(orig.fCenterBin), + //~ fConstBinAccess(orig.fConstBinAccess), + //~ fBinAccess(orig.fBinAccess), + //~ fNTimeBins(orig.fNTimeBins), + //~ fPointCache() + //~ { + //~ } - //Is the copy assignment operator necessary? I think the default one will do the right thing + //~ KTFrequencySpectrumFFTW::~KTFrequencySpectrumFFTW() + //~ { + //~ } //~ KTFrequencySpectrumFFTW& KTFrequencySpectrumFFTW::operator=(const KTFrequencySpectrumFFTW& rhs) //~ { diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh index 668e62e65..5c977019c 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh @@ -26,8 +26,9 @@ namespace Katydid KTFrequencySpectrumFFTW(); KTFrequencySpectrumFFTW(size_t nBins, double rangeMin=0., double rangeMax=1., bool arrayOrderIsFlipped=false); KTFrequencySpectrumFFTW(std::initializer_list value, size_t nBins, double rangeMin=0., double rangeMax=1., bool arrayOrderIsFlipped=false); - KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig); - virtual ~KTFrequencySpectrumFFTW(); + //see comment in cc file + //KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig); + virtual ~KTFrequencySpectrumFFTW() = default; public: bool GetIsArrayOrderFlipped() const; From 7c4e3f76c221d972d976517e732aef6b8413c5ea Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 20:50:48 +0200 Subject: [PATCH 12/48] Replace fft_complex in KTTimeSeriesFFTW --- Source/Data/Time/KTTimeSeriesFFTW.cc | 9 +++++---- Source/Data/Time/KTTimeSeriesFFTW.hh | 11 +++++------ 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Source/Data/Time/KTTimeSeriesFFTW.cc b/Source/Data/Time/KTTimeSeriesFFTW.cc index 03b59197f..658e07e4e 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.cc +++ b/Source/Data/Time/KTTimeSeriesFFTW.cc @@ -23,13 +23,13 @@ namespace Katydid KTTimeSeriesFFTW::KTTimeSeriesFFTW() : KTTimeSeries(), - KTPhysicalArray< 1, fftw_complex >() + KTPhysicalArray< 1, std::complex >() { } KTTimeSeriesFFTW::KTTimeSeriesFFTW(size_t nBins, double rangeMin, double rangeMax) : KTTimeSeries(), - KTPhysicalArray< 1, fftw_complex >(nBins, rangeMin, rangeMax) + KTPhysicalArray< 1, std::complex >(nBins, rangeMin, rangeMax) { } @@ -44,11 +44,12 @@ namespace Katydid { std::copy(value.begin(), value.end(), fData[iBin]); } + } KTTimeSeriesFFTW::KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig) : KTTimeSeries(), - KTPhysicalArray< 1, fftw_complex >(orig) + KTPhysicalArray< 1, std::complex >(orig) { } @@ -58,7 +59,7 @@ namespace Katydid KTTimeSeriesFFTW& KTTimeSeriesFFTW::operator=(const KTTimeSeriesFFTW& rhs) { - KTPhysicalArray< 1, fftw_complex >::operator=(rhs); + KTPhysicalArray< 1, std::complex >::operator=(rhs); return *this; } diff --git a/Source/Data/Time/KTTimeSeriesFFTW.hh b/Source/Data/Time/KTTimeSeriesFFTW.hh index 20e05cc8e..2ac07b785 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.hh +++ b/Source/Data/Time/KTTimeSeriesFFTW.hh @@ -8,7 +8,7 @@ #ifndef KTTIMESERIESFFTW_HH_ #define KTTIMESERIESFFTW_HH_ -#include "KTPhysicalArrayFFTW.hh" +#include "KTPhysicalArrayComplex.hh" #include "KTTimeSeries.hh" #include @@ -18,7 +18,7 @@ namespace Katydid - class KTTimeSeriesFFTW : public KTTimeSeries, public KTPhysicalArray< 1, fftw_complex > + class KTTimeSeriesFFTW : public KTTimeSeries, public KTPhysicalArray< 1, std::complex > { public: KTTimeSeriesFFTW(); @@ -49,7 +49,7 @@ namespace Katydid inline void KTTimeSeriesFFTW::Scale(double scale) { - this->KTPhysicalArray< 1, fftw_complex >::operator*=(scale); + this->KTPhysicalArray< 1, std::complex >::operator*=(scale); return; } @@ -65,14 +65,13 @@ namespace Katydid inline void KTTimeSeriesFFTW::SetValue(unsigned bin, double value) { - (*this)(bin)[0] = value; - (*this)(bin)[1] = 0.; + (*this)(bin) = std::complex {value, 0.}; return; } inline double KTTimeSeriesFFTW::GetValue(unsigned bin) const { - return (*this)(bin)[0]; + return (*this)(bin).real(); } } /* namespace Katydid */ From 53799f98bdf90f9d4083c198150495d3445ee16b Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 20:52:14 +0200 Subject: [PATCH 13/48] Modify KTTimeSeriesFFTW Applies the rule of zero --- Source/Data/Time/KTTimeSeriesFFTW.cc | 34 +++++++++++++++------------- Source/Data/Time/KTTimeSeriesFFTW.hh | 8 ++++--- 2 files changed, 23 insertions(+), 19 deletions(-) diff --git a/Source/Data/Time/KTTimeSeriesFFTW.cc b/Source/Data/Time/KTTimeSeriesFFTW.cc index 658e07e4e..3df56802e 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.cc +++ b/Source/Data/Time/KTTimeSeriesFFTW.cc @@ -46,22 +46,24 @@ namespace Katydid } } - - KTTimeSeriesFFTW::KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig) : - KTTimeSeries(), - KTPhysicalArray< 1, std::complex >(orig) - { - } - - KTTimeSeriesFFTW::~KTTimeSeriesFFTW() - { - } - - KTTimeSeriesFFTW& KTTimeSeriesFFTW::operator=(const KTTimeSeriesFFTW& rhs) - { - KTPhysicalArray< 1, std::complex >::operator=(rhs); - return *this; - } + + // rule of zero, the default ones should do the job + + //~ KTTimeSeriesFFTW::KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig) : + //~ KTTimeSeries(), + //~ KTPhysicalArray< 1, std::complex >(orig) + //~ { + //~ } + + //~ KTTimeSeriesFFTW::~KTTimeSeriesFFTW() + //~ { + //~ } + + //~ KTTimeSeriesFFTW& KTTimeSeriesFFTW::operator=(const KTTimeSeriesFFTW& rhs) + //~ { + //~ KTPhysicalArray< 1, std::complex >::operator=(rhs); + //~ return *this; + //~ } void KTTimeSeriesFFTW::Print(unsigned startPrint, unsigned nToPrint) const { diff --git a/Source/Data/Time/KTTimeSeriesFFTW.hh b/Source/Data/Time/KTTimeSeriesFFTW.hh index 2ac07b785..e5d839799 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.hh +++ b/Source/Data/Time/KTTimeSeriesFFTW.hh @@ -24,10 +24,12 @@ namespace Katydid KTTimeSeriesFFTW(); KTTimeSeriesFFTW(size_t nBins, double rangeMin=0., double rangeMax=1.); KTTimeSeriesFFTW(std::initializer_list value, size_t nBins, double rangeMin=0., double rangeMax=1.); - KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig); - virtual ~KTTimeSeriesFFTW(); + + //rule of zero + //KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig); + virtual ~KTTimeSeriesFFTW() = default; - KTTimeSeriesFFTW& operator=(const KTTimeSeriesFFTW& rhs); + //KTTimeSeriesFFTW& operator=(const KTTimeSeriesFFTW& rhs); virtual void Scale(double scale); From bf1f31c949cc5bbb0ac84489214bdeaad04b072b Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 21:41:21 +0200 Subject: [PATCH 14/48] Fix error in KTTimeSeriesFFTW and KTFrequencySpectrumFFTW Missing address operator --- Source/Data/Time/KTTimeSeriesFFTW.cc | 2 +- Source/Data/Transform/KTFrequencySpectrumFFTW.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Data/Time/KTTimeSeriesFFTW.cc b/Source/Data/Time/KTTimeSeriesFFTW.cc index 3df56802e..6db0aff14 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.cc +++ b/Source/Data/Time/KTTimeSeriesFFTW.cc @@ -42,7 +42,7 @@ namespace Katydid } for (unsigned iBin = 0; iBin < nBins; ++iBin) { - std::copy(value.begin(), value.end(), fData[iBin]); + std::copy(value.begin(), value.end(), &fData[iBin]); } } diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc index 657f9fa9b..d1bb95799 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc @@ -67,7 +67,7 @@ namespace Katydid { for (unsigned index = 0; index < nBins; ++index) { - std::copy(value.begin(), value.end(), fData[index]); + std::copy(value.begin(), value.end(), &fData[index]); } //std::copy(value.begin(), value.end(), this->begin()); } From 187cb5eaad2718db1ec88d15197c8b19b3c4bb1c Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 20 Apr 2021 21:44:55 +0200 Subject: [PATCH 15/48] Fix some remaining fftw_complex artifacts in KTTimeSeriesFFTW --- Source/Data/Time/KTTimeSeriesFFTW.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Data/Time/KTTimeSeriesFFTW.cc b/Source/Data/Time/KTTimeSeriesFFTW.cc index 6db0aff14..66d938ae2 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.cc +++ b/Source/Data/Time/KTTimeSeriesFFTW.cc @@ -84,7 +84,7 @@ namespace Katydid TH1D* hist = new TH1D(name.c_str(), "Time Series", (int)nBins, GetRangeMin(), GetRangeMax()); for (unsigned iBin=0; iBinSetBinContent((int)iBin+1, (*this)(iBin)[0]); + hist->SetBinContent((int)iBin+1, (*this)(iBin).real()); } hist->SetXTitle("Time (s)"); hist->SetYTitle("Voltage (V)"); @@ -99,7 +99,7 @@ namespace Katydid double value; for (unsigned iBin=0; iBin tMaxMag) tMaxMag = value; @@ -110,7 +110,7 @@ namespace Katydid { //value = (*this)(iBin)[0]; //hist->Fill(value*value); - hist->Fill((*this)(iBin)[0]); + hist->Fill((*this)(iBin).real()); } hist->SetXTitle("Voltage (V)"); return hist; From 2436623fdf69804329632be01fea7844596856bf Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 21 Apr 2021 15:02:18 +0200 Subject: [PATCH 16/48] Add GetNorm to the virtual KTFrequencySpectrum --- Source/Data/Transform/KTFrequencySpectrum.hh | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Data/Transform/KTFrequencySpectrum.hh b/Source/Data/Transform/KTFrequencySpectrum.hh index a2bd020b4..d3c8320d6 100644 --- a/Source/Data/Transform/KTFrequencySpectrum.hh +++ b/Source/Data/Transform/KTFrequencySpectrum.hh @@ -32,6 +32,7 @@ namespace Katydid virtual double GetAbs(unsigned bin) const = 0; virtual double GetArg(unsigned bin) const = 0; + virtual double GetNorm(unsigned bin) const = 0; virtual void SetPolar(unsigned bin, double abs, double arg) = 0; From f5277b71a869390016d2d29455610f727eeed6af Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 21 Apr 2021 15:32:48 +0200 Subject: [PATCH 17/48] Add GetNorm to the KTFrequencySpectrumPolar --- Source/Data/Transform/KTFrequencySpectrumPolar.hh | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Source/Data/Transform/KTFrequencySpectrumPolar.hh b/Source/Data/Transform/KTFrequencySpectrumPolar.hh index 9f0113827..476f99c67 100644 --- a/Source/Data/Transform/KTFrequencySpectrumPolar.hh +++ b/Source/Data/Transform/KTFrequencySpectrumPolar.hh @@ -41,6 +41,7 @@ namespace Katydid virtual double GetAbs(unsigned bin) const; virtual double GetArg(unsigned bin) const; + virtual double GetNorm(unsigned bin) const; virtual void SetPolar(unsigned bin, double abs, double arg); @@ -93,6 +94,12 @@ namespace Katydid { return (*this)(bin).arg(); } + + inline double KTFrequencySpectrumPolar::GetNorm(unsigned bin) const + { + double abs = (*this)(bin).abs(); + return abs*abs; + } inline void KTFrequencySpectrumPolar::SetPolar(unsigned bin, double abs, double arg) { From 63799be59dd816dce3367b9a9aa0cf0773aedbe9 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 21 Apr 2021 15:35:52 +0200 Subject: [PATCH 18/48] Add functions to KTTimeSeriesFFTW --- Source/Data/Time/KTTimeSeriesFFTW.hh | 49 ++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/Source/Data/Time/KTTimeSeriesFFTW.hh b/Source/Data/Time/KTTimeSeriesFFTW.hh index e5d839799..61eda4b5b 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.hh +++ b/Source/Data/Time/KTTimeSeriesFFTW.hh @@ -38,6 +38,17 @@ namespace Katydid virtual void SetValue(unsigned bin, double value); virtual double GetValue(unsigned bin) const; + + double GetReal(unsigned bin) const; + double GetImag(unsigned bin) const; + + void SetRect(unsigned bin, double real, double imag); + + double GetAbs(unsigned bin) const; + double GetArg(unsigned bin) const; + double GetNorm(unsigned bin) const; + + void SetPolar(unsigned bin, double abs, double arg); virtual void Print(unsigned startPrint, unsigned nToPrint) const; @@ -75,6 +86,44 @@ namespace Katydid { return (*this)(bin).real(); } + + inline double KTTimeSeriesFFTW::GetReal(unsigned bin) const + { + return (*this)(bin).real(); + } + + inline double KTTimeSeriesFFTW::GetImag(unsigned bin) const + { + return (*this)(bin).imag(); + } + + inline void KTTimeSeriesFFTW::SetRect(unsigned bin, double real, double imag) + { + + (*this)(bin) = std::complex(real, imag); + return; + } + + inline double KTTimeSeriesFFTW::GetAbs(unsigned bin) const + { + return std::abs((*this)(bin)); + } + + inline double KTTimeSeriesFFTW::GetNorm(unsigned bin) const + { + return std::norm((*this)(bin)); + } + + inline double KTTimeSeriesFFTW::GetArg(unsigned bin) const + { + return std::arg((*this)(bin)); + } + + inline void KTTimeSeriesFFTW::SetPolar(unsigned bin, double abs, double arg) + { + (*this)(bin) = std::polar(abs, arg); + return; + } } /* namespace Katydid */ #endif /* KTTIMESERIESFFTW_HH_ */ From f41e79c57468a5cf51c54cb68a77013a482be177 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 21 Apr 2021 20:18:20 +0200 Subject: [PATCH 19/48] Modify FFT processors They now also take std::complex for the timeseries objects --- Source/Transform/KTForwardFFTW.cc | 8 +- Source/Transform/KTFractionalFFT.cc | 116 +++++++++++++++++----------- Source/Transform/KTReverseFFTW.cc | 7 +- 3 files changed, 82 insertions(+), 49 deletions(-) diff --git a/Source/Transform/KTForwardFFTW.cc b/Source/Transform/KTForwardFFTW.cc index a18cdcf7a..a7c2fd644 100644 --- a/Source/Transform/KTForwardFFTW.cc +++ b/Source/Transform/KTForwardFFTW.cc @@ -542,7 +542,13 @@ namespace Katydid void KTForwardFFTW::DoTransform(const KTTimeSeriesFFTW* tsIn, KTFrequencySpectrumFFTW* fsOut) const { - fftw_execute_dft(fForwardPlan, tsIn->GetData(), reinterpret_cast(fsOut->GetData().data())); + fftw_complex *dataIn = + const_cast( + reinterpret_cast( + tsIn->GetData().data())); + fftw_complex *dataOut = reinterpret_cast( + fsOut->GetData().data()); + fftw_execute_dft(fForwardPlan, dataIn, dataOut); (*fsOut) *= sqrt(1. / (double) fTimeSize); return; } diff --git a/Source/Transform/KTFractionalFFT.cc b/Source/Transform/KTFractionalFFT.cc index 4d930ab4f..6daeb25b1 100644 --- a/Source/Transform/KTFractionalFFT.cc +++ b/Source/Transform/KTFractionalFFT.cc @@ -98,9 +98,7 @@ namespace Katydid KTDEBUG(evlog, "Processing component: " << iComponent); // Initialize vars - KTTimeSeriesFFTW* tsDup = new KTTimeSeriesFFTW( slice.GetSliceSize(), 0.0, slice.GetSliceLength() ); - double norm = 0.; // norm of the current TS value - double phase = 0.; // argument of current TS value + //KTTimeSeriesFFTW* tsDup = new KTTimeSeriesFFTW( slice.GetSliceSize(), 0.0, slice.GetSliceLength() ); // get TS from data object and make the new TS; remains owned by tsData KTTimeSeriesFFTW* ts = dynamic_cast< KTTimeSeriesFFTW* >(tsData.GetTimeSeries( iComponent )); @@ -110,27 +108,37 @@ namespace Katydid continue; } - // Loop through all bins - for( unsigned iBin = 0; iBin < ts->GetNTimeBins(); ++iBin ) - { - // Obtain norm and phase from components - norm = std::hypot( (*ts)(iBin)[0], (*ts)(iBin)[1] ); - phase = atan2( (*ts)(iBin)[1], (*ts)(iBin)[0] ); + //~ // Loop through all bins + //~ for( unsigned iBin = 0; iBin < ts->GetNTimeBins(); ++iBin ) + //~ { + //~ // Obtain norm and phase from components + //~ norm = std::hypot( (*ts)(iBin)[0], (*ts)(iBin)[1] ); + //~ phase = atan2( (*ts)(iBin)[1], (*ts)(iBin)[0] ); - // Shift phase - phase -= chirpRate * ((double)iBin + binOffset) * ((double)iBin + binOffset); + //~ // Shift phase + //~ phase -= chirpRate * ((double)iBin + binOffset) * ((double)iBin + binOffset); - // Assign components from norm and new phase - (*tsDup)(iBin)[0] = norm * cos( phase ); - (*tsDup)(iBin)[1] = norm * sin( phase ); - } + //~ // Assign components from norm and new phase + //~ (*tsDup)(iBin)[0] = norm * cos( phase ); + //~ (*tsDup)(iBin)[1] = norm * sin( phase ); + //~ } + + KTTimeSeriesFFTW* tsDup = new KTTimeSeriesFFTW( *ts ); //copy ts + + // Shift phase + auto phase = -chirpRate * square( + Eigen::ArrayXd::LinSpaced(ts->GetNTimeBins(), 0, + ts->GetNTimeBins()-1) + binOffset); + tsDup->GetData() *= exp(std::complex{0.0, 1.0}*phase); + newTSData.SetTimeSeries( tsDup, iComponent ); // tsDup now owned by newTSData } return true; } + // This function should be split in smaller pieces! bool KTFractionalFFT::ProcessTimeSeries( KTTimeSeriesData& tsData, KTTimeSeriesData& newTSData, KTFrequencySpectrumDataFFTW& newFSData, KTSliceHeader& slice ) { KTDEBUG(evlog, "Receiving time series for fractional FFT"); @@ -150,11 +158,7 @@ namespace Katydid for( unsigned iComponent = 0; iComponent < tsData.GetNComponents(); ++iComponent ) { - KTDEBUG(evlog, "Processing component: " << iComponent); - - // Initialize vars - double norm = 0.; - double phase = 0.; + KTDEBUG(evlog, "Processing component: " << iComponent); // get TS from data object and make the new TS; remains owned by tsData KTTimeSeriesFFTW* ts = dynamic_cast< KTTimeSeriesFFTW* >(tsData.GetTimeSeries( iComponent )); @@ -165,19 +169,29 @@ namespace Katydid } // First chirp transform - for( unsigned iBin = 0; iBin < ts->GetNTimeBins(); ++iBin ) - { - // Obtain norm and phase from components - norm = std::hypot( (*ts)(iBin)[0], (*ts)(iBin)[1] ); - phase = atan2( (*ts)(iBin)[1], (*ts)(iBin)[0] ); + //~ for( unsigned iBin = 0; iBin < ts->GetNTimeBins(); ++iBin ) + //~ { + //~ // Obtain norm and phase from components + //~ norm = std::hypot( (*ts)(iBin)[0], (*ts)(iBin)[1] ); + //~ phase = atan2( (*ts)(iBin)[1], (*ts)(iBin)[0] ); - // Shift phase - phase -= q1 * (double)iBin * (double)iBin; + //~ // Shift phase + //~ phase -= q1 * (double)iBin * (double)iBin; - // Assign components from norm and new phase - (tsDup)(iBin)[0] = norm * cos( phase ); - (tsDup)(iBin)[1] = norm * sin( phase ); - } + //~ // Assign components from norm and new phase + //~ (tsDup)(iBin)[0] = norm * cos( phase ); + //~ (tsDup)(iBin)[1] = norm * sin( phase ); + //~ } + + KTTimeSeriesFFTW tsDup{ *ts }; //copy ts + + // Shift phase + auto phaseTs = -q1 * square(Eigen::ArrayXd::LinSpaced( + ts->GetNTimeBins(), + 0, + ts->GetNTimeBins()-1)); + + tsDup.GetData() *= exp(std::complex{0.0, 1.0}*phaseTs); // Forward FFT KTFrequencySpectrumFFTW* fs = fForwardFFT.Transform( &tsDup ); @@ -200,10 +214,12 @@ namespace Katydid //~ } // Shift phase - auto bins = Eigen::ArrayXd::LinSpaced(fs->GetNFrequencyBins(), 0, - fs->GetNFrequencyBins()-1); + auto phaseFs = -q2 * square(Eigen::ArrayXd::LinSpaced( + fs->GetNFrequencyBins(), + 0, + fs->GetNFrequencyBins()-1)); - fs->GetData() *= exp(-std::complex{0.0, 1.0}*q2*bins*bins); + fs->GetData() *= exp(std::complex{0.0, 1.0}*phaseFs); // Reverse FFT KTTimeSeriesFFTW* newTS = fReverseFFT.TransformToComplex( fs ); @@ -213,20 +229,28 @@ namespace Katydid delete fs; // cleanup of fs // Third chirp transform - for( unsigned iBin = 0; iBin < newTS->GetNTimeBins(); ++iBin ) - { - // Obtain norm and phase from components - norm = std::hypot( (*newTS)(iBin)[0], (*newTS)(iBin)[1] ); - phase = atan2( (*newTS)(iBin)[1], (*newTS)(iBin)[0] ); + //~ for( unsigned iBin = 0; iBin < newTS->GetNTimeBins(); ++iBin ) + //~ { + //~ // Obtain norm and phase from components + //~ norm = std::hypot( (*newTS)(iBin)[0], (*newTS)(iBin)[1] ); + //~ phase = atan2( (*newTS)(iBin)[1], (*newTS)(iBin)[0] ); - // Shift phase - phase -= q1 * (double)iBin * (double)iBin; + //~ // Shift phase + //~ phase -= q1 * (double)iBin * (double)iBin; - // Assign components from norm and new phase - (*newTS)(iBin)[0] = norm * cos( phase ); - (*newTS)(iBin)[1] = norm * sin( phase ); - (*newFS)(iBin) = std::polar(norm, phase); - } + //~ // Assign components from norm and new phase + //~ (*newTS)(iBin)[0] = norm * cos( phase ); + //~ (*newTS)(iBin)[1] = norm * sin( phase ); + //~ (*newFS)(iBin) = std::polar(norm, phase); + //~ } + + // Shift phase + auto phaseNewTS = -q1 * square(Eigen::ArrayXd::LinSpaced( + newTS->GetNTimeBins(), + 0, + newTS->GetNTimeBins()-1)); + + newTS->GetData() *= exp(std::complex{0.0, 1.0}*phaseNewTS); newTSData.SetTimeSeries( newTS, iComponent ); // newTS now owned by newTSData newFSData.SetSpectrum( newFS, iComponent ); // newFS now owned by newFSData diff --git a/Source/Transform/KTReverseFFTW.cc b/Source/Transform/KTReverseFFTW.cc index 68c49b655..a3ee4926a 100644 --- a/Source/Transform/KTReverseFFTW.cc +++ b/Source/Transform/KTReverseFFTW.cc @@ -366,11 +366,14 @@ namespace Katydid void KTReverseFFTW::DoTransform(const KTFrequencySpectrumFFTW* fsIn, KTTimeSeriesFFTW* tsOut) const { - fftw_complex *data = + fftw_complex *dataIn = const_cast( reinterpret_cast( fsIn->GetData().data())); - fftw_execute_dft(fReversePlan, data, tsOut->GetData()); + + fftw_complex *dataOut = reinterpret_cast( + tsOut->GetData().data()); + fftw_execute_dft(fReversePlan, dataIn, dataOut); (*tsOut) *= sqrt(1. / double(fTimeSize)); return; } From 0642fe7ab9c9d54bb583115641778ffe97ab1286 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 21 Apr 2021 20:20:18 +0200 Subject: [PATCH 20/48] Modify files depending on KTTimeSeriesFFTW Mostly modifies a few lines per file that applied operations specific to fftw_complex --- Source/Executables/Main/RSAMatToEgg.cc | 7 ++++++- Source/Executables/Validation/TestComboFFTW.cc | 3 +-- Source/Executables/Validation/TestReverseFFTW.cc | 4 ++-- Source/IO/Conversions/KT2ROOT.cc | 6 +++--- Source/SpectrumAnalysis/KTDataAccumulator.cc | 12 +++++++----- Source/SpectrumAnalysis/KTWignerVille.cc | 6 ++++-- Source/SpectrumAnalysis/KTWignerVille.hh | 2 +- Source/Time/KTRSAMatReader.cc | 10 ++++++---- Source/Time/KTSingleChannelDAC.hh | 8 ++++---- Source/Transform/KTWindower.cc | 3 +-- 10 files changed, 35 insertions(+), 26 deletions(-) diff --git a/Source/Executables/Main/RSAMatToEgg.cc b/Source/Executables/Main/RSAMatToEgg.cc index db5c2f7c6..0ad67f94e 100644 --- a/Source/Executables/Main/RSAMatToEgg.cc +++ b/Source/Executables/Main/RSAMatToEgg.cc @@ -26,6 +26,8 @@ #include #include +#include + using namespace std; using namespace Katydid; @@ -182,7 +184,10 @@ int main(int argc, char** argv) // side note: if the time series array type were guaranteed to be contiguous in memory, this could be replaced with a memcpy for (unsigned iBin = 0; iBin < recSize; ++iBin) { - writer.set_at((*timeSeries)(iBin), iBin); + fftw_complex val; + val[0] = timeSeries->GetReal(iBin); + val[1] = timeSeries->GetImag(iBin); + writer.set_at(val, iBin); } stream->WriteRecord(firstRec); diff --git a/Source/Executables/Validation/TestComboFFTW.cc b/Source/Executables/Validation/TestComboFFTW.cc index 87b969d2f..6f435a7ae 100644 --- a/Source/Executables/Validation/TestComboFFTW.cc +++ b/Source/Executables/Validation/TestComboFFTW.cc @@ -50,8 +50,7 @@ int main() // The units are volts. for (unsigned iBin=0; iBinGetBinCenter(iBin) * mult); - (*timeSeries)(iBin)[1] = 0.0; + timeSeries->SetRect(iBin, cos(timeSeries->GetBinCenter(iBin) * mult), 0.0); } // Create and prepare the FFTs diff --git a/Source/Executables/Validation/TestReverseFFTW.cc b/Source/Executables/Validation/TestReverseFFTW.cc index fc9e98d70..8d365e77c 100644 --- a/Source/Executables/Validation/TestReverseFFTW.cc +++ b/Source/Executables/Validation/TestReverseFFTW.cc @@ -76,7 +76,7 @@ int main() for (unsigned iBin = 0; iBin < nTimeBins; iBin++) { - value = (*timeSeries)(iBin)[0]*(*timeSeries)(iBin)[0] + (*timeSeries)(iBin)[1]*(*timeSeries)(iBin)[1]; + value = timeSeries->GetAbs(iBin); if (value > maxValue) { maxValue = value; @@ -112,7 +112,7 @@ int main() double tsSum = 0.; // units: volts^2 for (unsigned iBin=0; iBinGetAbs(iBin); } KTINFO(vallog, "sum(timeSeries[i]^2) = " << tsSum << " V^2"); diff --git a/Source/IO/Conversions/KT2ROOT.cc b/Source/IO/Conversions/KT2ROOT.cc index 8beaa847d..c54a94c03 100644 --- a/Source/IO/Conversions/KT2ROOT.cc +++ b/Source/IO/Conversions/KT2ROOT.cc @@ -185,7 +185,7 @@ namespace Katydid TH1D* hist = new TH1D(histName.c_str(), "Time Series", (int) nBins, ts->GetRangeMin(), ts->GetRangeMax()); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, ::sqrt((*ts)(iBin)[0] * (*ts)(iBin)[0] + (*ts)(iBin)[1] * (*ts)(iBin)[1])); + hist->SetBinContent((int) iBin + 1, ts->GetAbs(iBin)); } hist->SetXTitle("Time (s)"); hist->SetYTitle("Voltage (V)"); @@ -198,7 +198,7 @@ namespace Katydid TH1D* hist = new TH1D(histName.c_str(), "Time Series (Real)", (int) nBins, ts->GetRangeMin(), ts->GetRangeMax()); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, (*ts)(iBin)[0]); + hist->SetBinContent((int) iBin + 1, ts->GetReal(iBin)); } hist->SetXTitle("Time (s)"); hist->SetYTitle("Voltage (V)"); @@ -211,7 +211,7 @@ namespace Katydid TH1D* hist = new TH1D(histName.c_str(), "Time Series (Imag)", (int) nBins, ts->GetRangeMin(), ts->GetRangeMax()); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, (*ts)(iBin)[1]); + hist->SetBinContent((int) iBin + 1, ts->GetImag(iBin)); } hist->SetXTitle("Time (s)"); hist->SetYTitle("Voltage (V)"); diff --git a/Source/SpectrumAnalysis/KTDataAccumulator.cc b/Source/SpectrumAnalysis/KTDataAccumulator.cc index 5f0ec951f..0e2777649 100644 --- a/Source/SpectrumAnalysis/KTDataAccumulator.cc +++ b/Source/SpectrumAnalysis/KTDataAccumulator.cc @@ -249,11 +249,13 @@ namespace Katydid { KTTimeSeriesFFTW* newTS = static_cast< KTTimeSeriesFFTW* >(data.GetTimeSeries(iComponent)); KTTimeSeriesFFTW* avTS = static_cast< KTTimeSeriesFFTW* >(accData.GetTimeSeries(iComponent)); - for (unsigned iBin = 0; iBin < arraySize; ++iBin) - { - (*avTS)(iBin)[0] = (*avTS)(iBin)[0] * remainingFrac + (*newTS)(iBin)[0] * fAveragingFrac; - (*avTS)(iBin)[1] = (*avTS)(iBin)[1] * remainingFrac + (*newTS)(iBin)[1] * fAveragingFrac; - } + //~ for (unsigned iBin = 0; iBin < arraySize; ++iBin) + //~ { + //~ (*avTS)(iBin)[0] = (*avTS)(iBin)[0] * remainingFrac + (*newTS)(iBin)[0] * fAveragingFrac; + //~ (*avTS)(iBin)[1] = (*avTS)(iBin)[1] * remainingFrac + (*newTS)(iBin)[1] * fAveragingFrac; + //~ } + + avTS->GetData() = avTS->GetData()*remainingFrac + newTS->GetData()*fAveragingFrac; } return true; diff --git a/Source/SpectrumAnalysis/KTWignerVille.cc b/Source/SpectrumAnalysis/KTWignerVille.cc index 1e3d72627..7abe30745 100644 --- a/Source/SpectrumAnalysis/KTWignerVille.cc +++ b/Source/SpectrumAnalysis/KTWignerVille.cc @@ -237,8 +237,10 @@ namespace Katydid t1_imag = data1It->imag(); t2_real = data2It->real(); t2_imag = data2It->imag(); - (*fInputArray)(fftBin)[0] = t1_real * t2_real + t1_imag * t2_imag; - (*fInputArray)(fftBin)[1] = t1_imag * t2_real - t1_real * t2_imag; + + fInputArray->SetRect(fftBin, + t1_real * t2_real + t1_imag * t2_imag, + t1_imag * t2_real - t1_real * t2_imag); //KTWARN(wvlog, " " << fftBin << " -- " << t1_real << " " << t1_imag << " -- " << t2_real << " " << t2_imag << " -- " << (*fInputArray)(fftBin)[0] << " " << (*fInputArray)(fftBin)[1]); ++data1It; } diff --git a/Source/SpectrumAnalysis/KTWignerVille.hh b/Source/SpectrumAnalysis/KTWignerVille.hh index 2b9c0af7e..a25cbda97 100644 --- a/Source/SpectrumAnalysis/KTWignerVille.hh +++ b/Source/SpectrumAnalysis/KTWignerVille.hh @@ -348,7 +348,7 @@ namespace Katydid unsigned tsSize = ts->size(); for (unsigned iBin = 0; iBin < tsSize; ++iBin) { - fBuffer[iComponent].push_back(std::complex< double >((*ts)(iBin)[0], (*ts)(iBin)[1])); + fBuffer[iComponent].push_back((*ts)(iBin)); } // we should only need to advance the start iterator if the start of this window // didn't fit in the last slice during the previous iteration diff --git a/Source/Time/KTRSAMatReader.cc b/Source/Time/KTRSAMatReader.cc index d134828e5..172748cfb 100644 --- a/Source/Time/KTRSAMatReader.cc +++ b/Source/Time/KTRSAMatReader.cc @@ -490,8 +490,9 @@ namespace Katydid float* dataImag = (float*)((mat_complex_split_t*)fTSArrayMat->data)->Im; for (unsigned iBin = 0; iBin < fSliceSize; iBin++) { - (*newSliceComplex)(iBin)[0] = double(dataReal[iBin + fSamplesRead]); - (*newSliceComplex)(iBin)[1] = double(dataImag[iBin + fSamplesRead]); + newSliceComplex->SetRect(iBin, + double(dataReal[iBin + fSamplesRead]), + double(dataImag[iBin + fSamplesRead])); } } else if (fDataPrecision == 2) @@ -500,8 +501,9 @@ namespace Katydid double* dataImag = (double*)((mat_complex_split_t*)fTSArrayMat->data)->Im; for (unsigned iBin = 0; iBin < fSliceSize; iBin++) { - (*newSliceComplex)(iBin)[0] = dataReal[iBin + fSamplesRead]; - (*newSliceComplex)(iBin)[1] = dataImag[iBin + fSamplesRead]; + newSliceComplex->SetRect(iBin, + dataReal[iBin + fSamplesRead], + dataImag[iBin + fSamplesRead]); } } else diff --git a/Source/Time/KTSingleChannelDAC.hh b/Source/Time/KTSingleChannelDAC.hh index e8ab243b0..251ad7171 100644 --- a/Source/Time/KTSingleChannelDAC.hh +++ b/Source/Time/KTSingleChannelDAC.hh @@ -206,8 +206,7 @@ namespace Katydid KTTimeSeriesFFTW* newTS = new KTTimeSeriesFFTW(nBins, ts.GetRangeMin(), ts.GetRangeMax()); for (unsigned bin = 0; bin < nBins; ++bin) { - (*newTS)(bin)[0] = Convert(ts(2 * bin)); - (*newTS)(bin)[1] = Convert(ts(2 * bin + 1)); + newTS->SetRect(bin, Convert(ts(2 * bin)), Convert(ts(2 * bin + 1))); } return newTS; } @@ -276,8 +275,9 @@ namespace Katydid avgValueImag += Convert(ts(2 * bin + 1)); ++bin; } - (*newTS)(oversampledBin)[0] = avgValueReal * fOversamplingScaleFactor; - (*newTS)(oversampledBin)[1] = avgValueImag * fOversamplingScaleFactor; + newTS->SetRect(oversampledBin, + avgValueReal * fOversamplingScaleFactor, + avgValueImag * fOversamplingScaleFactor); } #ifndef NDEBUG if (bin != ts.size() / 2) diff --git a/Source/Transform/KTWindower.cc b/Source/Transform/KTWindower.cc index 3e3461841..a85e06aa9 100644 --- a/Source/Transform/KTWindower.cc +++ b/Source/Transform/KTWindower.cc @@ -185,8 +185,7 @@ namespace Katydid for (unsigned iBin=0; iBin < nBins; ++iBin) { weight = fWindowFunction->GetWeight(iBin); - (*ts)(iBin)[0] = (*ts)(iBin)[0] * weight; - (*ts)(iBin)[1] = (*ts)(iBin)[1] * weight; + ts->SetRect(iBin, ts->GetReal(iBin)*weight, ts->GetImag(iBin)*weight); } return true; From 486e965a0b55f8a3b41f5558ffa4a7ae36991232 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 21 Apr 2021 20:26:37 +0200 Subject: [PATCH 21/48] Fix bug in TestReverseFFTW --- Source/Executables/Validation/TestReverseFFTW.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Executables/Validation/TestReverseFFTW.cc b/Source/Executables/Validation/TestReverseFFTW.cc index 8d365e77c..e88c6fa28 100644 --- a/Source/Executables/Validation/TestReverseFFTW.cc +++ b/Source/Executables/Validation/TestReverseFFTW.cc @@ -112,7 +112,7 @@ int main() double tsSum = 0.; // units: volts^2 for (unsigned iBin=0; iBinGetAbs(iBin); + tsSum += timeSeries->GetNorm(iBin); } KTINFO(vallog, "sum(timeSeries[i]^2) = " << tsSum << " V^2"); From ad47fefa93bc15af3c4c878b7062137db7963d52 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Thu, 22 Apr 2021 12:55:22 +0200 Subject: [PATCH 22/48] Fix some some wrong usage of absolute values in place of norm --- Source/Executables/Validation/TestReverseFFTW.cc | 2 +- Source/IO/Conversions/KT2ROOT.cc | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Executables/Validation/TestReverseFFTW.cc b/Source/Executables/Validation/TestReverseFFTW.cc index e88c6fa28..fbd9e30ae 100644 --- a/Source/Executables/Validation/TestReverseFFTW.cc +++ b/Source/Executables/Validation/TestReverseFFTW.cc @@ -76,7 +76,7 @@ int main() for (unsigned iBin = 0; iBin < nTimeBins; iBin++) { - value = timeSeries->GetAbs(iBin); + value = timeSeries->GetNorm(iBin); if (value > maxValue) { maxValue = value; diff --git a/Source/IO/Conversions/KT2ROOT.cc b/Source/IO/Conversions/KT2ROOT.cc index c54a94c03..a6f543348 100644 --- a/Source/IO/Conversions/KT2ROOT.cc +++ b/Source/IO/Conversions/KT2ROOT.cc @@ -421,7 +421,7 @@ namespace Katydid for (unsigned iBin = 0; iBin < nBins; ++iBin) { - hist->SetBinContent((int) iBin + 1, scaling * fs->GetAbs(iBin)); + hist->SetBinContent((int) iBin + 1, scaling * fs->GetNorm(iBin)); } hist->SetXTitle("Frequency (Hz)"); @@ -463,7 +463,7 @@ namespace Katydid // skip the DC bin; start at iBin = 1 for (unsigned iBin = 1; iBin < nBins; ++iBin) { - value = fs->GetAbs(iBin) * scaling; + value = fs->GetNorm(iBin) * scaling; if (value < tMinMag) tMinMag = value; if (value > tMaxMag) tMaxMag = value; } @@ -471,7 +471,7 @@ namespace Katydid TH1D* hist = new TH1D(name.c_str(), "Power Distribution", 100, tMinMag * 0.95, tMaxMag * 1.05); for (unsigned iBin = 0; iBin < nBins; ++iBin) { - value = fs->GetAbs(iBin); + value = fs->GetNorm(iBin); hist->Fill(value * scaling); } hist->SetXTitle("Power (W)"); From 3da235e50a6610b4fbde9df0e2c91919a4d58742 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 27 Apr 2021 16:46:01 +0200 Subject: [PATCH 23/48] Fix missing sign in KTConvolution.cc --- Source/SpectrumAnalysis/KTConvolution.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/SpectrumAnalysis/KTConvolution.cc b/Source/SpectrumAnalysis/KTConvolution.cc index e46b0baf5..fde95a2cc 100644 --- a/Source/SpectrumAnalysis/KTConvolution.cc +++ b/Source/SpectrumAnalysis/KTConvolution.cc @@ -446,7 +446,7 @@ namespace Katydid auto tmpImag = spectrum.GetImag(iBin); spectrum.SetRect(iBin, spectrum.GetReal(nBins - iBin - 1), -1. * spectrum.GetImag(nBins - iBin - 1)); - spectrum.SetRect(nBins - iBin -1, tmpReal, tmpImag); + spectrum.SetRect(nBins - iBin -1, tmpReal, -1. * tmpImag); } return; From bcc4863071cab5f3c4d67e151ee111aa63485e04 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 27 Apr 2021 16:53:32 +0200 Subject: [PATCH 24/48] Remove comments with old code --- Source/Data/Time/KTTimeSeriesFFTW.cc | 18 ---------- Source/Data/Time/KTTimeSeriesFFTW.hh | 6 +--- .../Data/Transform/KTFrequencySpectrumFFTW.cc | 34 ------------------- .../Data/Transform/KTFrequencySpectrumFFTW.hh | 23 +------------ .../KTFrequencyCandidateIdentifier.cc | 2 +- .../KTGainVariationProcessor.cc | 1 - .../KTSpectrumDiscriminator.cc | 4 +-- .../KTVariableSpectrumDiscriminator.cc | 6 ++-- 8 files changed, 8 insertions(+), 86 deletions(-) diff --git a/Source/Data/Time/KTTimeSeriesFFTW.cc b/Source/Data/Time/KTTimeSeriesFFTW.cc index 66d938ae2..3d20e779e 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.cc +++ b/Source/Data/Time/KTTimeSeriesFFTW.cc @@ -47,24 +47,6 @@ namespace Katydid } - // rule of zero, the default ones should do the job - - //~ KTTimeSeriesFFTW::KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig) : - //~ KTTimeSeries(), - //~ KTPhysicalArray< 1, std::complex >(orig) - //~ { - //~ } - - //~ KTTimeSeriesFFTW::~KTTimeSeriesFFTW() - //~ { - //~ } - - //~ KTTimeSeriesFFTW& KTTimeSeriesFFTW::operator=(const KTTimeSeriesFFTW& rhs) - //~ { - //~ KTPhysicalArray< 1, std::complex >::operator=(rhs); - //~ return *this; - //~ } - void KTTimeSeriesFFTW::Print(unsigned startPrint, unsigned nToPrint) const { stringstream printStream; diff --git a/Source/Data/Time/KTTimeSeriesFFTW.hh b/Source/Data/Time/KTTimeSeriesFFTW.hh index 61eda4b5b..aeb59b8f8 100644 --- a/Source/Data/Time/KTTimeSeriesFFTW.hh +++ b/Source/Data/Time/KTTimeSeriesFFTW.hh @@ -24,12 +24,8 @@ namespace Katydid KTTimeSeriesFFTW(); KTTimeSeriesFFTW(size_t nBins, double rangeMin=0., double rangeMax=1.); KTTimeSeriesFFTW(std::initializer_list value, size_t nBins, double rangeMin=0., double rangeMax=1.); - - //rule of zero - //KTTimeSeriesFFTW(const KTTimeSeriesFFTW& orig); - virtual ~KTTimeSeriesFFTW() = default; - //KTTimeSeriesFFTW& operator=(const KTTimeSeriesFFTW& rhs); + virtual ~KTTimeSeriesFFTW() = default; virtual void Scale(double scale); diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc index d1bb95799..eedee51f0 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.cc +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.cc @@ -72,40 +72,6 @@ namespace Katydid //std::copy(value.begin(), value.end(), this->begin()); } - //copy constructor, destructor and copy assignment operator - //shouldn't be necessary. I think the default ones will do the right thing - - //~ KTFrequencySpectrumFFTW::KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig) : - //~ KTPhysicalArray< 1, std::complex >(orig), - //~ KTFrequencySpectrum(), - //~ fIsArrayOrderFlipped(orig.fIsArrayOrderFlipped), - //~ fIsSizeEven(orig.fIsSizeEven), - //~ fLeftOfCenterOffset(orig.fLeftOfCenterOffset), - //~ fCenterBin(orig.fCenterBin), - //~ fConstBinAccess(orig.fConstBinAccess), - //~ fBinAccess(orig.fBinAccess), - //~ fNTimeBins(orig.fNTimeBins), - //~ fPointCache() - //~ { - //~ } - - //~ KTFrequencySpectrumFFTW::~KTFrequencySpectrumFFTW() - //~ { - //~ } - - //~ KTFrequencySpectrumFFTW& KTFrequencySpectrumFFTW::operator=(const KTFrequencySpectrumFFTW& rhs) - //~ { - //~ KTPhysicalArray< 1, std::complex >::operator=(rhs); - //~ fIsArrayOrderFlipped = rhs.fIsArrayOrderFlipped; - //~ fIsSizeEven = rhs.fIsSizeEven; - //~ fLeftOfCenterOffset = rhs.fLeftOfCenterOffset; - //~ fCenterBin = rhs.fCenterBin; - //~ fConstBinAccess = rhs.fConstBinAccess; - //~ fBinAccess = rhs.fBinAccess; - //~ fNTimeBins = rhs.fNTimeBins; - //~ return *this; - //~ } - const KTAxisProperties< 1 >& KTFrequencySpectrumFFTW::GetAxis() const { return *this; diff --git a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh index 5c977019c..baae9b402 100644 --- a/Source/Data/Transform/KTFrequencySpectrumFFTW.hh +++ b/Source/Data/Transform/KTFrequencySpectrumFFTW.hh @@ -26,8 +26,7 @@ namespace Katydid KTFrequencySpectrumFFTW(); KTFrequencySpectrumFFTW(size_t nBins, double rangeMin=0., double rangeMax=1., bool arrayOrderIsFlipped=false); KTFrequencySpectrumFFTW(std::initializer_list value, size_t nBins, double rangeMin=0., double rangeMax=1., bool arrayOrderIsFlipped=false); - //see comment in cc file - //KTFrequencySpectrumFFTW(const KTFrequencySpectrumFFTW& orig); + virtual ~KTFrequencySpectrumFFTW() = default; public: @@ -87,9 +86,6 @@ namespace Katydid public: // normal KTFrequencySpectrumPolar functions - //see comment in cc file - //virtual KTFrequencySpectrumFFTW& operator=(const KTFrequencySpectrumFFTW& rhs); - /// In-place calculation of the complex conjugate virtual KTFrequencySpectrumFFTW& CConjugate(); /// In-place calculation of the analytic associate @@ -176,44 +172,27 @@ namespace Katydid inline void KTFrequencySpectrumFFTW::SetRect(unsigned bin, double real, double imag) { - //~ fPointCache = &(*this)(bin); - //~ (*const_cast< fftw_complex* >(fPointCache))[0] = real; - //~ (*const_cast< fftw_complex* >(fPointCache))[1] = imag; - - //having the point cache for this really is not important - // for >=O2 the assembly code should be exactly the same - // can check it in https://godbolt.org/ (*this)(bin) = std::complex(real, imag); return; } inline double KTFrequencySpectrumFFTW::GetAbs(unsigned bin) const { - //fPointCache = &(*this)(bin); - //return sqrt((*fPointCache)[0]*(*fPointCache)[0] + (*fPointCache)[1]*(*fPointCache)[1]); return std::abs((*this)(bin)); } inline double KTFrequencySpectrumFFTW::GetNorm(unsigned bin) const { - //fPointCache = &(*this)(bin); - //return sqrt((*fPointCache)[0]*(*fPointCache)[0] + (*fPointCache)[1]*(*fPointCache)[1]); return std::norm((*this)(bin)); } inline double KTFrequencySpectrumFFTW::GetArg(unsigned bin) const { - //fPointCache = &(*this)(bin); - //return atan2((*fPointCache)[1], (*fPointCache)[0]); return std::arg((*this)(bin)); } inline void KTFrequencySpectrumFFTW::SetPolar(unsigned bin, double abs, double arg) { - //~ fPointCache = &(*this)(bin); - //~ (*const_cast< fftw_complex* >(fPointCache))[0] = abs * cos(arg); - //~ (*const_cast< fftw_complex* >(fPointCache))[1] = abs * sin(arg); - (*this)(bin) = std::polar(abs, arg); return; } diff --git a/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc b/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc index d6028b731..5905177de 100644 --- a/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc +++ b/Source/EventAnalysis/KTFrequencyCandidateIdentifier.cc @@ -150,7 +150,7 @@ namespace Katydid double peakValue = 0.; for (unsigned iBin=firstBin; iBin <= lastBin; iBin++) { - value = freqSpec->GetAbs(iBin); // sqrt((*freqSpec)(iBin)[0] * (*freqSpec)(iBin)[0] + (*freqSpec)(iBin)[1] * (*freqSpec)(iBin)[1]); + value = freqSpec->GetAbs(iBin); weightedMean += freqSpec->GetBinCenter(iBin) * value; integral += value; if (value > peakValue) peakValue = value; diff --git a/Source/SpectrumAnalysis/KTGainVariationProcessor.cc b/Source/SpectrumAnalysis/KTGainVariationProcessor.cc index 759195ce1..46736626a 100644 --- a/Source/SpectrumAnalysis/KTGainVariationProcessor.cc +++ b/Source/SpectrumAnalysis/KTGainVariationProcessor.cc @@ -462,7 +462,6 @@ namespace Katydid double mean = 0.; for (unsigned iBin=fitPointStartBin; iBinGetAbs(iBin); } mean /= (double)nBinsPerFitPoint; diff --git a/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc b/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc index 90e4300ed..461487389 100644 --- a/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc +++ b/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc @@ -223,7 +223,7 @@ namespace Katydid #pragma omp parallel for reduction(+:mean) for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { - magnitude[iBin] = spectrum->GetAbs(iBin); //sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); + magnitude[iBin] = spectrum->GetAbs(iBin); mean += magnitude[iBin]; variance += magnitude[iBin] * magnitude[iBin]; } @@ -472,7 +472,7 @@ namespace Katydid neighborhoodAmplitude = 0; for (unsigned jBin = iBin-fNeighborhoodRadius; jBin<= iBin+fNeighborhoodRadius; ++jBin) { - neighborhoodAmplitude += spectrum->GetAbs(jBin); //sqrt((*spectrum)(jBin)[0] * (*spectrum)(jBin)[0] + (*spectrum)(jBin)[1] * (*spectrum)(jBin)[1]); + neighborhoodAmplitude += spectrum->GetAbs(jBin); } } void KTSpectrumDiscriminator::SumAdjacentBinAmplitude(const KTFrequencySpectrumPolar* spectrum, double& neighborhoodAmplitude, const unsigned& iBin) diff --git a/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc b/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc index c60a15821..8a77e6834 100644 --- a/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc +++ b/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc @@ -536,7 +536,7 @@ namespace Katydid #pragma omp parallel for private(value) for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { - double value = spectrum->GetAbs(iBin); //sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); + double value = spectrum->GetAbs(iBin); double threshold = thresholdMult * (*splineImp)(iBin - fMinBin); double mean = (*splineImp)(iBin - fMinBin); double variance = (*varSplineImp)(iBin - fMinBin); @@ -572,7 +572,7 @@ namespace Katydid double mean = (*splineImp)(iBin - fMinBin); double variance = (*varSplineImp)(iBin - fMinBin); double threshold = mean + fSigmaThreshold * sqrt( variance ); - double value = spectrum->GetAbs(iBin); // sqrt((*spectrum)(iBin)[0] * (*spectrum)(iBin)[0] + (*spectrum)(iBin)[1] * (*spectrum)(iBin)[1]); + double value = spectrum->GetAbs(iBin); if (value >= threshold) { double neighborhoodAmplitude = 0.; @@ -714,7 +714,7 @@ namespace Katydid neighborhoodAmplitude = 0; for (unsigned jBin = iBin-fNeighborhoodRadius; jBin <= iBin+fNeighborhoodRadius; ++jBin) { - neighborhoodAmplitude += spectrum->GetAbs(jBin); //sqrt((*spectrum)(jBin)[0] * (*spectrum)(jBin)[0] + (*spectrum)(jBin)[1] * (*spectrum)(jBin)[1]); + neighborhoodAmplitude += spectrum->GetAbs(jBin); } } void KTVariableSpectrumDiscriminator::SumAdjacentBinAmplitude(const KTFrequencySpectrumPolar* spectrum, double& neighborhoodAmplitude, const unsigned& iBin) From ff3a8d5e76c0b8a832fefb53101f719108e231e8 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 28 Apr 2021 14:05:36 +0200 Subject: [PATCH 25/48] Add new files --- Source/Data/CMakeLists.txt | 2 + .../KTMultiChannelInnerProductData.cc | 188 ++++++++++++++++++ .../KTMultiChannelInnerProductData.hh | 106 ++++++++++ .../TestMultiChannelInnerProduct.cc | 78 ++++++++ Source/SpectrumAnalysis/CMakeLists.txt | 2 + .../KTMultiChannelInnerProduct.cc | 174 ++++++++++++++++ .../KTMultiChannelInnerProduct.hh | 104 ++++++++++ 7 files changed, 654 insertions(+) create mode 100644 Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc create mode 100644 Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh create mode 100644 Source/Executables/Validation/TestMultiChannelInnerProduct.cc create mode 100644 Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc create mode 100644 Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index 8bab20a44..260409b48 100644 --- a/Source/Data/CMakeLists.txt +++ b/Source/Data/CMakeLists.txt @@ -23,6 +23,7 @@ set (DATA_HEADERFILES SpectrumAnalysis/KTTimeSeriesDistData.hh SpectrumAnalysis/KTWV2DData.hh SpectrumAnalysis/KTWignerVilleData.hh + SpectrumAnalysis/KTMultiChannelInnerProductData.hh EventAnalysis/KTClassifierResultsData.hh EventAnalysis/KTFrequencyCandidate.hh EventAnalysis/KTFrequencyCandidateData.hh @@ -89,6 +90,7 @@ set (DATA_SOURCEFILES SpectrumAnalysis/KTTimeSeriesDistData.cc SpectrumAnalysis/KTWV2DData.cc SpectrumAnalysis/KTWignerVilleData.cc + SpectrumAnalysis/KTMultiChannelInnerProductData.cc EventAnalysis/KTClassifierResultsData.cc EventAnalysis/KTFrequencyCandidate.cc EventAnalysis/KTFrequencyCandidateData.cc diff --git a/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc b/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc new file mode 100644 index 000000000..d16c7bf08 --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc @@ -0,0 +1,188 @@ +/* + * KTConvolvedSpectrumData.cc + * + * Created on: Aug 25, 2017 + * Author: ezayas + */ + +#include "KTConvolvedSpectrumData.hh" + + +namespace Katydid +{ + const std::string KTConvolvedPowerSpectrumData::sName("convolved-power-spectrum"); + + KTConvolvedPowerSpectrumData::KTConvolvedPowerSpectrumData() : + KTPowerSpectrumDataCore(), + KTExtensibleData() + { + } + + KTConvolvedPowerSpectrumData::~KTConvolvedPowerSpectrumData() + { + } + + KTConvolvedPowerSpectrumData& KTConvolvedPowerSpectrumData::SetNComponents(unsigned num) + { + unsigned oldSize = fSpectra.size(); + // if num < oldSize + for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + fSpectra.resize(num); + // if num > oldSize + for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + return *this; + } + + const std::string KTConvolvedFrequencySpectrumDataFFTW::sName("convolved-frequency-spectrum-fftw"); + + KTConvolvedFrequencySpectrumDataFFTW::KTConvolvedFrequencySpectrumDataFFTW() : + KTFrequencySpectrumDataFFTWCore(), + KTExtensibleData() + { + } + + KTConvolvedFrequencySpectrumDataFFTW::~KTConvolvedFrequencySpectrumDataFFTW() + { + } + + KTConvolvedFrequencySpectrumDataFFTW& KTConvolvedFrequencySpectrumDataFFTW::SetNComponents(unsigned num) + { + unsigned oldSize = fSpectra.size(); + // if num < oldSize + for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + fSpectra.resize(num); + // if num > oldSize + for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + return *this; + } + + const std::string KTConvolvedFrequencySpectrumDataPolar::sName("convolved-frequency-spectrum-polar"); + + KTConvolvedFrequencySpectrumDataPolar::KTConvolvedFrequencySpectrumDataPolar() : + KTFrequencySpectrumDataPolarCore(), + KTExtensibleData() + { + } + + KTConvolvedFrequencySpectrumDataPolar::~KTConvolvedFrequencySpectrumDataPolar() + { + } + + KTConvolvedFrequencySpectrumDataPolar& KTConvolvedFrequencySpectrumDataPolar::SetNComponents(unsigned num) + { + unsigned oldSize = fSpectra.size(); + // if num < oldSize + for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + fSpectra.resize(num); + // if num > oldSize + for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + return *this; + } + + const std::string KTConvolvedPowerSpectrumVarianceData::sName("convolved-power-spectrum-variance"); + + KTConvolvedPowerSpectrumVarianceData::KTConvolvedPowerSpectrumVarianceData() : + KTFrequencySpectrumVarianceDataCore(), + KTExtensibleData() + { + } + + KTConvolvedPowerSpectrumVarianceData::~KTConvolvedPowerSpectrumVarianceData() + { + } + + KTConvolvedPowerSpectrumVarianceData& KTConvolvedPowerSpectrumVarianceData::SetNComponents(unsigned num) + { + unsigned oldSize = fSpectra.size(); + // if num < oldSize + for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + fSpectra.resize(num); + // if num > oldSize + for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + return *this; + } + + const std::string KTConvolvedFrequencySpectrumVarianceDataFFTW::sName("convolved-frequency-spectrum-variance-fftw"); + + KTConvolvedFrequencySpectrumVarianceDataFFTW::KTConvolvedFrequencySpectrumVarianceDataFFTW() : + KTFrequencySpectrumVarianceDataCore(), + KTExtensibleData() + { + } + + KTConvolvedFrequencySpectrumVarianceDataFFTW::~KTConvolvedFrequencySpectrumVarianceDataFFTW() + { + } + + KTConvolvedFrequencySpectrumVarianceDataFFTW& KTConvolvedFrequencySpectrumVarianceDataFFTW::SetNComponents(unsigned num) + { + unsigned oldSize = fSpectra.size(); + // if num < oldSize + for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + fSpectra.resize(num); + // if num > oldSize + for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + return *this; + } + + const std::string KTConvolvedFrequencySpectrumVarianceDataPolar::sName("convolved-frequency-spectrum-variance-polar"); + + KTConvolvedFrequencySpectrumVarianceDataPolar::KTConvolvedFrequencySpectrumVarianceDataPolar() : + KTFrequencySpectrumVarianceDataCore(), + KTExtensibleData() + { + } + + KTConvolvedFrequencySpectrumVarianceDataPolar::~KTConvolvedFrequencySpectrumVarianceDataPolar() + { + } + + KTConvolvedFrequencySpectrumVarianceDataPolar& KTConvolvedFrequencySpectrumVarianceDataPolar::SetNComponents(unsigned num) + { + unsigned oldSize = fSpectra.size(); + // if num < oldSize + for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + fSpectra.resize(num); + // if num > oldSize + for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + return *this; + } + + +} /* namespace Katydid */ diff --git a/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh b/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh new file mode 100644 index 000000000..1f9b8b170 --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh @@ -0,0 +1,106 @@ +/* + * KTMultiChannelInnerProductData.hh + * + * Created on: Apr 28, 2021 + * Author: F. Thomas + */ + +#ifndef KTMULTICHANNELINNERPRODUCTDATA_HH_ +#define KTMULTICHANNELINNERPRODUCTDATA_HH_ + +#include "KTData.hh" + +#include "KTPowerSpectrum.hh" +#include "KTPowerSpectrumData.hh" +#include "KTFrequencySpectrumFFTW.hh" +#include "KTFrequencySpectrumDataFFTW.hh" +#include "KTFrequencySpectrumPolar.hh" +#include "KTFrequencySpectrumDataPolar.hh" + +#include + +namespace Katydid +{ + + class KTConvolvedPowerSpectrumData : public KTPowerSpectrumDataCore, public Nymph::KTExtensibleData< KTConvolvedPowerSpectrumData > + { + public: + KTConvolvedPowerSpectrumData(); + virtual ~KTConvolvedPowerSpectrumData(); + + KTConvolvedPowerSpectrumData& SetNComponents(unsigned channels); + + public: + static const std::string sName; + + }; + + class KTConvolvedFrequencySpectrumDataFFTW : public KTFrequencySpectrumDataFFTWCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumDataFFTW > + { + public: + KTConvolvedFrequencySpectrumDataFFTW(); + virtual ~KTConvolvedFrequencySpectrumDataFFTW(); + + KTConvolvedFrequencySpectrumDataFFTW& SetNComponents(unsigned channels); + + public: + static const std::string sName; + + }; + + class KTConvolvedFrequencySpectrumDataPolar : public KTFrequencySpectrumDataPolarCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumDataPolar > + { + public: + KTConvolvedFrequencySpectrumDataPolar(); + virtual ~KTConvolvedFrequencySpectrumDataPolar(); + + KTConvolvedFrequencySpectrumDataPolar& SetNComponents(unsigned channels); + + public: + static const std::string sName; + + }; + + class KTConvolvedPowerSpectrumVarianceData : public KTFrequencySpectrumVarianceDataCore, public Nymph::KTExtensibleData< KTConvolvedPowerSpectrumVarianceData > + { + public: + KTConvolvedPowerSpectrumVarianceData(); + virtual ~KTConvolvedPowerSpectrumVarianceData(); + + KTConvolvedPowerSpectrumVarianceData& SetNComponents(unsigned channels); + + public: + static const std::string sName; + + }; + + class KTConvolvedFrequencySpectrumVarianceDataFFTW : public KTFrequencySpectrumVarianceDataCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumVarianceDataFFTW > + { + public: + KTConvolvedFrequencySpectrumVarianceDataFFTW(); + virtual ~KTConvolvedFrequencySpectrumVarianceDataFFTW(); + + KTConvolvedFrequencySpectrumVarianceDataFFTW& SetNComponents(unsigned channels); + + public: + static const std::string sName; + + }; + + class KTConvolvedFrequencySpectrumVarianceDataPolar : public KTFrequencySpectrumVarianceDataCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumVarianceDataPolar > + { + public: + KTConvolvedFrequencySpectrumVarianceDataPolar(); + virtual ~KTConvolvedFrequencySpectrumVarianceDataPolar(); + + KTConvolvedFrequencySpectrumVarianceDataPolar& SetNComponents(unsigned channels); + + public: + static const std::string sName; + + }; + + +} /* namespace Katydid */ + +#endif /* KTMULTICHANNELINNERPRODUCTDATA_HH_ */ diff --git a/Source/Executables/Validation/TestMultiChannelInnerProduct.cc b/Source/Executables/Validation/TestMultiChannelInnerProduct.cc new file mode 100644 index 000000000..03ef6cd7a --- /dev/null +++ b/Source/Executables/Validation/TestMultiChannelInnerProduct.cc @@ -0,0 +1,78 @@ +/* + * TestTemplate.cc + * + * Created on: Jan 26, 2021 + * Author: N.S. Oblath + * + * Tests performance of KTProcessorTemplate + * + * Usage: > TestTemplate + */ + +#include "KTProcessorTemplate.hh" + +#include "KTDummyDataClass1.hh" +#include "KTDummyDataClass2.hh" +#include "KTDummyDataClass3.hh" + +#include "KTNewDummyDataClass.hh" + +#include "KTLogger.hh" + +using namespace std; +using namespace Katydid; + +KTLOGGER(testlog, "TestTemplate"); + +int main() +{ + // Create and setup processor + KTProcessorTemplate tProc; + + // TODO: Set parameters in tProc + + + // Test processing of KTDummyDataClass1 + + KTDummyDataClass1 tDDC1; + + // TODO: Fill tDDC1 + + tProc.AnalyzeDummyClass1(tDDC1); + + // Check Results + KTNewDummyDataClass& tNDDC1 = tDDC1.Of< KTNewDummyDataClass >(); + + // TODO: Verify that the contents of tNDDC1 are as expected + + + // Test processing of KTDummyDataClass2 + + KTDummyDataClass2 tDDC2; + + // TODO: Fill tDDC2 + + tProc.AnalyzeDummyClass2(tDDC2); + + // Check Results + KTNewDummyDataClass& tNDDC2 = tDDC2.Of< KTNewDummyDataClass >(); + + // TODO: Verify that the contents of tNDDC2 are as expected + + + // Test processing of KTDummyDataClass3 + + KTDummyDataClass3 tDDC3; + + // TODO: Fill tDDC3 + + tProc.AnalyzeDummyClass3(tDDC3); + + // Check Results + KTNewDummyDataClass& tNDDC3 = tDDC3.Of< KTNewDummyDataClass >(); + + // TODO: Verify that the contents of tNDDC3 are as expected + + + return 0; +} diff --git a/Source/SpectrumAnalysis/CMakeLists.txt b/Source/SpectrumAnalysis/CMakeLists.txt index e1dec2264..febd78aa5 100644 --- a/Source/SpectrumAnalysis/CMakeLists.txt +++ b/Source/SpectrumAnalysis/CMakeLists.txt @@ -24,6 +24,7 @@ set (SPECTRUMANALYSIS_HEADERFILES KTSpectrogramStriper.hh KTSwitchFFTWPolar.hh KTVariableSpectrumDiscriminator.hh + KTMultiChannelInnerProduct.hh #KTAutoCorrMatrix.hh ) @@ -50,6 +51,7 @@ set (SPECTRUMANALYSIS_SOURCEFILES KTSpectrogramStriper.cc KTSwitchFFTWPolar.cc KTVariableSpectrumDiscriminator.cc + KTMultiChannelInnerProduct.cc #KTAutoCorrMatrix.cc ) diff --git a/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc b/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc new file mode 100644 index 000000000..e2c3f4e3d --- /dev/null +++ b/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc @@ -0,0 +1,174 @@ +/* + * KTMultiChannelInnerProduct.cc + * + * Created on: Apr 28, 2021 + * Author: F. Thomas + */ + +#include "KTLowPassFilter.hh" + +#include "KTFrequencySpectrumDataPolar.hh" +#include "KTFrequencySpectrumDataFFTW.hh" +#include "KTFrequencySpectrumFFTW.hh" +#include "KTFrequencySpectrumPolar.hh" +#include "KTLogger.hh" +#include "KTLowPassFilteredData.hh" +#include "KTMath.hh" +#include "KTPowerSpectrum.hh" +#include "KTPowerSpectrumData.hh" + +#include "param.hh" + + +namespace Katydid +{ + KTLOGGER(gclog, "KTLowPassFilter"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTLowPassFilter, "low-pass-filter"); + + KTLowPassFilter::KTLowPassFilter(const std::string& name) : + KTProcessor(name), + fRC(0.), + fFSPolarSignal("fs-polar", this), + fFSFFTWSignal("fs-fftw", this), + fPSSignal("ps", this), + fFSPolarSlot("fs-polar", this, &KTLowPassFilter::Filter, &fFSPolarSignal), + fFSFFTWSlot("fs-fftw", this, &KTLowPassFilter::Filter, &fFSFFTWSignal), + fPSSlot("ps", this, &KTLowPassFilter::Filter, &fPSSignal) + { + } + + KTLowPassFilter::~KTLowPassFilter() + { + } + + bool KTLowPassFilter::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + SetRC(node->get_value("rc", GetRC())); + + return true; + } + + bool KTLowPassFilter::Filter(KTFrequencySpectrumDataPolar& fsData) + { + unsigned nComponents = fsData.GetNComponents(); + KTLowPassFilteredFSDataPolar& newData = fsData.Of< KTLowPassFilteredFSDataPolar >().SetNComponents(nComponents); + + for (unsigned iComponent=0; iComponent().SetNComponents(nComponents); + + for (unsigned iComponent=0; iComponent().SetNComponents(nComponents); + + for (unsigned iComponent=0; iComponentsize(); + KTFrequencySpectrumPolar* newSpectrum = new KTFrequencySpectrumPolar(nBins, frequencySpectrum->GetRangeMin(), frequencySpectrum->GetRangeMax()); + newSpectrum->SetNTimeBins(frequencySpectrum->GetNTimeBins()); + + double twoPiRCf, abs, arg; + for (unsigned iBin = 0; iBin < nBins; ++iBin) + { + abs = (*frequencySpectrum)(iBin).abs(); + arg = (*frequencySpectrum)(iBin).arg(); + twoPiRCf = KTMath::TwoPi() * fRC * frequencySpectrum->GetBinCenter(iBin); + (*newSpectrum)(iBin).set_polar(abs / sqrt(1 + twoPiRCf*twoPiRCf), arg - atan(twoPiRCf)); + } + + return newSpectrum; + } + + KTFrequencySpectrumFFTW* KTLowPassFilter::Filter(const KTFrequencySpectrumFFTW* frequencySpectrum) const + { + KTDEBUG(gclog, "Creating new FS for filtered data"); + unsigned nBins = frequencySpectrum->size(); + KTFrequencySpectrumFFTW* newSpectrum = new KTFrequencySpectrumFFTW(nBins, frequencySpectrum->GetRangeMin(), frequencySpectrum->GetRangeMax()); + newSpectrum->SetNTimeBins(frequencySpectrum->GetNTimeBins()); + + double twoPiRCf, real, imag, denom, newReal, newImag; + for (unsigned iBin = 0; iBin < nBins; ++iBin) + { + real = (*frequencySpectrum)(iBin)[0]; + imag = (*frequencySpectrum)(iBin)[1]; + twoPiRCf = KTMath::TwoPi() * fRC * frequencySpectrum->GetBinCenter(iBin); + denom = 1 + twoPiRCf * twoPiRCf; + (*newSpectrum)(iBin)[0] = (real + imag * twoPiRCf) / denom; + (*newSpectrum)(iBin)[1] = (imag - real * twoPiRCf) / denom; + } + + return newSpectrum; + } + + KTPowerSpectrum* KTLowPassFilter::Filter(const KTPowerSpectrum* powerSpectrum) const + { + KTDEBUG(gclog, "Creating new PS for filtered data"); + unsigned nBins = powerSpectrum->size(); + KTPowerSpectrum* newSpectrum = new KTPowerSpectrum(nBins, powerSpectrum->GetRangeMin(), powerSpectrum->GetRangeMax()); + + double twoPiRCf; + for (unsigned iBin = 0; iBin < nBins; ++iBin) + { + twoPiRCf = KTMath::TwoPi() * fRC * powerSpectrum->GetBinCenter(iBin); + (*newSpectrum)(iBin) = (*powerSpectrum)(iBin) / (1 + twoPiRCf * twoPiRCf); + } + + return newSpectrum; + } + +} /* namespace Katydid */ diff --git a/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh b/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh new file mode 100644 index 000000000..b37dfc1eb --- /dev/null +++ b/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh @@ -0,0 +1,104 @@ +/** + @file KTMultiChannelInnerProduct.hh + @brief Contains KTMultiChannelInnerProduct + @details Applies the Matched Filter inner product for multiple channels + @author: F. Thomas + @date: Apr 28, 2021 + */ + +#ifndef KTMULTICHANNELINNERPRODUCT_HH_ +#define KTMULTICHANNELINNERPRODUCT_HH_ + +#include "KTProcessor.hh" + +#include "KTMemberVariable.hh" +#include "KTSlot.hh" + +#include + +namespace scarab +{ + class param_node; +} + +namespace Katydid +{ + class KTFrequencySpectrumDataFFTW; + class KTFrequencySpectrumDataPolar; + class KTFrequencySpectrumFFTW; + class KTFrequencySpectrumPolar; + class KTPowerSpectrum; + class KTPowerSpectrumData; + + /*! + @class KTMultiChannelInnerProduct + @author F. Thomas + + @brief Applies the Matched Filter inner product for multiple channels + + @details + Frequency-domain implementation of a single-pole RC low-pass filter. + This is no-doubt a non-ideal implementation, but demonstrates the features of a processor very well. + + The relationship between the cutoff frequency, f_c, and the RC constant is: + ``` + RC = 1 / 2*pi*f_c + ``` + + Configuration name: "low-pass-filter" + + Available configuration values: + - "rc": double -- RC time constant of the filter + + Slots: + - "fs-polar": void (KTDataPtr) -- Applies a low-pass filter; Requires KTFrequencySpectrumDataPolar; Adds KTLowPassFilteredFSDataPolar; Emits signal "fs-polar" + - "fs-fftw": void (KTDataPtr) -- Applies a low-pass filter; Requires KTFrequencySpectrumDataFFTW; Adds KTLowPassFilteredFSDataFFTW; Emits signal "fs-fftw" + - "ps": void (KTDataPtr) -- Applies a low-pass filter; Requires KTPowerSpectrum; Adds KTLowPassFilteredPSData; Emits signal "ps" + + Signals: + - "fs-polar": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredFSDataPolar. + - "fs-fftw": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredFSDataFFTW. + - "ps": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredPSData. + */ + class KTMultiChannelInnerProduct : public Nymph::KTProcessor + { + public: + KTMultiChannelInnerProduct(const std::string& name = "multi-channel-inner-product"); + virtual ~KTMultiChannelInnerProduct(); + + bool Configure(const scarab::param_node* node); + + MEMBERVARIABLE(double, RC); + + + public: + bool Filter(KTFrequencySpectrumDataPolar& fsData); + bool Filter(KTFrequencySpectrumDataFFTW& fsData); + bool Filter(KTPowerSpectrumData& psData); + + KTFrequencySpectrumPolar* Filter(const KTFrequencySpectrumPolar* frequencySpectrum) const; + KTFrequencySpectrumFFTW* Filter(const KTFrequencySpectrumFFTW* frequencySpectrum) const; + KTPowerSpectrum* Filter(const KTPowerSpectrum* powerSpectrum) const; + + //*************** + // Signals + //*************** + + private: + Nymph::KTSignalData fFSPolarSignal; + Nymph::KTSignalData fFSFFTWSignal; + Nymph::KTSignalData fPSSignal; + + //*************** + // Slots + //*************** + + private: + Nymph::KTSlotDataOneType< KTFrequencySpectrumDataPolar > fFSPolarSlot; + Nymph::KTSlotDataOneType< KTFrequencySpectrumDataFFTW > fFSFFTWSlot; + Nymph::KTSlotDataOneType< KTPowerSpectrumData > fPSSlot; + + }; + +} /* namespace Katydid */ +#endif /* KTMULTICHANNELINNERPRODUCT_HH_ */ From fd65779347220a486099c54d023364383de42d0d Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Sat, 1 May 2021 13:34:09 +0200 Subject: [PATCH 26/48] Add move operations to KTPhysicalArray --- Source/Utility/KTPhysicalArrayComplex.cc | 2 -- Source/Utility/KTPhysicalArrayComplex.hh | 15 +++++++++++++-- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/Source/Utility/KTPhysicalArrayComplex.cc b/Source/Utility/KTPhysicalArrayComplex.cc index b74271a46..8b919b262 100644 --- a/Source/Utility/KTPhysicalArrayComplex.cc +++ b/Source/Utility/KTPhysicalArrayComplex.cc @@ -39,8 +39,6 @@ namespace Katydid fData.fill(value); } - KTPhysicalArray< 1, value_type >::~KTPhysicalArray() - {} const Eigen::Array< std::complex, Eigen::Dynamic, 1, Eigen::ColMajor >& KTPhysicalArray< 1, std::complex >::GetData() const { diff --git a/Source/Utility/KTPhysicalArrayComplex.hh b/Source/Utility/KTPhysicalArrayComplex.hh index 37e8cce77..0c39c4611 100644 --- a/Source/Utility/KTPhysicalArrayComplex.hh +++ b/Source/Utility/KTPhysicalArrayComplex.hh @@ -92,7 +92,12 @@ namespace Katydid explicit KTPhysicalArray(size_t nBins, double rangeMin=0., double rangeMax=1.); explicit KTPhysicalArray(value_type value, size_t nBins, double rangeMin=0., double rangeMax=1.); - virtual ~KTPhysicalArray(); + virtual ~KTPhysicalArray() = default; + KTPhysicalArray( KTPhysicalArray< 1, value_type> &&) = default; + KTPhysicalArray& operator=( KTPhysicalArray< 1, value_type> &&) = default; + + KTPhysicalArray( const KTPhysicalArray< 1, value_type>& ) = default; + KTPhysicalArray operator=( const KTPhysicalArray< 1, value_type>& ) = default; public: const array_type& GetData() const; @@ -220,7 +225,13 @@ namespace Katydid KTPhysicalArray(); KTPhysicalArray(size_t xNBins, double xRangeMin, double xRangeMax, size_t yNBins, double yRangeMin, double yRangeMax); KTPhysicalArray(value_type value, size_t xNBins, double xRangeMin, double xRangeMax, size_t yNBins, double yRangeMin, double yRangeMax); - virtual ~KTPhysicalArray(); + virtual ~KTPhysicalArray() = default; + + KTPhysicalArray( KTPhysicalArray< 2, value_type> &&) = default; + KTPhysicalArray& operator=( KTPhysicalArray< 2, value_type> &&) = default; + + KTPhysicalArray( const KTPhysicalArray< 2, value_type>& ) = default; + KTPhysicalArray operator=( const KTPhysicalArray< 2, value_type>& ) = default; public: const matrix_type& GetData() const; From b45c46da92668d023b9e1bab225f035c35d7cfa6 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Sat, 1 May 2021 16:45:49 +0200 Subject: [PATCH 27/48] Add new processor KTMatrixAggregator --- Source/Data/CMakeLists.txt | 6 +- .../Transform/KTAggregatedTSMatrixData.cc | 15 ++ .../Transform/KTAggregatedTSMatrixData.hh | 35 +++++ Source/SpectrumAnalysis/CMakeLists.txt | 4 +- Source/Transform/CMakeLists.txt | 2 + Source/Transform/KTMatrixAggregator.cc | 141 ++++++++++++++++++ Source/Transform/KTMatrixAggregator.hh | 89 +++++++++++ 7 files changed, 288 insertions(+), 4 deletions(-) create mode 100644 Source/Data/Transform/KTAggregatedTSMatrixData.cc create mode 100644 Source/Data/Transform/KTAggregatedTSMatrixData.hh create mode 100644 Source/Transform/KTMatrixAggregator.cc create mode 100644 Source/Transform/KTMatrixAggregator.hh diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index 260409b48..b630c782f 100644 --- a/Source/Data/CMakeLists.txt +++ b/Source/Data/CMakeLists.txt @@ -23,7 +23,7 @@ set (DATA_HEADERFILES SpectrumAnalysis/KTTimeSeriesDistData.hh SpectrumAnalysis/KTWV2DData.hh SpectrumAnalysis/KTWignerVilleData.hh - SpectrumAnalysis/KTMultiChannelInnerProductData.hh + # SpectrumAnalysis/KTMultiChannelInnerProductData.hh EventAnalysis/KTClassifierResultsData.hh EventAnalysis/KTFrequencyCandidate.hh EventAnalysis/KTFrequencyCandidateData.hh @@ -69,6 +69,7 @@ set (DATA_HEADERFILES Transform/KTTimeFrequency.hh Transform/KTTimeFrequencyDataPolar.hh Transform/KTTimeFrequencyPolar.hh + Transform/KTAggregatedTSMatrixData.hh ) set (DATA_SOURCEFILES @@ -90,7 +91,7 @@ set (DATA_SOURCEFILES SpectrumAnalysis/KTTimeSeriesDistData.cc SpectrumAnalysis/KTWV2DData.cc SpectrumAnalysis/KTWignerVilleData.cc - SpectrumAnalysis/KTMultiChannelInnerProductData.cc + # SpectrumAnalysis/KTMultiChannelInnerProductData.cc EventAnalysis/KTClassifierResultsData.cc EventAnalysis/KTFrequencyCandidate.cc EventAnalysis/KTFrequencyCandidateData.cc @@ -136,6 +137,7 @@ set (DATA_SOURCEFILES Transform/KTTimeFrequency.cc Transform/KTTimeFrequencyDataPolar.cc Transform/KTTimeFrequencyPolar.cc + Transform/KTAggregatedTSMatrixData.cc ) if (NOT FFTW_FOUND) diff --git a/Source/Data/Transform/KTAggregatedTSMatrixData.cc b/Source/Data/Transform/KTAggregatedTSMatrixData.cc new file mode 100644 index 000000000..ad904d4a2 --- /dev/null +++ b/Source/Data/Transform/KTAggregatedTSMatrixData.cc @@ -0,0 +1,15 @@ +/* + * KTTimeSeriesMatrixData.cc + * + * Created on: Apr 29, 2021 + * Author: F. Thomas + */ + +#include "KTAggregatedTSMatrixData.hh" + + +namespace Katydid +{ + const std::string KTAggregatedTSMatrixData::sName("aggregated-ts-matrix"); + +} /* namespace Katydid */ diff --git a/Source/Data/Transform/KTAggregatedTSMatrixData.hh b/Source/Data/Transform/KTAggregatedTSMatrixData.hh new file mode 100644 index 000000000..6b9d0e9f6 --- /dev/null +++ b/Source/Data/Transform/KTAggregatedTSMatrixData.hh @@ -0,0 +1,35 @@ +/* + * KTAggregatedTSMatrixData.hh + * + * Created on: Apr 29, 2021 + * Author: F. Thomas + */ + +#ifndef KTAGGREGATEDTSMATRIXDATA_HH_ +#define KTAGGREGATEDTSMATRIXDATA_HH_ + +#include "KTData.hh" +#include "KTPhysicalArrayComplex.hh" + +namespace Katydid +{ + + class KTAggregatedTSMatrixData: public KTPhysicalArray<2, std::complex>, public Nymph::KTExtensibleData< KTAggregatedTSMatrixData > + { + public: + KTAggregatedTSMatrixData() = default; + + virtual ~KTAggregatedTSMatrixData() = default; + KTAggregatedTSMatrixData (KTAggregatedTSMatrixData &&) = default; + KTAggregatedTSMatrixData ( const KTAggregatedTSMatrixData &) = default; + + KTAggregatedTSMatrixData& operator=( const KTAggregatedTSMatrixData &) = default; + KTAggregatedTSMatrixData& operator=( KTAggregatedTSMatrixData && ) = default; + + static const std::string sName; + + }; + +} /* namespace Katydid */ + +#endif /* KTAGGREGATEDTSMATRIXDATA_HH_ */ diff --git a/Source/SpectrumAnalysis/CMakeLists.txt b/Source/SpectrumAnalysis/CMakeLists.txt index febd78aa5..ffa291964 100644 --- a/Source/SpectrumAnalysis/CMakeLists.txt +++ b/Source/SpectrumAnalysis/CMakeLists.txt @@ -24,7 +24,7 @@ set (SPECTRUMANALYSIS_HEADERFILES KTSpectrogramStriper.hh KTSwitchFFTWPolar.hh KTVariableSpectrumDiscriminator.hh - KTMultiChannelInnerProduct.hh + #KTMultiChannelInnerProduct.hh #KTAutoCorrMatrix.hh ) @@ -51,7 +51,7 @@ set (SPECTRUMANALYSIS_SOURCEFILES KTSpectrogramStriper.cc KTSwitchFFTWPolar.cc KTVariableSpectrumDiscriminator.cc - KTMultiChannelInnerProduct.cc + #KTMultiChannelInnerProduct.cc #KTAutoCorrMatrix.cc ) diff --git a/Source/Transform/CMakeLists.txt b/Source/Transform/CMakeLists.txt index 4127a62f5..688521c22 100644 --- a/Source/Transform/CMakeLists.txt +++ b/Source/Transform/CMakeLists.txt @@ -14,6 +14,7 @@ set (TRANSFORM_NODICT_HEADERFILES KTSincWindow.hh KTWindower.hh KTWindowFunction.hh + KTMatrixAggregator.hh ) if (FFTW_FOUND) @@ -37,6 +38,7 @@ set (TRANSFORM_SOURCEFILES KTSincWindow.cc KTWindower.cc KTWindowFunction.cc + KTMatrixAggregator.cc ) if (FFTW_FOUND) diff --git a/Source/Transform/KTMatrixAggregator.cc b/Source/Transform/KTMatrixAggregator.cc new file mode 100644 index 000000000..0458748a4 --- /dev/null +++ b/Source/Transform/KTMatrixAggregator.cc @@ -0,0 +1,141 @@ +/* + * KTMultiChannelInnerProduct.cc + * + * Created on: Apr 28, 2021 + * Author: F. Thomas + */ + +#include "KTMatrixAggregator.hh" + +#include "KTLogger.hh" +#include "KTMath.hh" +#include "KTTimeSeries.hh" +#include "KTTimeSeriesData.hh" +#include "KTTimeSeriesFFTW.hh" + +#include "KTDemangle.hh" + +#include "param.hh" + + +namespace Katydid +{ + KTLOGGER(magglog, "KTMatrixAggregator"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTMatrixAggregator, "matrix-aggregator"); + + KTMatrixAggregator::KTMatrixAggregator(const std::string& name) : + KTProcessor{name}, + fMaxCols{1}, + fSignalCount{0}, + fNRows{0}, + fBufferMat{}, + fMatrixSignal{"matrix", this} + //why does this not work? + //fTSSlot("ts-fftw", this, &KTMatrixAggregator::SlotFunction, &fMatrixSignal) + { + RegisterSlot("ts-fftw", this, &KTMatrixAggregator::SlotFunction); + } + + bool KTMatrixAggregator::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + fMaxCols = node->get_value< unsigned >("max-cols", fMaxCols); + + return true; + } + + void KTMatrixAggregator::ShrinkMatrix() + { + fBufferMat.conservativeResize(Eigen::NoChange_t{}, fSignalCount); + } + + bool KTMatrixAggregator::AddCol(KTTimeSeriesData& tsData) + { + const KTTimeSeriesFFTW* ts = dynamic_cast< const KTTimeSeriesFFTW* >(tsData.GetTimeSeries(0)); + //KTTimeSeries* ts = tsData.GetTimeSeries(0); + + unsigned nComponents = tsData.GetNComponents(); + unsigned nTimebins = ts->GetNTimeBins(); + unsigned rows = nComponents * nTimebins; + + if(fBufferMat.rows()== 0) + { + fBufferMat.resize(rows, fMaxCols); + } + + if(fBufferMat.rows() != rows) + { + KTERROR(magglog, "The input slice has not the same number of timebins or components as the first."); + return false; + } + + for (unsigned iComponent=0; iComponent(tsData.GetTimeSeries(iComponent)); + + if (ts->GetNTimeBins() != nTimebins) + { + KTERROR(magglog, "TS " << iComponent << " has not the same length as TS 0."); + return false; + } + + fBufferMat.block(nTimebins*iComponent, fSignalCount, nTimebins, 1) = ts->GetData(); + KTDEBUG(magglog, "Added TS " << iComponent << " to the matrix."); + + } + + + return true; + } + + void KTMatrixAggregator::SlotFunction(Nymph::KTDataPtr data) + { + + if (! data->Has< KTTimeSeriesData >()) + { + KTERROR(magglog, "Data not found with type <" << DemangledName(typeid(KTTimeSeriesData)) << ">"); + return; + } + + // Call the function + if (! AddCol(data->Of< KTTimeSeriesData >())) + { + KTERROR(magglog, "Could not add TS to matrix"); + return; + } + // Increase signal counter + fSignalCount++; + bool emitSignal = false; + if (fSignalCount == fMaxCols) + { + fSignalCount = 0; + KTINFO(magglog, "Completed the matrix."); + emitSignal = true; + } + + if (data->GetLastData()) + { + KTDEBUG(magglog, "Emitting last-data signal"); + //shrink the matrix to current signal count if end of data is reached + + ShrinkMatrix(); + KTINFO(magglog, "Completed the matrix."); + emitSignal = true; + } + + if(emitSignal) { + KTAggregatedTSMatrixData& aggMatrix = data->Of< KTAggregatedTSMatrixData >(); + //adjust labels and KTAxis things? + aggMatrix.GetData() = std::move(fBufferMat); + fMatrixSignal(data); + fBufferMat = Eigen::ArrayXcd(); + + } + + return; + } + +} /* namespace Katydid */ diff --git a/Source/Transform/KTMatrixAggregator.hh b/Source/Transform/KTMatrixAggregator.hh new file mode 100644 index 000000000..28018bd91 --- /dev/null +++ b/Source/Transform/KTMatrixAggregator.hh @@ -0,0 +1,89 @@ +/** + @file KTMatrixAggregator.hh + @brief Contains KTMatrixAggregator + @details Aggregates data from several egg files and components to a single matrix + @author: F. Thomas + @date: Apr 29, 2021 + */ + +#ifndef KTMATRIXAGGREGATOR_HH_ +#define KTMATRIXAGGREGATOR_HH_ + +#include "KTProcessor.hh" + +#include "KTMemberVariable.hh" +#include "KTSlot.hh" + +#include "KTAggregatedTSMatrixData.hh" + +#include +#include + +namespace scarab +{ + class param_node; +} + +namespace Katydid +{ + + class KTTimeSeriesData; + + /*! + @class KTMatrixAggregator + @author F. Thomas + + @brief Aggregates data from several egg files and components to a single matrix + + @details + Frequency-domain implementation of a single-pole RC low-pass filter. + This is no-doubt a non-ideal implementation, but demonstrates the features of a processor very well. + + The relationship between the cutoff frequency, f_c, and the RC constant is: + ``` + RC = 1 / 2*pi*f_c + ``` + + Configuration name: "low-pass-filter" + + Available configuration values: + - "rc": double -- RC time constant of the filter + + Slots: + - "ts-fftw": void (Nymph::KTDataPtr) -- Combine multiple multi-channel complex time series to a matrix; Requires KTTimeSeriesData; Adds KTAggregatedTSMatrixData; Emits signal "matrix" + Signals: + - "matrix": void (Nymph::KTDataPtr) -- Emitted upon matrix conversion; Guarantees KTLowPassFilteredPSData. + */ + class KTMatrixAggregator : public Nymph::KTProcessor + { + public: + KTMatrixAggregator(const std::string& name = "matrix-aggregator"); + virtual ~KTMatrixAggregator() = default; + + bool Configure(const scarab::param_node* node); + + MEMBERVARIABLE(int, MaxCols); + + private: + bool AddCol(KTTimeSeriesData& tsData); + void SlotFunction(Nymph::KTDataPtr data); + void ShrinkMatrix(); + + Eigen::ArrayXcd fBufferMat; + unsigned fNRows; + + unsigned fSignalCount; + + //*************** + // Signals + //*************** + Nymph::KTSignalData fMatrixSignal; + + //*************** + // Slots + //*************** + //Nymph::KTSlotDataOneType< Nymph::KTDataPtr > fTSSlot; + }; + +} /* namespace Katydid */ +#endif /* KTMATRIXAGGREGATOR_HH_ */ From 47dbf9dbdbbdbf6490ea3215957507df5660c876 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Sat, 1 May 2021 16:47:21 +0200 Subject: [PATCH 28/48] Fix small issues in KTPhysicalArrayComplex --- Source/Utility/KTPhysicalArrayComplex.cc | 5 ----- Source/Utility/KTPhysicalArrayComplex.hh | 8 ++++---- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/Source/Utility/KTPhysicalArrayComplex.cc b/Source/Utility/KTPhysicalArrayComplex.cc index 8b919b262..adfa33ab9 100644 --- a/Source/Utility/KTPhysicalArrayComplex.cc +++ b/Source/Utility/KTPhysicalArrayComplex.cc @@ -278,11 +278,6 @@ namespace Katydid { fData.fill(value); } - - - KTPhysicalArray< 2, std::complex >::~KTPhysicalArray() - { - } size_t KTPhysicalArray< 2, std::complex >::cols() const { diff --git a/Source/Utility/KTPhysicalArrayComplex.hh b/Source/Utility/KTPhysicalArrayComplex.hh index 0c39c4611..6cf09692c 100644 --- a/Source/Utility/KTPhysicalArrayComplex.hh +++ b/Source/Utility/KTPhysicalArrayComplex.hh @@ -94,10 +94,10 @@ namespace Katydid virtual ~KTPhysicalArray() = default; KTPhysicalArray( KTPhysicalArray< 1, value_type> &&) = default; - KTPhysicalArray& operator=( KTPhysicalArray< 1, value_type> &&) = default; + KTPhysicalArray<1, value_type>& operator=( KTPhysicalArray< 1, value_type> &&) = default; KTPhysicalArray( const KTPhysicalArray< 1, value_type>& ) = default; - KTPhysicalArray operator=( const KTPhysicalArray< 1, value_type>& ) = default; + KTPhysicalArray<1, value_type>& operator=( const KTPhysicalArray< 1, value_type>& ) = default; public: const array_type& GetData() const; @@ -228,10 +228,10 @@ namespace Katydid virtual ~KTPhysicalArray() = default; KTPhysicalArray( KTPhysicalArray< 2, value_type> &&) = default; - KTPhysicalArray& operator=( KTPhysicalArray< 2, value_type> &&) = default; + KTPhysicalArray<2, value_type>& operator=( KTPhysicalArray< 2, value_type> &&) = default; KTPhysicalArray( const KTPhysicalArray< 2, value_type>& ) = default; - KTPhysicalArray operator=( const KTPhysicalArray< 2, value_type>& ) = default; + KTPhysicalArray<2, value_type>& operator=( const KTPhysicalArray< 2, value_type>& ) = default; public: const matrix_type& GetData() const; From c7ce37c5c2885af8bda2cee31fb872002497a743 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 4 May 2021 15:23:50 +0200 Subject: [PATCH 29/48] Add Test for MatrixAggregator --- .../Validation/TestMatrixAggregator.cc | 73 +++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 Source/Executables/Validation/TestMatrixAggregator.cc diff --git a/Source/Executables/Validation/TestMatrixAggregator.cc b/Source/Executables/Validation/TestMatrixAggregator.cc new file mode 100644 index 000000000..bcf62adcb --- /dev/null +++ b/Source/Executables/Validation/TestMatrixAggregator.cc @@ -0,0 +1,73 @@ +/* + * TestMatrixAggregator.cc + * + * Created on: May 4, 2021 + * Author: F. Thomas + * + * Tests performance of KTMatrixAggregator + * + * Usage: > TestMatrixAggregator + */ +#include +#include + +#include "KTMatrixAggregator.hh" + +#include "KTTimeSeriesData.hh" +#include "KTTimeSeriesFFTW.hh" + +#include "KTAggregatedTSMatrixData.hh" + +#include "KTLogger.hh" + +using namespace std; +using namespace Katydid; + +KTLOGGER(testlog, "TestMatrixAggregator"); + +int main() +{ + // Create and setup processor + KTMatrixAggregator tAgg {}; + + unsigned maxCols = 3; + unsigned nCols = 5; + unsigned nChannels = 3; + unsigned nTimeBins = 4; + + tAgg.SetMaxCols(maxCols); + + // Test processing of KTTimeSeriesData + + Nymph::KTDataPtr data = boost::make_shared(); + for (int i=0; iHas()); + KTTimeSeriesData& tsData = data->Of< KTTimeSeriesData >(); + KTDEBUG("Make tsData success"); + tsData.SetNComponents(nChannels); + for (int j=0; jSetRect(k, i+j+2*k, i-j+k); + } + KTDEBUG("TS: " << *ts); + tsData.SetTimeSeries(ts, i); + } + KTDEBUG("Call the slot function"); + tAgg.SlotFunction(data); + } + + // Check Results + KTAggregatedTSMatrixData& dataMatrix = data->Of< KTAggregatedTSMatrixData >(); + + KTDEBUG("Aggregated: " << dataMatrix); + + return 0; +} From 2ac4eb02abbf4e9218e8e750d4f599f71c72fa01 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 4 May 2021 17:01:14 +0200 Subject: [PATCH 30/48] Fix issues in TestMatrixAggregator --- Source/Executables/Validation/CMakeLists.txt | 2 +- Source/Executables/Validation/TestMatrixAggregator.cc | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/Source/Executables/Validation/CMakeLists.txt b/Source/Executables/Validation/CMakeLists.txt index 5e1b8156f..8fccc5f20 100644 --- a/Source/Executables/Validation/CMakeLists.txt +++ b/Source/Executables/Validation/CMakeLists.txt @@ -88,7 +88,7 @@ if (Katydid_ENABLE_TESTING) TestSpectrumDiscriminator TestTrackProcessing TestWindowFunction - + TestMatrixAggregator ) pbuilder_executables( PROGRAMS LIB_DEPENDENCIES ) diff --git a/Source/Executables/Validation/TestMatrixAggregator.cc b/Source/Executables/Validation/TestMatrixAggregator.cc index bcf62adcb..4aa44ca0d 100644 --- a/Source/Executables/Validation/TestMatrixAggregator.cc +++ b/Source/Executables/Validation/TestMatrixAggregator.cc @@ -42,10 +42,11 @@ int main() Nymph::KTDataPtr data = boost::make_shared(); for (int i=0; iHas()); + if(i==nCols-1) + { + data->SetLastData(true); + } KTTimeSeriesData& tsData = data->Of< KTTimeSeriesData >(); - KTDEBUG("Make tsData success"); tsData.SetNComponents(nChannels); for (int j=0; jSetRect(k, i+j+2*k, i-j+k); } KTDEBUG("TS: " << *ts); - tsData.SetTimeSeries(ts, i); + tsData.SetTimeSeries(ts, j); } KTDEBUG("Call the slot function"); tAgg.SlotFunction(data); + tAgg.PrintBuffer(); } // Check Results From 4b898e041c68c578d0aeda2701abf37903ea84c9 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 4 May 2021 17:03:47 +0200 Subject: [PATCH 31/48] Add better Debug capabilities to KTMatrixAggregator --- Source/Transform/KTMatrixAggregator.cc | 19 +++++++++++++++---- Source/Transform/KTMatrixAggregator.hh | 7 +++++-- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/Source/Transform/KTMatrixAggregator.cc b/Source/Transform/KTMatrixAggregator.cc index 0458748a4..f181b4f03 100644 --- a/Source/Transform/KTMatrixAggregator.cc +++ b/Source/Transform/KTMatrixAggregator.cc @@ -52,6 +52,12 @@ namespace Katydid fBufferMat.conservativeResize(Eigen::NoChange_t{}, fSignalCount); } + void KTMatrixAggregator::PrintBuffer() + { + KTDEBUG(magglog, "Buffer matrix has size (" << fBufferMat.rows() << "," << fBufferMat.cols() << ")"); + KTDEBUG(magglog, "\n"<< fBufferMat); + } + bool KTMatrixAggregator::AddCol(KTTimeSeriesData& tsData) { const KTTimeSeriesFFTW* ts = dynamic_cast< const KTTimeSeriesFFTW* >(tsData.GetTimeSeries(0)); @@ -61,9 +67,12 @@ namespace Katydid unsigned nTimebins = ts->GetNTimeBins(); unsigned rows = nComponents * nTimebins; + KTDEBUG(magglog, "Buffer matrix has size (" << fBufferMat.rows() << "," << fBufferMat.cols() << ")"); if(fBufferMat.rows()== 0) { + KTDEBUG(magglog, "Resizing to (" << rows << "," << fMaxCols << ")" ); fBufferMat.resize(rows, fMaxCols); + KTDEBUG(magglog, "Resizing successful"); } if(fBufferMat.rows() != rows) @@ -83,8 +92,7 @@ namespace Katydid } fBufferMat.block(nTimebins*iComponent, fSignalCount, nTimebins, 1) = ts->GetData(); - KTDEBUG(magglog, "Added TS " << iComponent << " to the matrix."); - + KTDEBUG(magglog, "Added data successfully"); } @@ -111,6 +119,7 @@ namespace Katydid bool emitSignal = false; if (fSignalCount == fMaxCols) { + KTDEBUG(magglog, "Matrix full"); fSignalCount = 0; KTINFO(magglog, "Completed the matrix."); emitSignal = true; @@ -118,7 +127,7 @@ namespace Katydid if (data->GetLastData()) { - KTDEBUG(magglog, "Emitting last-data signal"); + KTDEBUG(magglog, "Reached last data"); //shrink the matrix to current signal count if end of data is reached ShrinkMatrix(); @@ -129,9 +138,11 @@ namespace Katydid if(emitSignal) { KTAggregatedTSMatrixData& aggMatrix = data->Of< KTAggregatedTSMatrixData >(); //adjust labels and KTAxis things? + KTDEBUG(magglog, "Emitting the signal"); aggMatrix.GetData() = std::move(fBufferMat); fMatrixSignal(data); - fBufferMat = Eigen::ArrayXcd(); + KTDEBUG(magglog, "Reset buffer matrix"); + fBufferMat = Eigen::ArrayXXcd(); } diff --git a/Source/Transform/KTMatrixAggregator.hh b/Source/Transform/KTMatrixAggregator.hh index 28018bd91..022c99138 100644 --- a/Source/Transform/KTMatrixAggregator.hh +++ b/Source/Transform/KTMatrixAggregator.hh @@ -64,12 +64,15 @@ namespace Katydid MEMBERVARIABLE(int, MaxCols); + void PrintBuffer(); + + void SlotFunction(Nymph::KTDataPtr data); + private: bool AddCol(KTTimeSeriesData& tsData); - void SlotFunction(Nymph::KTDataPtr data); void ShrinkMatrix(); - Eigen::ArrayXcd fBufferMat; + Eigen::ArrayXXcd fBufferMat; unsigned fNRows; unsigned fSignalCount; From ccbd8190a3b4735211f9053ad44bc6c96b086cdb Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 5 May 2021 14:22:33 +0200 Subject: [PATCH 32/48] Remove unnecessary destructors --- Source/Data/Transform/KTAggregatedTSMatrixData.cc | 2 +- Source/Data/Transform/KTAggregatedTSMatrixData.hh | 4 +++- Source/Transform/KTMatrixAggregator.hh | 1 - 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Source/Data/Transform/KTAggregatedTSMatrixData.cc b/Source/Data/Transform/KTAggregatedTSMatrixData.cc index ad904d4a2..132b2b722 100644 --- a/Source/Data/Transform/KTAggregatedTSMatrixData.cc +++ b/Source/Data/Transform/KTAggregatedTSMatrixData.cc @@ -1,5 +1,5 @@ /* - * KTTimeSeriesMatrixData.cc + * KTAggregatedTSMatrixData.cc * * Created on: Apr 29, 2021 * Author: F. Thomas diff --git a/Source/Data/Transform/KTAggregatedTSMatrixData.hh b/Source/Data/Transform/KTAggregatedTSMatrixData.hh index 6b9d0e9f6..a0df62742 100644 --- a/Source/Data/Transform/KTAggregatedTSMatrixData.hh +++ b/Source/Data/Transform/KTAggregatedTSMatrixData.hh @@ -17,7 +17,8 @@ namespace Katydid class KTAggregatedTSMatrixData: public KTPhysicalArray<2, std::complex>, public Nymph::KTExtensibleData< KTAggregatedTSMatrixData > { public: - KTAggregatedTSMatrixData() = default; + + /* KTAggregatedTSMatrixData() = default; virtual ~KTAggregatedTSMatrixData() = default; KTAggregatedTSMatrixData (KTAggregatedTSMatrixData &&) = default; @@ -25,6 +26,7 @@ namespace Katydid KTAggregatedTSMatrixData& operator=( const KTAggregatedTSMatrixData &) = default; KTAggregatedTSMatrixData& operator=( KTAggregatedTSMatrixData && ) = default; + */ static const std::string sName; diff --git a/Source/Transform/KTMatrixAggregator.hh b/Source/Transform/KTMatrixAggregator.hh index 022c99138..a6e3e18da 100644 --- a/Source/Transform/KTMatrixAggregator.hh +++ b/Source/Transform/KTMatrixAggregator.hh @@ -58,7 +58,6 @@ namespace Katydid { public: KTMatrixAggregator(const std::string& name = "matrix-aggregator"); - virtual ~KTMatrixAggregator() = default; bool Configure(const scarab::param_node* node); From daf8f93647c5f613dd36b405b03fac0045a7d1bd Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 5 May 2021 14:23:23 +0200 Subject: [PATCH 33/48] Add processor to convert timeseries to normalized MF templates --- Source/Data/CMakeLists.txt | 2 + .../KTAggregatedTemplateMatrixData.cc | 15 +++ .../KTAggregatedTemplateMatrixData.hh | 37 +++++++ Source/Transform/CMakeLists.txt | 2 + Source/Transform/KTConvertToTemplate.cc | 67 ++++++++++++ Source/Transform/KTConvertToTemplate.hh | 100 ++++++++++++++++++ 6 files changed, 223 insertions(+) create mode 100644 Source/Data/Transform/KTAggregatedTemplateMatrixData.cc create mode 100644 Source/Data/Transform/KTAggregatedTemplateMatrixData.hh create mode 100644 Source/Transform/KTConvertToTemplate.cc create mode 100644 Source/Transform/KTConvertToTemplate.hh diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index b630c782f..a17160444 100644 --- a/Source/Data/CMakeLists.txt +++ b/Source/Data/CMakeLists.txt @@ -70,6 +70,7 @@ set (DATA_HEADERFILES Transform/KTTimeFrequencyDataPolar.hh Transform/KTTimeFrequencyPolar.hh Transform/KTAggregatedTSMatrixData.hh + Transform/KTAggregatedTemplateMatrixData.hh ) set (DATA_SOURCEFILES @@ -138,6 +139,7 @@ set (DATA_SOURCEFILES Transform/KTTimeFrequencyDataPolar.cc Transform/KTTimeFrequencyPolar.cc Transform/KTAggregatedTSMatrixData.cc + Transform/KTAggregatedTemplateMatrixData.cc ) if (NOT FFTW_FOUND) diff --git a/Source/Data/Transform/KTAggregatedTemplateMatrixData.cc b/Source/Data/Transform/KTAggregatedTemplateMatrixData.cc new file mode 100644 index 000000000..1ad33a8f8 --- /dev/null +++ b/Source/Data/Transform/KTAggregatedTemplateMatrixData.cc @@ -0,0 +1,15 @@ +/* + * KTAggregatedTemplateMatrixData.cc + * + * Created on: May 4, 2021 + * Author: F. Thomas + */ + +#include "KTAggregatedTemplateMatrixData.hh" + + +namespace Katydid +{ + const std::string KTAggregatedTemplateMatrixData::sName("aggregated-template-matrix"); + +} /* namespace Katydid */ diff --git a/Source/Data/Transform/KTAggregatedTemplateMatrixData.hh b/Source/Data/Transform/KTAggregatedTemplateMatrixData.hh new file mode 100644 index 000000000..1e4c983d3 --- /dev/null +++ b/Source/Data/Transform/KTAggregatedTemplateMatrixData.hh @@ -0,0 +1,37 @@ +/* + * KTAggregatedTemplateMatrixData.hh + * + * Created on: May 4, 2021 + * Author: F. Thomas + */ + +#ifndef KTAGGREGATEDTEMPLATEMATRIXDATA_HH_ +#define KTAGGREGATEDTEMPLATEMATRIXDATA_HH_ + +#include "KTData.hh" +#include "KTPhysicalArrayComplex.hh" + +namespace Katydid +{ + + class KTAggregatedTemplateMatrixData: public KTPhysicalArray<2, std::complex>, public Nymph::KTExtensibleData< KTAggregatedTemplateMatrixData > + { + public: + /* + KTAggregatedTemplateMatrixData() = default; + + virtual ~KTAggregatedTemplateMatrixData() = default; + KTAggregatedTemplateMatrixData (KTAggregatedTemplateMatrixData &&) = default; + KTAggregatedTemplateMatrixData ( const KTAggregatedTemplateMatrixData &) = default; + + KTAggregatedTemplateMatrixData& operator=( const KTAggregatedTemplateMatrixData &) = default; + KTAggregatedTemplateMatrixData& operator=( KTAggregatedTemplateMatrixData && ) = default; + */ + + static const std::string sName; + + }; + +} /* namespace Katydid */ + +#endif /* KTAGGREGATEDTEMPLATEMATRIXDATA_HH_ */ diff --git a/Source/Transform/CMakeLists.txt b/Source/Transform/CMakeLists.txt index 688521c22..9daa0141c 100644 --- a/Source/Transform/CMakeLists.txt +++ b/Source/Transform/CMakeLists.txt @@ -15,6 +15,7 @@ set (TRANSFORM_NODICT_HEADERFILES KTWindower.hh KTWindowFunction.hh KTMatrixAggregator.hh + KTConvertToTemplate.hh ) if (FFTW_FOUND) @@ -39,6 +40,7 @@ set (TRANSFORM_SOURCEFILES KTWindower.cc KTWindowFunction.cc KTMatrixAggregator.cc + KTConvertToTemplate.cc ) if (FFTW_FOUND) diff --git a/Source/Transform/KTConvertToTemplate.cc b/Source/Transform/KTConvertToTemplate.cc new file mode 100644 index 000000000..aa5ea2219 --- /dev/null +++ b/Source/Transform/KTConvertToTemplate.cc @@ -0,0 +1,67 @@ +/* + * KTConvertToTemplate.cc + * + * Created on: May 4, 2021 + * Author: F. Thomas + */ + +#include "KTConvertToTemplate.hh" + +#include "KTLogger.hh" +#include "KTAggregatedTemplateMatrixData.hh" +#include "KTAggregatedTSMatrixData.hh" +#include "KTMath.hh" + +#include "param.hh" + + +namespace Katydid +{ + KTLOGGER(ctemplatelog, "KTConvertToTemplate"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTConvertToTemplate, "template-converter"); + + KTConvertToTemplate::KTConvertToTemplate(const std::string& name) : + KTProcessor(name), + fNoiseTemperature(0.), + fTSSignal("template-matrix", this), + fTSSlot("ts-matrix", this, &KTConvertToTemplate::Convert, &fTSSignal) + { + } + + void KTConvertToTemplate::CalcNoiseStd() + { + double kb = 1.380649e-23; + double R = 50.0; + // Johnson-Nyquist noise + fNoiseStd = sqrt(4*kb*R*fNoiseTemperature*fBandwidth); + } + + bool KTConvertToTemplate::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + SetNoiseTemperature(node->get_value("T", GetNoiseTemperature())); + SetBandwidth(node->get_value("bandwidth", GetBandwidth())); + + CalcNoiseStd(); + + return true; + } + + bool KTConvertToTemplate::Convert(KTAggregatedTSMatrixData& fData) + { + + KTAggregatedTemplateMatrixData& newData = fData.Of< KTAggregatedTemplateMatrixData >(); + + auto energy = sqrt((fData.GetData()*conj(fData.GetData())).colwise().sum()); + + auto normalized = fData.GetData()*sqrt(2)/(energy*fNoiseStd); + + newData.GetData() = normalized.transpose(); + + return true; + } + +} /* namespace Katydid */ diff --git a/Source/Transform/KTConvertToTemplate.hh b/Source/Transform/KTConvertToTemplate.hh new file mode 100644 index 000000000..48a99a5d3 --- /dev/null +++ b/Source/Transform/KTConvertToTemplate.hh @@ -0,0 +1,100 @@ +/** + @file KTConvertToTemplate.hh + @brief Contains KTConvertToTemplate + @details Converts timeseries data to a matched filter template + @author: F. Thomas + @date: May 4, 2021 + */ + +#ifndef KTCONVERTTOTEMPLATE_HH_ +#define KTCONVERTTOTEMPLATE_HH_ + +#include "KTProcessor.hh" + +#include "KTMemberVariable.hh" +#include "KTSlot.hh" + +#include + +namespace scarab +{ + class param_node; +} + +namespace Katydid +{ + + class KTAggregatedTSMatrixData; + + /*! + @class KTConvertToTemplate + @author F. Thomas + + @brief Converts data to Matched Filter templates + + @details + Converts a matrix of timeseries data to a matrix of normalized Matched Filter templates. + The normalization includes a noise contribution, which assumes white Gaussian noise. + Using this normalization the MF output will be a dimensionless SNR. + + For an input timeseries X the normalized template is + ``` + Y = 1/sqrt(X^H*C^-1*X) C^-1 X + ``` + where C is the noise covariance matrix. + With the white Gaussian assumption this simplifies to + ``` + Y = sqrt(2)/sqrt(X^2*V) * X + ``` + where V is the Johnson-Nyquist noise variance. + + Configuration name: "template-converter" + + Available configuration values: + - "T": double -- Noise temperature + - "bandwidth": double -- Bandwidth of the input data, i.e. half the sampling rate + + Slots: + - "ts-matrix": void (KTDataPtr) -- Converts timeseries to normalized MF template; Requires KTAggregatedTSMatrixData; Adds KTAggregatedTemplateMatrixDatar; Emits signal "template-matrix" + + Signals: + - "template-matrix": void (KTDataPtr) -- Emitted conversion; Guarantees KTAggregatedTemplateMatrixData. + + */ + class KTConvertToTemplate : public Nymph::KTProcessor + { + public: + KTConvertToTemplate(const std::string& name = "template-converter"); + + bool Configure(const scarab::param_node* node); + + MEMBERVARIABLE(double, NoiseTemperature); + MEMBERVARIABLE(double, Bandwidth); + + + public: + bool Convert(KTAggregatedTSMatrixData& fData); + + private: + + double fNoiseStd; + void CalcNoiseStd(); + + //*************** + // Signals + //*************** + + private: + Nymph::KTSignalData fTSSignal; + + //*************** + // Slots + //*************** + + private: + Nymph::KTSlotDataOneType< KTAggregatedTSMatrixData > fTSSlot; + + }; + +} /* namespace Katydid */ +#endif /* KTCONVERTTOTEMPLATE_HH_ */ From 0465a5d831b5a83fa772095236b4cbd1e302df07 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 5 May 2021 15:46:01 +0200 Subject: [PATCH 34/48] Add test for KTConvertToTemplate --- ...trixAggregator.cc => TestMatchedFilter.cc} | 34 ++++++++++++++----- 1 file changed, 26 insertions(+), 8 deletions(-) rename Source/Executables/Validation/{TestMatrixAggregator.cc => TestMatchedFilter.cc} (59%) diff --git a/Source/Executables/Validation/TestMatrixAggregator.cc b/Source/Executables/Validation/TestMatchedFilter.cc similarity index 59% rename from Source/Executables/Validation/TestMatrixAggregator.cc rename to Source/Executables/Validation/TestMatchedFilter.cc index 4aa44ca0d..a6c0a4239 100644 --- a/Source/Executables/Validation/TestMatrixAggregator.cc +++ b/Source/Executables/Validation/TestMatchedFilter.cc @@ -1,22 +1,24 @@ /* - * TestMatrixAggregator.cc + * TestMatchedFilter.cc * * Created on: May 4, 2021 * Author: F. Thomas * - * Tests performance of KTMatrixAggregator + * Tests KTMatrixAggregator, KTConvertToTemplate * - * Usage: > TestMatrixAggregator + * Usage: > TestMatchedFilter */ #include #include #include "KTMatrixAggregator.hh" +#include "KTConvertToTemplate.hh" #include "KTTimeSeriesData.hh" #include "KTTimeSeriesFFTW.hh" #include "KTAggregatedTSMatrixData.hh" +#include "KTAggregatedTemplateMatrixData.hh" #include "KTLogger.hh" @@ -27,16 +29,21 @@ KTLOGGER(testlog, "TestMatrixAggregator"); int main() { - // Create and setup processor + // Create and setup processors KTMatrixAggregator tAgg {}; + KTConvertToTemplate tConvert {}; unsigned maxCols = 3; unsigned nCols = 5; unsigned nChannels = 3; unsigned nTimeBins = 4; + double noiseStd = sqrt(2.0); + tAgg.SetMaxCols(maxCols); + tConvert.SetNoiseStd(noiseStd); + // Test processing of KTTimeSeriesData Nymph::KTDataPtr data = boost::make_shared(); @@ -50,7 +57,7 @@ int main() tsData.SetNComponents(nChannels); for (int j=0; jSetRect(k, i+j+2*k, i-j+k); } - KTDEBUG("TS: " << *ts); + KTDEBUG(testlog, "TS: " << *ts); tsData.SetTimeSeries(ts, j); } - KTDEBUG("Call the slot function"); + KTDEBUG(testlog, "Call the slot function"); tAgg.SlotFunction(data); tAgg.PrintBuffer(); } @@ -69,7 +76,18 @@ int main() // Check Results KTAggregatedTSMatrixData& dataMatrix = data->Of< KTAggregatedTSMatrixData >(); - KTDEBUG("Aggregated: " << dataMatrix); + KTDEBUG(testlog, "Aggregated: " << dataMatrix); + + // Test conversion to template + tConvert.Convert(dataMatrix); + + KTAggregatedTemplateMatrixData& templateMatrix = data->Of< KTAggregatedTemplateMatrixData >(); + + KTDEBUG(testlog, "Template Matrix: " << templateMatrix); + + KTDEBUG(testlog, "Norm" << sqrt((dataMatrix.GetData()*conj(dataMatrix.GetData())).colwise().sum())); + + KTDEBUG(testlog, " Product: " << (templateMatrix.GetData().transpose()*conj(dataMatrix.GetData())).colwise().sum()); return 0; } From 4f7973650dab0eb1693fbb552b6bd488aaaa5187 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 5 May 2021 15:46:53 +0200 Subject: [PATCH 35/48] Add TestMatchedFilter.cc to build --- Source/Executables/Validation/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Executables/Validation/CMakeLists.txt b/Source/Executables/Validation/CMakeLists.txt index 8fccc5f20..6e4d4a639 100644 --- a/Source/Executables/Validation/CMakeLists.txt +++ b/Source/Executables/Validation/CMakeLists.txt @@ -88,7 +88,7 @@ if (Katydid_ENABLE_TESTING) TestSpectrumDiscriminator TestTrackProcessing TestWindowFunction - TestMatrixAggregator + TestMatchedFilter ) pbuilder_executables( PROGRAMS LIB_DEPENDENCIES ) From 1ab57e033dbe4519ce80ceae80d2a3d18352204d Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 5 May 2021 15:49:07 +0200 Subject: [PATCH 36/48] Add small fixes and debug output to KTConvertToTemplate --- Source/Transform/KTConvertToTemplate.cc | 12 +++++++++++- Source/Transform/KTConvertToTemplate.hh | 2 +- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/Source/Transform/KTConvertToTemplate.cc b/Source/Transform/KTConvertToTemplate.cc index aa5ea2219..659b28743 100644 --- a/Source/Transform/KTConvertToTemplate.cc +++ b/Source/Transform/KTConvertToTemplate.cc @@ -36,6 +36,7 @@ namespace Katydid double R = 50.0; // Johnson-Nyquist noise fNoiseStd = sqrt(4*kb*R*fNoiseTemperature*fBandwidth); + KTDEBUG(ctemplatelog, "Noise standard deviation is: " << fNoiseStd); } bool KTConvertToTemplate::Configure(const scarab::param_node* node) @@ -55,12 +56,21 @@ namespace Katydid KTAggregatedTemplateMatrixData& newData = fData.Of< KTAggregatedTemplateMatrixData >(); + KTDEBUG(ctemplatelog, "Calculating energy"); auto energy = sqrt((fData.GetData()*conj(fData.GetData())).colwise().sum()); - auto normalized = fData.GetData()*sqrt(2)/(energy*fNoiseStd); + auto normalization = sqrt(2)/(energy*fNoiseStd); + KTDEBUG(ctemplatelog, "Calculating normalized matrix"); + auto normalized = fData.GetData().rowwise()*normalization; + + KTDEBUG(ctemplatelog, "Store transposed matrix in newData"); newData.GetData() = normalized.transpose(); + KTDEBUG(ctemplatelog, "Template matrix has shape (" + << newData.GetData().rows() + << "," << newData.GetData().cols() << ")" ); + return true; } diff --git a/Source/Transform/KTConvertToTemplate.hh b/Source/Transform/KTConvertToTemplate.hh index 48a99a5d3..e2575f0a8 100644 --- a/Source/Transform/KTConvertToTemplate.hh +++ b/Source/Transform/KTConvertToTemplate.hh @@ -70,6 +70,7 @@ namespace Katydid MEMBERVARIABLE(double, NoiseTemperature); MEMBERVARIABLE(double, Bandwidth); + MEMBERVARIABLE(double, NoiseStd); public: @@ -77,7 +78,6 @@ namespace Katydid private: - double fNoiseStd; void CalcNoiseStd(); //*************** From caf06bb541edba55cae0e805917dae9a581a0381 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 5 May 2021 23:35:35 +0200 Subject: [PATCH 37/48] Small fix for KTConvertToTemplate --- Source/Transform/KTConvertToTemplate.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Source/Transform/KTConvertToTemplate.cc b/Source/Transform/KTConvertToTemplate.cc index 659b28743..3571d78e4 100644 --- a/Source/Transform/KTConvertToTemplate.cc +++ b/Source/Transform/KTConvertToTemplate.cc @@ -57,14 +57,15 @@ namespace Katydid KTAggregatedTemplateMatrixData& newData = fData.Of< KTAggregatedTemplateMatrixData >(); KTDEBUG(ctemplatelog, "Calculating energy"); - auto energy = sqrt((fData.GetData()*conj(fData.GetData())).colwise().sum()); + auto energy = fData.GetData().abs2().colwise().sum(); - auto normalization = sqrt(2)/(energy*fNoiseStd); + auto normalization = sqrt(2)/(sqrt(energy)*fNoiseStd); KTDEBUG(ctemplatelog, "Calculating normalized matrix"); auto normalized = fData.GetData().rowwise()*normalization; KTDEBUG(ctemplatelog, "Store transposed matrix in newData"); + // Eigen does not calculate anything before this assignment newData.GetData() = normalized.transpose(); KTDEBUG(ctemplatelog, "Template matrix has shape (" From f88c8c111b78f71a05d5e61c2e09c4606b83cd2b Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Thu, 6 May 2021 13:51:20 +0200 Subject: [PATCH 38/48] Fix issue in KTConvertToTemplate --- Source/Transform/KTConvertToTemplate.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/Transform/KTConvertToTemplate.cc b/Source/Transform/KTConvertToTemplate.cc index 3571d78e4..95d9b9bdb 100644 --- a/Source/Transform/KTConvertToTemplate.cc +++ b/Source/Transform/KTConvertToTemplate.cc @@ -25,6 +25,8 @@ namespace Katydid KTConvertToTemplate::KTConvertToTemplate(const std::string& name) : KTProcessor(name), fNoiseTemperature(0.), + fBandwidth(0.), + fNoiseStd(0.), fTSSignal("template-matrix", this), fTSSlot("ts-matrix", this, &KTConvertToTemplate::Convert, &fTSSignal) { @@ -57,7 +59,7 @@ namespace Katydid KTAggregatedTemplateMatrixData& newData = fData.Of< KTAggregatedTemplateMatrixData >(); KTDEBUG(ctemplatelog, "Calculating energy"); - auto energy = fData.GetData().abs2().colwise().sum(); + auto energy = (fData.GetData()*conj(fData.GetData())).colwise().sum(); auto normalization = sqrt(2)/(sqrt(energy)*fNoiseStd); From 3e0825971c187c8fd7839b61f67ff558faefce20 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Thu, 6 May 2021 13:53:58 +0200 Subject: [PATCH 39/48] Add KTInnerProductOptimizer --- Source/Data/CMakeLists.txt | 6 +- .../SpectrumAnalysis/KTInnerProductData.cc | 15 ++ .../SpectrumAnalysis/KTInnerProductData.hh | 28 +++ .../KTInnerProductOptimizerData.cc | 15 ++ .../KTInnerProductOptimizerData.hh | 33 +++ .../KTMultiChannelInnerProductData.cc | 188 ------------------ .../KTMultiChannelInnerProductData.hh | 106 ---------- Source/SpectrumAnalysis/CMakeLists.txt | 2 + .../KTInnerProductOptimizer.cc | 57 ++++++ .../KTInnerProductOptimizer.hh | 95 +++++++++ 10 files changed, 249 insertions(+), 296 deletions(-) create mode 100644 Source/Data/SpectrumAnalysis/KTInnerProductData.cc create mode 100644 Source/Data/SpectrumAnalysis/KTInnerProductData.hh create mode 100644 Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.cc create mode 100644 Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.hh delete mode 100644 Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc delete mode 100644 Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh create mode 100644 Source/SpectrumAnalysis/KTInnerProductOptimizer.cc create mode 100644 Source/SpectrumAnalysis/KTInnerProductOptimizer.hh diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index a17160444..6bdee2392 100644 --- a/Source/Data/CMakeLists.txt +++ b/Source/Data/CMakeLists.txt @@ -23,7 +23,8 @@ set (DATA_HEADERFILES SpectrumAnalysis/KTTimeSeriesDistData.hh SpectrumAnalysis/KTWV2DData.hh SpectrumAnalysis/KTWignerVilleData.hh - # SpectrumAnalysis/KTMultiChannelInnerProductData.hh + SpectrumAnalysis/KTInnerProductData.hh + SpectrumAnalysis/KTInnerProductOptimizerData.hh EventAnalysis/KTClassifierResultsData.hh EventAnalysis/KTFrequencyCandidate.hh EventAnalysis/KTFrequencyCandidateData.hh @@ -92,7 +93,8 @@ set (DATA_SOURCEFILES SpectrumAnalysis/KTTimeSeriesDistData.cc SpectrumAnalysis/KTWV2DData.cc SpectrumAnalysis/KTWignerVilleData.cc - # SpectrumAnalysis/KTMultiChannelInnerProductData.cc + SpectrumAnalysis/KTInnerProductOptimizerData.cc + SpectrumAnalysis/KTInnerProductData.cc EventAnalysis/KTClassifierResultsData.cc EventAnalysis/KTFrequencyCandidate.cc EventAnalysis/KTFrequencyCandidateData.cc diff --git a/Source/Data/SpectrumAnalysis/KTInnerProductData.cc b/Source/Data/SpectrumAnalysis/KTInnerProductData.cc new file mode 100644 index 000000000..740358010 --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTInnerProductData.cc @@ -0,0 +1,15 @@ +/* + * KTInnerProductData.cc + * + * Created on: Apr 28, 2021 + * Author: F. Thomas + */ + +#include "KTInnerProductData.hh" + + +namespace Katydid +{ + const std::string KTInnerProductData::sName("snr-matrix"); + +} /* namespace Katydid */ diff --git a/Source/Data/SpectrumAnalysis/KTInnerProductData.hh b/Source/Data/SpectrumAnalysis/KTInnerProductData.hh new file mode 100644 index 000000000..51d413536 --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTInnerProductData.hh @@ -0,0 +1,28 @@ +/* + * KTInnerProductData.hh + * + * Created on: Apr 28, 2021 + * Author: F. Thomas + */ + +#ifndef KTINNERPRODUCTDATA_HH_ +#define KTINNERPRODUCTDATA_HH_ + +#include "KTData.hh" +#include "KTPhysicalArrayComplex.hh" + + +namespace Katydid +{ + + class KTInnerProductData: public KTPhysicalArray<2, std::complex>, public Nymph::KTExtensibleData< KTInnerProductData > + { + public: + + static const std::string sName; + + }; + +} /* namespace Katydid */ + +#endif /* KTINNERPRODUCTDATA_HH_ */ diff --git a/Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.cc b/Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.cc new file mode 100644 index 000000000..7b4a9930f --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.cc @@ -0,0 +1,15 @@ +/* + * KTInnerProductOptimizerData.cc + * + * Created on: May 5, 2021 + * Author: F. Thomas + */ + +#include "KTInnerProductOptimizerData.hh" + + +namespace Katydid +{ + const std::string KTInnerProductOptimizerData::sName("match-result"); + +} /* namespace Katydid */ diff --git a/Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.hh b/Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.hh new file mode 100644 index 000000000..457f8d00e --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTInnerProductOptimizerData.hh @@ -0,0 +1,33 @@ +/* + * KTInnerProductOptimizerData.hh + * + * Created on: May 5, 2021 + * Author: F. Thomas + */ + +#ifndef KTINNERPRODUCTOPTIMIZERDATA_HH_ +#define KTINNERPRODUCTOPTIMIZERDATA_HH_ + +#include "KTData.hh" +#include "KTPhysicalArrayComplex.hh" + +#include + + +namespace Katydid +{ + + class KTInnerProductOptimizerData: public Nymph::KTExtensibleData< KTInnerProductOptimizerData > + { + public: + + static const std::string sName; + + Eigen::ArrayXd fMaxVals; + Eigen::Array< Eigen::MatrixXf::Index, Eigen::Dynamic, 1> fMaxInds; + + }; + +} /* namespace Katydid */ + +#endif /* KTINNERPRODUCTOPTIMIZERDATA_HH_ */ diff --git a/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc b/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc deleted file mode 100644 index d16c7bf08..000000000 --- a/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.cc +++ /dev/null @@ -1,188 +0,0 @@ -/* - * KTConvolvedSpectrumData.cc - * - * Created on: Aug 25, 2017 - * Author: ezayas - */ - -#include "KTConvolvedSpectrumData.hh" - - -namespace Katydid -{ - const std::string KTConvolvedPowerSpectrumData::sName("convolved-power-spectrum"); - - KTConvolvedPowerSpectrumData::KTConvolvedPowerSpectrumData() : - KTPowerSpectrumDataCore(), - KTExtensibleData() - { - } - - KTConvolvedPowerSpectrumData::~KTConvolvedPowerSpectrumData() - { - } - - KTConvolvedPowerSpectrumData& KTConvolvedPowerSpectrumData::SetNComponents(unsigned num) - { - unsigned oldSize = fSpectra.size(); - // if num < oldSize - for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) - { - delete fSpectra[iComponent]; - } - fSpectra.resize(num); - // if num > oldSize - for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) - { - fSpectra[iComponent] = NULL; - } - return *this; - } - - const std::string KTConvolvedFrequencySpectrumDataFFTW::sName("convolved-frequency-spectrum-fftw"); - - KTConvolvedFrequencySpectrumDataFFTW::KTConvolvedFrequencySpectrumDataFFTW() : - KTFrequencySpectrumDataFFTWCore(), - KTExtensibleData() - { - } - - KTConvolvedFrequencySpectrumDataFFTW::~KTConvolvedFrequencySpectrumDataFFTW() - { - } - - KTConvolvedFrequencySpectrumDataFFTW& KTConvolvedFrequencySpectrumDataFFTW::SetNComponents(unsigned num) - { - unsigned oldSize = fSpectra.size(); - // if num < oldSize - for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) - { - delete fSpectra[iComponent]; - } - fSpectra.resize(num); - // if num > oldSize - for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) - { - fSpectra[iComponent] = NULL; - } - return *this; - } - - const std::string KTConvolvedFrequencySpectrumDataPolar::sName("convolved-frequency-spectrum-polar"); - - KTConvolvedFrequencySpectrumDataPolar::KTConvolvedFrequencySpectrumDataPolar() : - KTFrequencySpectrumDataPolarCore(), - KTExtensibleData() - { - } - - KTConvolvedFrequencySpectrumDataPolar::~KTConvolvedFrequencySpectrumDataPolar() - { - } - - KTConvolvedFrequencySpectrumDataPolar& KTConvolvedFrequencySpectrumDataPolar::SetNComponents(unsigned num) - { - unsigned oldSize = fSpectra.size(); - // if num < oldSize - for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) - { - delete fSpectra[iComponent]; - } - fSpectra.resize(num); - // if num > oldSize - for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) - { - fSpectra[iComponent] = NULL; - } - return *this; - } - - const std::string KTConvolvedPowerSpectrumVarianceData::sName("convolved-power-spectrum-variance"); - - KTConvolvedPowerSpectrumVarianceData::KTConvolvedPowerSpectrumVarianceData() : - KTFrequencySpectrumVarianceDataCore(), - KTExtensibleData() - { - } - - KTConvolvedPowerSpectrumVarianceData::~KTConvolvedPowerSpectrumVarianceData() - { - } - - KTConvolvedPowerSpectrumVarianceData& KTConvolvedPowerSpectrumVarianceData::SetNComponents(unsigned num) - { - unsigned oldSize = fSpectra.size(); - // if num < oldSize - for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) - { - delete fSpectra[iComponent]; - } - fSpectra.resize(num); - // if num > oldSize - for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) - { - fSpectra[iComponent] = NULL; - } - return *this; - } - - const std::string KTConvolvedFrequencySpectrumVarianceDataFFTW::sName("convolved-frequency-spectrum-variance-fftw"); - - KTConvolvedFrequencySpectrumVarianceDataFFTW::KTConvolvedFrequencySpectrumVarianceDataFFTW() : - KTFrequencySpectrumVarianceDataCore(), - KTExtensibleData() - { - } - - KTConvolvedFrequencySpectrumVarianceDataFFTW::~KTConvolvedFrequencySpectrumVarianceDataFFTW() - { - } - - KTConvolvedFrequencySpectrumVarianceDataFFTW& KTConvolvedFrequencySpectrumVarianceDataFFTW::SetNComponents(unsigned num) - { - unsigned oldSize = fSpectra.size(); - // if num < oldSize - for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) - { - delete fSpectra[iComponent]; - } - fSpectra.resize(num); - // if num > oldSize - for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) - { - fSpectra[iComponent] = NULL; - } - return *this; - } - - const std::string KTConvolvedFrequencySpectrumVarianceDataPolar::sName("convolved-frequency-spectrum-variance-polar"); - - KTConvolvedFrequencySpectrumVarianceDataPolar::KTConvolvedFrequencySpectrumVarianceDataPolar() : - KTFrequencySpectrumVarianceDataCore(), - KTExtensibleData() - { - } - - KTConvolvedFrequencySpectrumVarianceDataPolar::~KTConvolvedFrequencySpectrumVarianceDataPolar() - { - } - - KTConvolvedFrequencySpectrumVarianceDataPolar& KTConvolvedFrequencySpectrumVarianceDataPolar::SetNComponents(unsigned num) - { - unsigned oldSize = fSpectra.size(); - // if num < oldSize - for (unsigned iComponent = num; iComponent < oldSize; ++iComponent) - { - delete fSpectra[iComponent]; - } - fSpectra.resize(num); - // if num > oldSize - for (unsigned iComponent = oldSize; iComponent < num; ++iComponent) - { - fSpectra[iComponent] = NULL; - } - return *this; - } - - -} /* namespace Katydid */ diff --git a/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh b/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh deleted file mode 100644 index 1f9b8b170..000000000 --- a/Source/Data/SpectrumAnalysis/KTMultiChannelInnerProductData.hh +++ /dev/null @@ -1,106 +0,0 @@ -/* - * KTMultiChannelInnerProductData.hh - * - * Created on: Apr 28, 2021 - * Author: F. Thomas - */ - -#ifndef KTMULTICHANNELINNERPRODUCTDATA_HH_ -#define KTMULTICHANNELINNERPRODUCTDATA_HH_ - -#include "KTData.hh" - -#include "KTPowerSpectrum.hh" -#include "KTPowerSpectrumData.hh" -#include "KTFrequencySpectrumFFTW.hh" -#include "KTFrequencySpectrumDataFFTW.hh" -#include "KTFrequencySpectrumPolar.hh" -#include "KTFrequencySpectrumDataPolar.hh" - -#include - -namespace Katydid -{ - - class KTConvolvedPowerSpectrumData : public KTPowerSpectrumDataCore, public Nymph::KTExtensibleData< KTConvolvedPowerSpectrumData > - { - public: - KTConvolvedPowerSpectrumData(); - virtual ~KTConvolvedPowerSpectrumData(); - - KTConvolvedPowerSpectrumData& SetNComponents(unsigned channels); - - public: - static const std::string sName; - - }; - - class KTConvolvedFrequencySpectrumDataFFTW : public KTFrequencySpectrumDataFFTWCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumDataFFTW > - { - public: - KTConvolvedFrequencySpectrumDataFFTW(); - virtual ~KTConvolvedFrequencySpectrumDataFFTW(); - - KTConvolvedFrequencySpectrumDataFFTW& SetNComponents(unsigned channels); - - public: - static const std::string sName; - - }; - - class KTConvolvedFrequencySpectrumDataPolar : public KTFrequencySpectrumDataPolarCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumDataPolar > - { - public: - KTConvolvedFrequencySpectrumDataPolar(); - virtual ~KTConvolvedFrequencySpectrumDataPolar(); - - KTConvolvedFrequencySpectrumDataPolar& SetNComponents(unsigned channels); - - public: - static const std::string sName; - - }; - - class KTConvolvedPowerSpectrumVarianceData : public KTFrequencySpectrumVarianceDataCore, public Nymph::KTExtensibleData< KTConvolvedPowerSpectrumVarianceData > - { - public: - KTConvolvedPowerSpectrumVarianceData(); - virtual ~KTConvolvedPowerSpectrumVarianceData(); - - KTConvolvedPowerSpectrumVarianceData& SetNComponents(unsigned channels); - - public: - static const std::string sName; - - }; - - class KTConvolvedFrequencySpectrumVarianceDataFFTW : public KTFrequencySpectrumVarianceDataCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumVarianceDataFFTW > - { - public: - KTConvolvedFrequencySpectrumVarianceDataFFTW(); - virtual ~KTConvolvedFrequencySpectrumVarianceDataFFTW(); - - KTConvolvedFrequencySpectrumVarianceDataFFTW& SetNComponents(unsigned channels); - - public: - static const std::string sName; - - }; - - class KTConvolvedFrequencySpectrumVarianceDataPolar : public KTFrequencySpectrumVarianceDataCore, public Nymph::KTExtensibleData< KTConvolvedFrequencySpectrumVarianceDataPolar > - { - public: - KTConvolvedFrequencySpectrumVarianceDataPolar(); - virtual ~KTConvolvedFrequencySpectrumVarianceDataPolar(); - - KTConvolvedFrequencySpectrumVarianceDataPolar& SetNComponents(unsigned channels); - - public: - static const std::string sName; - - }; - - -} /* namespace Katydid */ - -#endif /* KTMULTICHANNELINNERPRODUCTDATA_HH_ */ diff --git a/Source/SpectrumAnalysis/CMakeLists.txt b/Source/SpectrumAnalysis/CMakeLists.txt index ffa291964..547e8fcde 100644 --- a/Source/SpectrumAnalysis/CMakeLists.txt +++ b/Source/SpectrumAnalysis/CMakeLists.txt @@ -24,6 +24,7 @@ set (SPECTRUMANALYSIS_HEADERFILES KTSpectrogramStriper.hh KTSwitchFFTWPolar.hh KTVariableSpectrumDiscriminator.hh + KTInnerProductOptimizer.hh #KTMultiChannelInnerProduct.hh #KTAutoCorrMatrix.hh ) @@ -51,6 +52,7 @@ set (SPECTRUMANALYSIS_SOURCEFILES KTSpectrogramStriper.cc KTSwitchFFTWPolar.cc KTVariableSpectrumDiscriminator.cc + KTInnerProductOptimizer.cc #KTMultiChannelInnerProduct.cc #KTAutoCorrMatrix.cc ) diff --git a/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc b/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc new file mode 100644 index 000000000..2b11948b4 --- /dev/null +++ b/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc @@ -0,0 +1,57 @@ +/* + * KTInnerProductOptimizer.cc + * + * Created on: May 5, 2021 + * Author: F. Thomas + */ + +#include "KTInnerProductOptimizer.hh" +#include "KTInnerProductOptimizerData.hh" +#include "KTInnerProductData.hh" + +#include "KTLogger.hh" + +#include + +#include "param.hh" + + +namespace Katydid +{ + KTLOGGER(ipolog, "KTInnerProductOptimizer"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTInnerProductOptimizer, "inner-product-optimizer"); + + KTInnerProductOptimizer::KTInnerProductOptimizer(const std::string& name) : + KTProcessor(name), + fOptSignal("opt", this), + fIPSlot("snr-matrix", this, &KTInnerProductOptimizer::FindOptimum, &fOptSignal) + { + } + + + bool KTInnerProductOptimizer::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + + return true; + } + + bool KTInnerProductOptimizer::FindOptimum(KTInnerProductData& fData) + { + KTInnerProductOptimizerData& fOpt = fData.Of(); + + auto snr = fData.GetData().abs().eval(); + + fOpt.fMaxInds.resize(snr.cols()); + fOpt.fMaxVals.resize(snr.cols()); + + for(int i=0;i + +namespace scarab +{ + class param_node; +} + +namespace Katydid +{ + + class KTInnerProductData; + + /*! + @class KTConvertToTemplate + @author F. Thomas + + @brief Converts data to Matched Filter templates + + @details + Converts a matrix of timeseries data to a matrix of normalized Matched Filter templates. + The normalization includes a noise contribution, which assumes white Gaussian noise. + Using this normalization the MF output will be a dimensionless SNR. + + For an input timeseries X the normalized template is + ``` + Y = 1/sqrt(X^H*C^-1*X) C^-1 X + ``` + where C is the noise covariance matrix. + With the white Gaussian assumption this simplifies to + ``` + Y = sqrt(2)/sqrt(X^2*V) * X + ``` + where V is the Johnson-Nyquist noise variance. + + Configuration name: "template-converter" + + Available configuration values: + - "T": double -- Noise temperature + - "bandwidth": double -- Bandwidth of the input data, i.e. half the sampling rate + + Slots: + - "ts-matrix": void (KTDataPtr) -- Converts timeseries to normalized MF template; Requires KTAggregatedTSMatrixData; Adds KTAggregatedTemplateMatrixDatar; Emits signal "template-matrix" + + Signals: + - "template-matrix": void (KTDataPtr) -- Emitted conversion; Guarantees KTAggregatedTemplateMatrixData. + + */ + class KTInnerProductOptimizer : public Nymph::KTProcessor + { + public: + KTInnerProductOptimizer(const std::string& name = "inner-product-optimizer"); + + bool Configure(const scarab::param_node* node); + + + + public: + bool FindOptimum(KTInnerProductData& fData); + + private: + + //*************** + // Signals + //*************** + + private: + Nymph::KTSignalData fOptSignal; + + //*************** + // Slots + //*************** + + private: + Nymph::KTSlotDataOneType< KTInnerProductData > fIPSlot; + + }; + +} /* namespace Katydid */ +#endif /* KTINNERPRODUCT_OPTIMIZER_HH_ */ From 8885117220a4a14ded0d387694083c53f9c39826 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 11 May 2021 11:27:50 +0200 Subject: [PATCH 40/48] Add proper inner product processor --- Source/SpectrumAnalysis/CMakeLists.txt | 4 +- Source/SpectrumAnalysis/KTInnerProduct.cc | 59 ++++++ Source/SpectrumAnalysis/KTInnerProduct.hh | 88 +++++++++ .../KTMultiChannelInnerProduct.cc | 174 ------------------ .../KTMultiChannelInnerProduct.hh | 104 ----------- 5 files changed, 149 insertions(+), 280 deletions(-) create mode 100644 Source/SpectrumAnalysis/KTInnerProduct.cc create mode 100644 Source/SpectrumAnalysis/KTInnerProduct.hh delete mode 100644 Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc delete mode 100644 Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh diff --git a/Source/SpectrumAnalysis/CMakeLists.txt b/Source/SpectrumAnalysis/CMakeLists.txt index 547e8fcde..a8c4d0495 100644 --- a/Source/SpectrumAnalysis/CMakeLists.txt +++ b/Source/SpectrumAnalysis/CMakeLists.txt @@ -25,7 +25,7 @@ set (SPECTRUMANALYSIS_HEADERFILES KTSwitchFFTWPolar.hh KTVariableSpectrumDiscriminator.hh KTInnerProductOptimizer.hh - #KTMultiChannelInnerProduct.hh + KTInnerProduct.hh #KTAutoCorrMatrix.hh ) @@ -53,7 +53,7 @@ set (SPECTRUMANALYSIS_SOURCEFILES KTSwitchFFTWPolar.cc KTVariableSpectrumDiscriminator.cc KTInnerProductOptimizer.cc - #KTMultiChannelInnerProduct.cc + KTInnerProduct.cc #KTAutoCorrMatrix.cc ) diff --git a/Source/SpectrumAnalysis/KTInnerProduct.cc b/Source/SpectrumAnalysis/KTInnerProduct.cc new file mode 100644 index 000000000..d5b16141e --- /dev/null +++ b/Source/SpectrumAnalysis/KTInnerProduct.cc @@ -0,0 +1,59 @@ +/* + * KTInnerProduct.cc + * + * Created on: Apr 28, 2021 + * Author: F. Thomas + */ + +#include "KTInnerProduct.hh" +#include "KTAggregatedTemplateMatrixData.hh" +#include "KTAggregatedTSMatrixData.hh" +#include "KTInnerProductData.hh" + +#include "KTLogger.hh" + +#include "param.hh" + + +namespace Katydid +{ + KTLOGGER(iprodlog, "KTInnerProduct"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTInnerProduct, "inner-product"); + + KTInnerProduct::KTInnerProduct(const std::string& name) : + KTProcessor(name), + fSnrSignal("snr-matrix", this), + fDataSlot("data-matrix", this, &KTInnerProduct::Multiply, &fSnrSignal), + fTemplateSlot("template-matrix", this, &KTInnerProduct::SetTemplates) //does not emit a signal on purpose + { + } + + + bool KTInnerProduct::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + return true; + } + + bool KTInnerProduct::SetTemplates(KTAggregatedTemplateMatrixData& fTemplate) + { + KTDEBUG(iprodlog, "Setting Template matrix"); + fTemplateMatrix = std::move(fTemplate.GetData()); + + return true; + } + + bool KTInnerProduct::Multiply(KTAggregatedTSMatrixData& fData) + { + KTInnerProductData& newData = fData.Of< KTInnerProductData >(); + + KTDEBUG(iprodlog, "Run inner product"); + newData.GetData() = fTemplateMatrix.matrix() * fData.GetData().matrix(); + + return true; + } + +} /* namespace Katydid */ diff --git a/Source/SpectrumAnalysis/KTInnerProduct.hh b/Source/SpectrumAnalysis/KTInnerProduct.hh new file mode 100644 index 000000000..60584d0cc --- /dev/null +++ b/Source/SpectrumAnalysis/KTInnerProduct.hh @@ -0,0 +1,88 @@ +/** + @file KTInnerProduct.hh + @brief Contains KTInnerProduct + @details Applies the Matched Filter inner product for multiple channels + @author: F. Thomas + @date: Apr 28, 2021 + */ + +#ifndef KTINNERPRODUCT_HH_ +#define KTINNERPRODUCT_HH_ + +#include "KTProcessor.hh" + +#include "KTMemberVariable.hh" +#include "KTSlot.hh" + +#include + +namespace scarab +{ + class param_node; +} + +namespace Katydid +{ + + class KTAggregatedTemplateMatrixData; + class KTAggregatedTSMatrixData; + + /*! + @class KTInnerProduct + @author F. Thomas + + @brief Applies the Matched Filter inner product for multiple channels + + @details + Frequency-domain implementation of a single-pole RC low-pass filter. + This is no-doubt a non-ideal implementation, but demonstrates the features of a processor very well. + + The relationship between the cutoff frequency, f_c, and the RC constant is: + ``` + RC = 1 / 2*pi*f_c + ``` + + Configuration name: "low-pass-filter" + + Available configuration values: + - "rc": double -- RC time constant of the filter + + Slots: + - "fs-polar": void (KTDataPtr) -- Applies a low-pass filter; Requires KTFrequencySpectrumDataPolar; Adds KTLowPassFilteredFSDataPolar; Emits signal "fs-polar" + - "ts-fftw": void (KTDataPtr) -- Applies a low-pass filter; Requires KTFrequencySpectrumDataFFTW; Adds KTLowPassFilteredFSDataFFTW; Emits signal "fs-fftw" + + Signals: + - "snr": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredPSData. + */ + class KTInnerProduct : public Nymph::KTProcessor + { + public: + KTInnerProduct(const std::string& name = "inner-product"); + bool Configure(const scarab::param_node* node); + + + public: + bool SetTemplates(KTAggregatedTemplateMatrixData& fTemplate); + bool Multiply(KTAggregatedTSMatrixData& fData); + + //*************** + // Signals + //*************** + + private: + Nymph::KTSignalData fSnrSignal; + + Eigen::ArrayXXcd fTemplateMatrix; + + //*************** + // Slots + //*************** + + private: + Nymph::KTSlotDataOneType< KTAggregatedTSMatrixData > fDataSlot; + Nymph::KTSlotDataOneType< KTAggregatedTemplateMatrixData > fTemplateSlot; + + }; + +} /* namespace Katydid */ +#endif /* KTINNERPRODUCT_HH_ */ diff --git a/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc b/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc deleted file mode 100644 index e2c3f4e3d..000000000 --- a/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.cc +++ /dev/null @@ -1,174 +0,0 @@ -/* - * KTMultiChannelInnerProduct.cc - * - * Created on: Apr 28, 2021 - * Author: F. Thomas - */ - -#include "KTLowPassFilter.hh" - -#include "KTFrequencySpectrumDataPolar.hh" -#include "KTFrequencySpectrumDataFFTW.hh" -#include "KTFrequencySpectrumFFTW.hh" -#include "KTFrequencySpectrumPolar.hh" -#include "KTLogger.hh" -#include "KTLowPassFilteredData.hh" -#include "KTMath.hh" -#include "KTPowerSpectrum.hh" -#include "KTPowerSpectrumData.hh" - -#include "param.hh" - - -namespace Katydid -{ - KTLOGGER(gclog, "KTLowPassFilter"); - - // Register the processor - KT_REGISTER_PROCESSOR(KTLowPassFilter, "low-pass-filter"); - - KTLowPassFilter::KTLowPassFilter(const std::string& name) : - KTProcessor(name), - fRC(0.), - fFSPolarSignal("fs-polar", this), - fFSFFTWSignal("fs-fftw", this), - fPSSignal("ps", this), - fFSPolarSlot("fs-polar", this, &KTLowPassFilter::Filter, &fFSPolarSignal), - fFSFFTWSlot("fs-fftw", this, &KTLowPassFilter::Filter, &fFSFFTWSignal), - fPSSlot("ps", this, &KTLowPassFilter::Filter, &fPSSignal) - { - } - - KTLowPassFilter::~KTLowPassFilter() - { - } - - bool KTLowPassFilter::Configure(const scarab::param_node* node) - { - if (node == NULL) return false; - - SetRC(node->get_value("rc", GetRC())); - - return true; - } - - bool KTLowPassFilter::Filter(KTFrequencySpectrumDataPolar& fsData) - { - unsigned nComponents = fsData.GetNComponents(); - KTLowPassFilteredFSDataPolar& newData = fsData.Of< KTLowPassFilteredFSDataPolar >().SetNComponents(nComponents); - - for (unsigned iComponent=0; iComponent().SetNComponents(nComponents); - - for (unsigned iComponent=0; iComponent().SetNComponents(nComponents); - - for (unsigned iComponent=0; iComponentsize(); - KTFrequencySpectrumPolar* newSpectrum = new KTFrequencySpectrumPolar(nBins, frequencySpectrum->GetRangeMin(), frequencySpectrum->GetRangeMax()); - newSpectrum->SetNTimeBins(frequencySpectrum->GetNTimeBins()); - - double twoPiRCf, abs, arg; - for (unsigned iBin = 0; iBin < nBins; ++iBin) - { - abs = (*frequencySpectrum)(iBin).abs(); - arg = (*frequencySpectrum)(iBin).arg(); - twoPiRCf = KTMath::TwoPi() * fRC * frequencySpectrum->GetBinCenter(iBin); - (*newSpectrum)(iBin).set_polar(abs / sqrt(1 + twoPiRCf*twoPiRCf), arg - atan(twoPiRCf)); - } - - return newSpectrum; - } - - KTFrequencySpectrumFFTW* KTLowPassFilter::Filter(const KTFrequencySpectrumFFTW* frequencySpectrum) const - { - KTDEBUG(gclog, "Creating new FS for filtered data"); - unsigned nBins = frequencySpectrum->size(); - KTFrequencySpectrumFFTW* newSpectrum = new KTFrequencySpectrumFFTW(nBins, frequencySpectrum->GetRangeMin(), frequencySpectrum->GetRangeMax()); - newSpectrum->SetNTimeBins(frequencySpectrum->GetNTimeBins()); - - double twoPiRCf, real, imag, denom, newReal, newImag; - for (unsigned iBin = 0; iBin < nBins; ++iBin) - { - real = (*frequencySpectrum)(iBin)[0]; - imag = (*frequencySpectrum)(iBin)[1]; - twoPiRCf = KTMath::TwoPi() * fRC * frequencySpectrum->GetBinCenter(iBin); - denom = 1 + twoPiRCf * twoPiRCf; - (*newSpectrum)(iBin)[0] = (real + imag * twoPiRCf) / denom; - (*newSpectrum)(iBin)[1] = (imag - real * twoPiRCf) / denom; - } - - return newSpectrum; - } - - KTPowerSpectrum* KTLowPassFilter::Filter(const KTPowerSpectrum* powerSpectrum) const - { - KTDEBUG(gclog, "Creating new PS for filtered data"); - unsigned nBins = powerSpectrum->size(); - KTPowerSpectrum* newSpectrum = new KTPowerSpectrum(nBins, powerSpectrum->GetRangeMin(), powerSpectrum->GetRangeMax()); - - double twoPiRCf; - for (unsigned iBin = 0; iBin < nBins; ++iBin) - { - twoPiRCf = KTMath::TwoPi() * fRC * powerSpectrum->GetBinCenter(iBin); - (*newSpectrum)(iBin) = (*powerSpectrum)(iBin) / (1 + twoPiRCf * twoPiRCf); - } - - return newSpectrum; - } - -} /* namespace Katydid */ diff --git a/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh b/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh deleted file mode 100644 index b37dfc1eb..000000000 --- a/Source/SpectrumAnalysis/KTMultiChannelInnerProduct.hh +++ /dev/null @@ -1,104 +0,0 @@ -/** - @file KTMultiChannelInnerProduct.hh - @brief Contains KTMultiChannelInnerProduct - @details Applies the Matched Filter inner product for multiple channels - @author: F. Thomas - @date: Apr 28, 2021 - */ - -#ifndef KTMULTICHANNELINNERPRODUCT_HH_ -#define KTMULTICHANNELINNERPRODUCT_HH_ - -#include "KTProcessor.hh" - -#include "KTMemberVariable.hh" -#include "KTSlot.hh" - -#include - -namespace scarab -{ - class param_node; -} - -namespace Katydid -{ - class KTFrequencySpectrumDataFFTW; - class KTFrequencySpectrumDataPolar; - class KTFrequencySpectrumFFTW; - class KTFrequencySpectrumPolar; - class KTPowerSpectrum; - class KTPowerSpectrumData; - - /*! - @class KTMultiChannelInnerProduct - @author F. Thomas - - @brief Applies the Matched Filter inner product for multiple channels - - @details - Frequency-domain implementation of a single-pole RC low-pass filter. - This is no-doubt a non-ideal implementation, but demonstrates the features of a processor very well. - - The relationship between the cutoff frequency, f_c, and the RC constant is: - ``` - RC = 1 / 2*pi*f_c - ``` - - Configuration name: "low-pass-filter" - - Available configuration values: - - "rc": double -- RC time constant of the filter - - Slots: - - "fs-polar": void (KTDataPtr) -- Applies a low-pass filter; Requires KTFrequencySpectrumDataPolar; Adds KTLowPassFilteredFSDataPolar; Emits signal "fs-polar" - - "fs-fftw": void (KTDataPtr) -- Applies a low-pass filter; Requires KTFrequencySpectrumDataFFTW; Adds KTLowPassFilteredFSDataFFTW; Emits signal "fs-fftw" - - "ps": void (KTDataPtr) -- Applies a low-pass filter; Requires KTPowerSpectrum; Adds KTLowPassFilteredPSData; Emits signal "ps" - - Signals: - - "fs-polar": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredFSDataPolar. - - "fs-fftw": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredFSDataFFTW. - - "ps": void (KTDataPtr) -- Emitted upon low-pass filtering; Guarantees KTLowPassFilteredPSData. - */ - class KTMultiChannelInnerProduct : public Nymph::KTProcessor - { - public: - KTMultiChannelInnerProduct(const std::string& name = "multi-channel-inner-product"); - virtual ~KTMultiChannelInnerProduct(); - - bool Configure(const scarab::param_node* node); - - MEMBERVARIABLE(double, RC); - - - public: - bool Filter(KTFrequencySpectrumDataPolar& fsData); - bool Filter(KTFrequencySpectrumDataFFTW& fsData); - bool Filter(KTPowerSpectrumData& psData); - - KTFrequencySpectrumPolar* Filter(const KTFrequencySpectrumPolar* frequencySpectrum) const; - KTFrequencySpectrumFFTW* Filter(const KTFrequencySpectrumFFTW* frequencySpectrum) const; - KTPowerSpectrum* Filter(const KTPowerSpectrum* powerSpectrum) const; - - //*************** - // Signals - //*************** - - private: - Nymph::KTSignalData fFSPolarSignal; - Nymph::KTSignalData fFSFFTWSignal; - Nymph::KTSignalData fPSSignal; - - //*************** - // Slots - //*************** - - private: - Nymph::KTSlotDataOneType< KTFrequencySpectrumDataPolar > fFSPolarSlot; - Nymph::KTSlotDataOneType< KTFrequencySpectrumDataFFTW > fFSFFTWSlot; - Nymph::KTSlotDataOneType< KTPowerSpectrumData > fPSSlot; - - }; - -} /* namespace Katydid */ -#endif /* KTMULTICHANNELINNERPRODUCT_HH_ */ From 503b8a4dbe5038607a2161111f41a648d983a61c Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Tue, 11 May 2021 11:29:11 +0200 Subject: [PATCH 41/48] Add complex conjugate to template conversion --- Source/Transform/KTConvertToTemplate.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Transform/KTConvertToTemplate.cc b/Source/Transform/KTConvertToTemplate.cc index 95d9b9bdb..b4e27b9c3 100644 --- a/Source/Transform/KTConvertToTemplate.cc +++ b/Source/Transform/KTConvertToTemplate.cc @@ -68,7 +68,7 @@ namespace Katydid KTDEBUG(ctemplatelog, "Store transposed matrix in newData"); // Eigen does not calculate anything before this assignment - newData.GetData() = normalized.transpose(); + newData.GetData() = conj(normalized.transpose()); KTDEBUG(ctemplatelog, "Template matrix has shape (" << newData.GetData().rows() From 569fd8f227a33330ca83e01cefbc91b332edebbd Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 14:12:14 +0200 Subject: [PATCH 42/48] Add config file for Matched Filter --- Examples/ConfigFiles/MatchedFilterConfig.yaml | 97 +++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100644 Examples/ConfigFiles/MatchedFilterConfig.yaml diff --git a/Examples/ConfigFiles/MatchedFilterConfig.yaml b/Examples/ConfigFiles/MatchedFilterConfig.yaml new file mode 100644 index 000000000..afa5437bf --- /dev/null +++ b/Examples/ConfigFiles/MatchedFilterConfig.yaml @@ -0,0 +1,97 @@ +processor-toolbox: + + processors: + + - type: egg-processor + name: egg1 + - type: matrix-aggregator + name: matrix1 + - type: template-converter + name: to-template + + - type: egg-processor + name: egg2 + - type: matrix-aggregator + name: matrix2 + + - type: inner-product + name: prod + - type: inner-product-optimizer + name: optimize + + connections: + + # First egg processing + - signal: "egg1:ts" + slot: "matrix1:ts-fftw" + + - signal: "matrix1:matrix" + slot: "to-template:ts-matrix" + + - signal: "to-template:template-matrix" + slot: "prod:template-matrix" + + # Second egg processing + - signal: "egg2:ts" + slot: "matrix2:ts-fftw" + + - signal: "matrix2:matrix" + slot: "prod:data-matrix" + + - signal: "prod:snr-matrix" + slot: "optimize:snr-matrix" + + + run-queue: + - egg1 + - egg2 + + +egg1: + filenames: ["locust_mc_Seed395_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed390_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed395_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed390_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed395_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg"] + egg-reader: egg3 + slice-size: 8192 + +egg2: + filenames: ["locust_mc_Seed395_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed390_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed395_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed390_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed395_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg", + "locust_mc_Seed385_LO25.8781G_Radius0.100_Pos0.009.egg"] + egg-reader: egg3 + slice-size: 8192 + +matrix1: + max-cols: 400 + +matrix2: + max-cols: 400 + +to-template: + + T: 10 + bandwidth: 200.0e6 From a0369ea6c9b2716318fc3c8c663ef16f3c64289a Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 14:13:37 +0200 Subject: [PATCH 43/48] Add better debug and progress messages --- Source/SpectrumAnalysis/KTInnerProduct.cc | 2 +- Source/SpectrumAnalysis/KTInnerProductOptimizer.cc | 2 ++ Source/Transform/KTMatrixAggregator.cc | 8 +++++--- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/Source/SpectrumAnalysis/KTInnerProduct.cc b/Source/SpectrumAnalysis/KTInnerProduct.cc index d5b16141e..6ee421e1d 100644 --- a/Source/SpectrumAnalysis/KTInnerProduct.cc +++ b/Source/SpectrumAnalysis/KTInnerProduct.cc @@ -50,7 +50,7 @@ namespace Katydid { KTInnerProductData& newData = fData.Of< KTInnerProductData >(); - KTDEBUG(iprodlog, "Run inner product"); + KTPROG(iprodlog, "Running inner product"); newData.GetData() = fTemplateMatrix.matrix() * fData.GetData().matrix(); return true; diff --git a/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc b/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc index 2b11948b4..1a0093fe9 100644 --- a/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc +++ b/Source/SpectrumAnalysis/KTInnerProductOptimizer.cc @@ -41,6 +41,7 @@ namespace Katydid bool KTInnerProductOptimizer::FindOptimum(KTInnerProductData& fData) { + KTInnerProductOptimizerData& fOpt = fData.Of(); auto snr = fData.GetData().abs().eval(); @@ -48,6 +49,7 @@ namespace Katydid fOpt.fMaxInds.resize(snr.cols()); fOpt.fMaxVals.resize(snr.cols()); + KTPROG(ipolog, "Finding best matches"); for(int i=0;iGetData(); - KTDEBUG(magglog, "Added data successfully"); + KTDEBUG(magglog, "Added component " << iComponent << " successfully"); } @@ -116,12 +116,13 @@ namespace Katydid } // Increase signal counter fSignalCount++; + KTDEBUG(magglog, "Signal count is " << fSignalCount); bool emitSignal = false; if (fSignalCount == fMaxCols) { KTDEBUG(magglog, "Matrix full"); fSignalCount = 0; - KTINFO(magglog, "Completed the matrix."); + emitSignal = true; } @@ -131,14 +132,15 @@ namespace Katydid //shrink the matrix to current signal count if end of data is reached ShrinkMatrix(); - KTINFO(magglog, "Completed the matrix."); emitSignal = true; } if(emitSignal) { + KTPROG(magglog, "Completed the matrix."); KTAggregatedTSMatrixData& aggMatrix = data->Of< KTAggregatedTSMatrixData >(); //adjust labels and KTAxis things? KTDEBUG(magglog, "Emitting the signal"); + KTDEBUG(magglog, "Final matrix has size (" << fBufferMat.rows() << "," << fBufferMat.cols() << ")"); aggMatrix.GetData() = std::move(fBufferMat); fMatrixSignal(data); KTDEBUG(magglog, "Reset buffer matrix"); From 9301f9e31a94c0aca5614049b560af0719a0fc12 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 14:20:29 +0200 Subject: [PATCH 44/48] Fix for KTEgg3Reader This fixes an issue for reading multiple egg files. After the first egg file the reader did only read a single time slice per file. --- Source/Time/KTEgg3Reader.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Time/KTEgg3Reader.cc b/Source/Time/KTEgg3Reader.cc index 08669309a..e19da1cdc 100644 --- a/Source/Time/KTEgg3Reader.cc +++ b/Source/Time/KTEgg3Reader.cc @@ -478,6 +478,7 @@ namespace Katydid return Nymph::KTDataPtr(); } inNewFile = true; + fReadState.fStatus = MonarchReadState::kContinueReading; } ++fRecordsProcessed; fReadState.fCurrentRecord = fM3Stream->GetAcqFirstRecordId() + fM3Stream->GetRecordCountInAcq(); From 8285a9e9d46235955ed72ae138392eb6deb784d1 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 15:01:26 +0200 Subject: [PATCH 45/48] Add OpenMP to build --- CMakeLists.txt | 3 ++- Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc | 4 ++-- .../KTVariableSpectrumDiscriminator.cc | 12 ++++++------ 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 52482df37..3f0b5c2f7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -212,7 +212,8 @@ else (DLIB_FOUND) endif (DLIB_FOUND) # OpenMP -#find_package (OpenMP) +find_package (OpenMP) + if (OPENMP_FOUND AND NOT Katydid_SINGLETHREADED) set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") diff --git a/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc b/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc index 461487389..d9a0f2dfd 100644 --- a/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc +++ b/Source/SpectrumAnalysis/KTSpectrumDiscriminator.cc @@ -251,7 +251,7 @@ namespace Katydid } // loop over bins, checking against the threshold -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double value = magnitude[iBin]; @@ -438,7 +438,7 @@ namespace Katydid // loop over bins, checking against the threshold -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double value = (*spectrum)(iBin); diff --git a/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc b/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc index 8a77e6834..3fbfbae8f 100644 --- a/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc +++ b/Source/SpectrumAnalysis/KTVariableSpectrumDiscriminator.cc @@ -429,7 +429,7 @@ namespace Katydid } // loop over bins, checking against the threshold -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double value = (*spectrum)(iBin).abs(); @@ -462,7 +462,7 @@ namespace Katydid //************** else if (fThresholdMode == eSigma) { -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double mean = (*splineImp)(iBin - fMinBin); @@ -533,7 +533,7 @@ namespace Katydid } // loop over bins, checking against the threshold -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double value = spectrum->GetAbs(iBin); @@ -566,7 +566,7 @@ namespace Katydid //************** else if (fThresholdMode == eSigma) { -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double mean = (*splineImp)(iBin - fMinBin); @@ -636,7 +636,7 @@ namespace Katydid } // loop over bins, checking against the threshold -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double value = (*spectrum)(iBin); @@ -669,7 +669,7 @@ namespace Katydid //************** else if (fThresholdMode == eSigma) { -#pragma omp parallel for private(value) +#pragma omp parallel for for (unsigned iBin=fMinBin; iBin<=fMaxBin; ++iBin) { double mean = (*splineImp)(iBin - fMinBin); From 9570e68f90d9258cc3d7a56b2216e2d5eb5d9081 Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 15:14:15 +0200 Subject: [PATCH 46/48] Add compiler optimizations for SIMD instructions --- CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3f0b5c2f7..71486f477 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,10 @@ # Minimum cmake verison 3.1 required for the variable CMAKE_CXX_STANDARD cmake_minimum_required (VERSION 3.1) +#probably shouldn't go here like that but I don't know where it has to go +#feel free to fix it +add_compile_options(-march=native -mfma) + # Define the project cmake_policy( SET CMP0048 NEW ) # version in project() From 13c6bb89832a0f1b6389e63c6825b0bfecc534ee Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 21:43:12 +0200 Subject: [PATCH 47/48] Fix for KTPhysicalArrayComplex Makes functions public that were accidentally protected --- Source/Utility/KTPhysicalArrayComplex.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Utility/KTPhysicalArrayComplex.hh b/Source/Utility/KTPhysicalArrayComplex.hh index 6cf09692c..fca488e42 100644 --- a/Source/Utility/KTPhysicalArrayComplex.hh +++ b/Source/Utility/KTPhysicalArrayComplex.hh @@ -243,11 +243,11 @@ namespace Katydid protected: matrix_type fData; std::string fLabel; + + public: size_t cols() const; size_t rows() const; - - public: const value_type& operator()(unsigned i, unsigned j) const; value_type& operator()(unsigned i, unsigned j); From 62b870d5d308adfb914608b692fe9635f86dab2f Mon Sep 17 00:00:00 2001 From: Florian Thomas Date: Wed, 12 May 2021 21:44:58 +0200 Subject: [PATCH 48/48] Update MatchedFilter test --- .../Validation/TestMatchedFilter.cc | 42 ++++++++++++++++--- 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/Source/Executables/Validation/TestMatchedFilter.cc b/Source/Executables/Validation/TestMatchedFilter.cc index a6c0a4239..d661f78b4 100644 --- a/Source/Executables/Validation/TestMatchedFilter.cc +++ b/Source/Executables/Validation/TestMatchedFilter.cc @@ -13,12 +13,16 @@ #include "KTMatrixAggregator.hh" #include "KTConvertToTemplate.hh" +#include "KTInnerProduct.hh" +#include "KTInnerProductOptimizer.hh" #include "KTTimeSeriesData.hh" #include "KTTimeSeriesFFTW.hh" #include "KTAggregatedTSMatrixData.hh" #include "KTAggregatedTemplateMatrixData.hh" +#include "KTInnerProductData.hh" +#include "KTInnerProductOptimizerData.hh" #include "KTLogger.hh" @@ -32,9 +36,11 @@ int main() // Create and setup processors KTMatrixAggregator tAgg {}; KTConvertToTemplate tConvert {}; + KTInnerProduct tInnerProd {}; + KTInnerProductOptimizer tOpt {}; - unsigned maxCols = 3; - unsigned nCols = 5; + unsigned maxCols = 5; + unsigned nCols = 19; unsigned nChannels = 3; unsigned nTimeBins = 4; @@ -80,14 +86,38 @@ int main() // Test conversion to template tConvert.Convert(dataMatrix); - KTAggregatedTemplateMatrixData& templateMatrix = data->Of< KTAggregatedTemplateMatrixData >(); - KTDEBUG(testlog, "Template Matrix: " << templateMatrix); - KTDEBUG(testlog, "Norm" << sqrt((dataMatrix.GetData()*conj(dataMatrix.GetData())).colwise().sum())); + tInnerProd.SetTemplates(templateMatrix); + tInnerProd.Multiply(dataMatrix); + + KTInnerProductData& product = data->Of< KTInnerProductData >(); + KTDEBUG(testlog, "Product matrix: " << product); + + Eigen::ArrayXd norm = dataMatrix.GetData().matrix().colwise().norm(); + + tOpt.FindOptimum(product); + + KTInnerProductOptimizerData& res = data->Of< KTInnerProductOptimizerData >(); + + KTDEBUG(testlog, "Arg: " << res.fMaxInds); + KTDEBUG(testlog, "Vals: " << res.fMaxVals); + KTDEBUG(testlog, "Norm: " << norm); + + //checking results + Eigen::ArrayXi indices = Eigen::ArrayXi::LinSpaced(dataMatrix.cols(), 0, dataMatrix.cols()-1); + + if(indices.matrix() == res.fMaxInds.matrix().cast() && + norm.isApprox(res.fMaxVals, 1e-7)) + { + KTPROG(testlog, "Result matches!"); + } else { + KTPROG(testlog, "Result does not match!"); + return 1; + } - KTDEBUG(testlog, " Product: " << (templateMatrix.GetData().transpose()*conj(dataMatrix.GetData())).colwise().sum()); + KTPROG(testlog, "Test passed!"); return 0; }