Commit 14d86b6a authored by thomas.forbriger's avatar thomas.forbriger
Browse files

libpsdxx [FIX][!!!][API]: compute magnitude squared coherency

The common measure of coherency in the literature appears to be magnitude
squared coherency. For this reason the library now implements this definition.
The previous scaling is provided as complex values of normalized cross-power
spectral density.
parent 34375a21
......@@ -74,13 +74,15 @@ $(call CHECKVARS,TF_BROWSER TF_WWWBASEDIR)
LIBHEADERS=debug.h error.h function_template_log_sampling.h psd.h helper.h
LIBCCSRC= \
dpsdcomputer_coherency.cc \
dpsdcomputer_mscoherency.cc \
dpsdcomputer_normalized_cpsd.cc \
dpsdcomputer_cross_psd.cc \
dpsdcomputer_parameters.cc \
dpsdcomputer_processor.cc \
dpsdcomputer_psd.cc \
error.cc \
function_abs.cc \
function_abssqr.cc \
function_arg.cc \
function_conj.cc \
function_lin_frequency.cc \
......
/*! \file dpsdcomputer_coherency.cc
* \brief compute coherency (implementation)
/*! \file dpsdcomputer_mscoherency.cc
* \brief compute magnitude squared coherency (implementation)
*
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 03/01/2019
*
* compute coherency (implementation)
* compute magnitude squared coherency (implementation)
*
* Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach)
*
......@@ -30,8 +30,8 @@
*
* ============================================================================
*/
#define PSDXX_DPSDCOMPUTER_COHERENCY_CC_VERSION \
"PSDXX_DPSDCOMPUTER_COHERENCY_CC V1.0"
#define PSDXX_DPSDCOMPUTER_MSCOHERENCY_CC_VERSION \
"PSDXX_DPSDCOMPUTER_MSCOHERENCY_CC V1.0"
#include <psdxx/psd.h>
#include <aff/iterator.h>
......@@ -39,27 +39,27 @@
namespace psd {
TDISeries DPSDComputer::coherency(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const
TDISeries DPSDComputer::mscoherency(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const
{
SpectralValues sv=this->processor(s1, s2, true, true);
TDISeries retval=psd::abs(sv.cpsd);
TDISeries retval=psd::abssqr(sv.cpsd);
aff::Iterator<TDISeries::Tseries> S(retval.data);
double psdmax1=aff::func::max(sv.psd1.data);
double psdmax2=aff::func::max(sv.psd2.data);
double psdlim=std::sqrt(psdmax1*psdmax2)*1.e-10;
double psdlim=std::sqrt(psdmax1*psdmax2)*1.e-20;
aff::Browser<TDISeries::Tseries> P1(sv.psd1.data);
aff::Browser<TDISeries::Tseries> P2(sv.psd2.data);
while(S.valid() && P1.valid() && P2.valid())
{
double psd=std::sqrt(*P1)*std::sqrt(*P2);
double psd=std::abs((*P1)*(*P2));
psd = psd > psdlim? psd : psdlim;
*S /= psd;
++S; ++P1; ++P2;
}
return(retval);
} // TDISeries DPSDComputer::coherency(...) const
} // TDISeries DPSDComputer::mscoherency(...) const
} // namespace psd
/* ----- END OF dpsdcomputer_coherency.cc ----- */
/* ----- END OF dpsdcomputer_mscoherency.cc ----- */
/*! \file dpsdcomputer_normalized_cpsd.cc
* \brief compute normalized cross-power spectral density (implementation)
*
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 03/03/2019
*
* compute normalized cross-power spectral density (implementation)
*
* Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
* ----
*
* REVISIONS and CHANGES
* - 03/03/2019 V1.0 Thomas Forbriger
*
* ============================================================================
*/
#define PSDXX_DPSDCOMPUTER_NORMALIZED_CPSD_CC_VERSION \
"PSDXX_DPSDCOMPUTER_NORMALIZED_CPSD_CC V1.0"
#include <psdxx/psd.h>
#include <aff/iterator.h>
#include <aff/functions/max.h>
namespace psd {
TDCISeries DPSDComputer::normalized_cpsd(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const
{
SpectralValues sv=this->processor(s1, s2, true, true);
TDCISeries retval=sv.cpsd;
aff::Iterator<TDCISeries::Tseries> S(retval.data);
double psdmax1=aff::func::max(sv.psd1.data);
double psdmax2=aff::func::max(sv.psd2.data);
double psdlim=std::sqrt(psdmax1*psdmax2)*1.e-10;
aff::Browser<TDISeries::Tseries> P1(sv.psd1.data);
aff::Browser<TDISeries::Tseries> P2(sv.psd2.data);
while(S.valid() && P1.valid() && P2.valid())
{
double psd=std::sqrt((*P1)*(*P2));
psd = psd > psdlim? psd : psdlim;
*S /= psd;
++S; ++P1; ++P2;
}
return(retval);
} // TDCISeries DPSDComputer::normalized_cpsd(...) const
} // namespace psd
/* ----- END OF dpsdcomputer_normalized_cpsd.cc ----- */
/*! \file function_abssqr.cc
* \brief square of magnitude (implementation)
*
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 03/03/2019
*
* square of magnitude (implementation)
*
* Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
* ----
*
* REVISIONS and CHANGES
* - 03/03/2019 V1.0 Thomas Forbriger
*
* ============================================================================
*/
#define PSDXX_FUNCTION_ABSSQR_CC_VERSION \
"PSDXX_FUNCTION_ABSSQR_CC V1.0"
#include <psdxx/psd.h>
#include <aff/iterator.h>
namespace psd {
TDseries abssqr(const TDCseries::Tcoc& s)
{
TDseries retval(s.shape());
aff::Iterator<TDseries> S(retval);
aff::Browser<TDCseries::Tcoc> C(s);
while(S.valid() && C.valid())
{
*S = C->real()*C->real()+C->imag()*C->imag();
++S; ++C;
}
return(retval);
} // TDseries abssqr(const TDCseries::Tcoc& s)
/* ---------------------------------------------------------------------- */
TDISeries abssqr(const TDCISeries::Tcoc& s)
{
TDISeries retval;
retval.interval=s.interval;
retval.data=abs(s.data);
return(retval);
} // TDISeries abssqr(const TDCISeries::Tcoc& s)
} // namespace psd
/* ----- END OF function_abssqr.cc ----- */
......@@ -137,6 +137,9 @@ namespace psd {
/* ====================================================================== */
//! return magnitude squared of complex values
TDISeries abssqr(const TDCISeries::Tcoc& s);
//! return absolute value (magnitude) of complex values
TDISeries abs(const TDCISeries::Tcoc& s);
......@@ -149,6 +152,9 @@ namespace psd {
//! return sqrt values
TDISeries sqrt(const TDISeries::Tcoc& s);
//! return magnitude squared of complex values
TDseries abssqr(const TDCseries::Tcoc& s);
//! return absolute value (magnitude) of complex values
TDseries abs(const TDCseries::Tcoc& s);
......@@ -201,9 +207,12 @@ namespace psd {
//! compute cross power spectrum
TDCISeries cross_psd(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const;
//! compute coherency
TDISeries coherency(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const;
//! compute normalized cross-power spectral density
TDCISeries normalized_cpsd(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const;
//! compute magnitude squared coherency
TDISeries mscoherency(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const;
private:
//! number of segments split input time series
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment