Commit 326ccad7 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

croposp [WP][FEATURE]: compute coherency

parent b51d7ca5
......@@ -41,6 +41,7 @@
#include <iostream>
#include <tsioxx/cmdlinefiles.h>
#include <tfxx/misc.h>
#include <aff/seriesoperators.h>
using std::cout;
using std::cerr;
......@@ -525,6 +526,55 @@ int main(int iargc, char* argv[])
} // if (opt.compute_transfer || opt.compute_phase || opt.compute_npsd ||
// opt.compute_coherency)
/* ---------------------------------------------------------------------- */
/* compute coherency
* -----------------
*/
if (opt.compute_coherency)
{
croposp::TNamedSeriesVector coherency_vector;
croposp::Pairs pairs(collection_of_series.size());
TFXX_assert(PSD_vector.size()==pairs.size(),
"inconsistency; report this as a bug");
TFXX_assert(CPSD_vector.size()==pairs.pairs(),
"inconsistency; report this as a bug");
croposp::TNamedSeriesVector sqrtpsd_vector;
for (unsigned int k=0; k<PSD_vector.size(); ++k)
{
sqrtpsd_vector[k]=PSD_vector[k];
psd::TDseries::Tcoc psd=PSD_vector[k].series;
sqrtpsd_vector[k].series=psd::sqrt(psd);
} // for (unsigned int k=0; k<PSD_vector.size(); ++k)
for (unsigned int k=1; k<pairs.size(); ++k)
{
for (unsigned int l=0; l<k; ++l)
{
psd::TDseries::Tcoc acpsd=psd::abs(CPSD_vector[pairs(k,l)].series);
coherency_vector[pairs(k,l)].series
=acpsd/(sqrtpsd_vector[k].series*sqrtpsd_vector[l].series);
coherency_vector[pairs(k,l)].label=
opt.prefix_coherency+" "+CPSD_vector[pairs(k,l)].label;
} // for (unsigned int l=0; l<k; ++l)
} // for (unsigned int l=0; l<k; ++l)
write_named_series(opt.outfile_coherency,
"coherency of pairs of signals",
frequencies,
coherency_vector,
opt.verbose,
opt.overwrite);
} // if (opt.compute_coherency)
} // main()
/* ----- END OF croposp.cc ----- */
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