Commit 2dc7add6 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

croposp [FIX][!!!]: compute magntitude squared coherence

The magnitude squared coherence appears to be the more common quantity in
literature.
parent 8b8082b6
......@@ -28,6 +28,10 @@
* REVISIONS and CHANGES
* - 25/12/2018 V1.0 Thomas Forbriger
* - 14/02/2019 V1.1 first operational version (coarsely tested)
* - 03/03/2019 V1.2 compute magnitude squared coherency rather than
* its suqare root
* commonly in literature coherency refers to the
* magnitude squared coherency
*
* ============================================================================
*/
......@@ -35,7 +39,7 @@
* version string is set in croposp.h
*
#define CROPOSP_VERSION \
"CROPOSP V1.1 Cross power spectral density"
"CROPOSP V1.2 Cross power spectral density"
*/
#include "croposp.h"
......@@ -148,7 +152,7 @@ int main(int iargc, char* argv[])
"-phase f compute phase transfer function for pairs of\n"
" signals and write to file \"f\"\n"
" Sleeman et al (2006, eq. 13)\n"
"-coherency f compute coherecy for pairs of signals\n"
"-coherency f compute magnitude squared coherecy for pairs of signals\n"
" and write to file \"f\"\n"
"\n"
"keys to be appended to file names if required:\n"
......@@ -206,7 +210,7 @@ int main(int iargc, char* argv[])
// 21: set phase label pattern
{"patphase",arg_yes,"PHA %1-%2 (ref: %3)"},
// 22: set coherency label pattern
{"patcoherency",arg_yes,"COH %1 %2"},
{"patcoherency",arg_yes,"MSC %1 %2"},
// 23: average upon output
{"avgout",arg_no,"-"},
// 24: average upon output
......@@ -619,8 +623,8 @@ int main(int iargc, char* argv[])
// opt.compute_coherency)
/* ---------------------------------------------------------------------- */
/* compute coherency
* -----------------
/* compute magnitude squared coherency
* -----------------------------------
*/
if (opt.compute_coherency)
......@@ -639,25 +643,14 @@ int main(int iargc, char* argv[])
TFXX_assert(CPSD_vector.size()==pairs.pairs(),
"inconsistency; report this as a bug");
croposp::TNamedSeriesVector sqrtpsd_vector(pairs.size());
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);
psd::TDseries::Tcoc mscpsd=psd::abssqr(CPSD_vector[pairs(k,l)].series);
coherency_vector[pairs(k,l)].series
=acpsd/(sqrtpsd_vector[k].series*sqrtpsd_vector[l].series);
=mscpsd/(PSD_vector[k].series*PSD_vector[l].series);
coherency_vector[pairs(k,l)].label=opt.pattern_coherency;
coherency_vector[pairs(k,l)].label
=croposp::patsubst(coherency_vector[pairs(k,l)].label, "%1",
......
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