croposp [FEATURE]: provide option to average complex values

src/ts/croposp/croposp.cc

@@ -76,7 +76,7 @@ int main(int iargc, char* argv[])
     CROPOSP_VERSION "\n"
     "usage: croposp [-verbose] [-itype f] [-trim] [-datetolerance t]" "\n"
     "               [-pattern p] [-overwrite] [-DEBUG]" "\n"
-    "               [-log n] [-avgout]\n"
+    "               [-log n] [-avgout] [-avgcomplex]\n"
     "               [-nsegments n] [-divisor n] [-padfactor n] [-overlap f]\n"
     "               [-psd f] [-transfer f] [-phase f] [-coherence f]" "\n"
     "               [-npsd f] [-npsdmode f]\n"
@@ -133,6 +133,10 @@ int main(int iargc, char* argv[])
     "             phase response of input channels strongly varies with\n"
     "             frequency, results may differ; using this option\n"
     "             increases computation time\n"
+    "-avgcomplex  only effective together with -avgout\n"
+    "             if computing coherence results, average the complex\n"
+    "             fraction of cross-power spectral density values\n"
+    "             along frequency, before taking the absolute value, etc\n"
     "-nsegments n set number of segments to be used for each input signal\n"
     "-divisor n   set number of samples to an integer multiple of \"n\"\n"
     "-padfactor n pad with zeroes to increase length of series by a factor\n"
@@ -215,6 +219,8 @@ int main(int iargc, char* argv[])
     {"avgout",arg_no,"-"},
     // 24: average upon output
     {"npsdmode",arg_yes,"0"},
+    // 25: average complex value of fraction of cross-power spectral density
+    {"avgcomplex",arg_no,"-"},
     {NULL}
   };
 
@@ -305,6 +311,7 @@ int main(int iargc, char* argv[])
 
   opt.avg_on_output=cmdline.optset(23);
   opt.npsd_mode=cmdline.int_arg(24);
+  opt.avg_complex=cmdline.optset(25);
 
   TFXX_assert((opt.npsd_mode>=0) && (opt.npsd_mode<4),
               "npsd mode out of range");
@@ -321,6 +328,11 @@ int main(int iargc, char* argv[])
               "overlap fraction is out of accepted range");
   TFXX_assert(((opt.datetolerance >= 0.) && (opt.datetolerance <=1.)),
               "datetolerance is out of accepted range");
+  if (opt.avg_complex)
+  {
+    TFXX_assert(opt.avg_on_output,
+                "use -avgcomplex only together with -avgout");
+  }
   
   // extract commandline arguments
   TFXX_assert(cmdline.extra(), "missing input file");
@@ -815,8 +827,21 @@ int main(int iargc, char* argv[])
     TFXX_assert(CPSD_vector.size()==pairs.pairs(),
                 "inconsistency; report this as a bug");
 
+    // compute logarithmic frequency sampling, if required
+    psd::TDseries logfreq=frequencies;
+    if (opt.avg_complex) 
+    {
+      logfreq=psd::log_frequency(frequencies, opt.n_per_decade);
+    }
+
     for (unsigned int i=0; i<pairs.size(); ++i)
     {
+      psd::TDseries::Tcoc psd=PSD_vector[i].series.copyout();
+      if (opt.avg_complex)
+      {
+        psd=psd::log_sampling(psd, frequencies, logfreq);
+      }
+
       for (unsigned int k=1; k<pairs.size(); ++k)
       {
         for (unsigned int l=0; l<k; ++l)
@@ -840,10 +865,14 @@ int main(int iargc, char* argv[])
             else
             { Ckl=CPSD_vector[pairs(k,l)].series; }
 
-            psd::TDCseries::Tcoc coherent_psd=(Cki*Cil/Ckl);
+            psd::TDCseries coherent_psd=(Cki*Cil/Ckl);
+            if (opt.avg_complex)
+            {
+              coherent_psd=::psd::log_sampling(coherent_psd,
+                                             frequencies, logfreq);
+            }
 
             psd::TDseries npsd(coherent_psd.shape());
-            psd::TDseries::Tcoc psd(PSD_vector[i].series);
 
             aff::Iterator<psd::TDseries> iNPSD(npsd);
             aff::Browser<psd::TDCseries::Tcoc> iCPSD(coherent_psd);
@@ -891,9 +920,20 @@ int main(int iargc, char* argv[])
       }
     }
 
-    write_named_series(opt.outfile_npsd,
-                       "power spectral density of incoherent component",
-                       frequencies, npsd_vector, opt);
+    if (opt.avg_complex)
+    {
+      croposp::Options outopt=opt;
+      outopt.avg_on_output=false;
+      write_named_series(opt.outfile_npsd,
+                         "power spectral density of incoherent component",
+                         logfreq, npsd_vector, outopt);
+    }
+    else
+    {
+      write_named_series(opt.outfile_npsd,
+                         "power spectral density of incoherent component",
+                         frequencies, npsd_vector, opt);
+    }
 
   } // if (opt.compute_npsd)
 

src/ts/croposp/croposp.h

@@ -28,17 +28,18 @@
  * REVISIONS and CHANGES 
  *  - 06/02/2019   V1.0   Thomas Forbriger
  *  - 14/02/2019   V1.1   set program version 1.1
+ *  - 08/04/2019   V1.2   introduce option avgcomplex
  * 
  * ============================================================================
  */
 #define CROPOSP_VERSION \
-  "CROPOSP   V1.2   Cross power spectral density"
+  "CROPOSP   V1.3   Cross power spectral density"
 
 // include guard
 #ifndef CROPOSP_H_VERSION
 
 #define CROPOSP_H_VERSION \
-  "CROPOSP_H   V1.1"
+  "CROPOSP_H   V1.2"
 
 #include <string>
 #include <vector>
@@ -55,7 +56,7 @@ namespace croposp {
     bool verbose, trim, debug, overwrite;
     std::string inputformat, labelpattern;
     double datetolerance;
-    bool logscale, avg_on_output;
+    bool logscale, avg_on_output, avg_complex;
     unsigned int n_per_decade;
     unsigned int nsegments, divisor, padfactor;
     double overlap;