croposp [WP][FEATURE]: compute PSD of incoherent component

src/ts/croposp/croposp.cc

@@ -629,8 +629,8 @@ int main(int iargc, char* argv[])
   } // if (opt.compute_coherency)
 
   /* ---------------------------------------------------------------------- */
-  /* compute amplitude transfer function
-   * -----------------------------------
+  /* compute amplitude phase function
+   * --------------------------------
    */
 
   if (opt.compute_phase)
@@ -738,6 +738,92 @@ int main(int iargc, char* argv[])
 
   } // if (opt.compute_transfer)
 
+  /* ---------------------------------------------------------------------- */
+  /* compute incoherent component
+   * ----------------------------
+   */
+
+  if (opt.compute_npsd)
+  {
+    if (opt.verbose)
+    {
+      cout << endl 
+        << "Compute incoherent component." 
+        << endl;
+    }
+
+    TFXX_assert(pairs.size()>2, "requires at least three input time series");
+
+    croposp::Triples triples(pairs, false);
+    croposp::TNamedSeriesVector npsd_vector(triples.triples());
+
+    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");
+
+    for (unsigned int i=0; i<pairs.size(); ++i)
+    {
+      for (unsigned int k=1; k<pairs.size(); ++k)
+      {
+        for (unsigned int l=0; l<k; ++l)
+        {
+          if ((k!=i) && (l!=i))
+          {
+            psd::TDCseries::Tcoc Cki,Cil,Ckl;
+
+            if (pairs.swap(k,i))
+            { Cki=psd::conj(CPSD_vector[pairs(k,i)].series); }
+            else
+            { Cki=CPSD_vector[pairs(k,i)].series; }
+
+            if (pairs.swap(i,l))
+            { Cil=psd::conj(CPSD_vector[pairs(i,l)].series); }
+            else
+            { Cil=CPSD_vector[pairs(i,l)].series; }
+
+            if (pairs.swap(k,l))
+            { Ckl=psd::conj(CPSD_vector[pairs(k,l)].series); }
+            else
+            { Ckl=CPSD_vector[pairs(k,l)].series; }
+
+            psd::TDCseries::Tcoc coherent_psd=(Cki*Cil/Ckl);
+
+            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);
+            aff::Browser<psd::TDseries::Tcoc> iPSD(psd);
+
+            while(iNPSD.valid() && iCPSD.valid() && iPSD.valid())
+            {
+              *iNPSD = std::abs(*iPSD - (*iCPSD));
+              ++iNPSD;
+              ++iCPSD;
+              ++iPSD;
+            }
+
+            npsd_vector[triples(i,k,l)].series=npsd;
+            npsd_vector[triples(i,k,l)].label
+              =opt.prefix_npsd+" "
+              +vector_of_metadata[i].label+" (ref: "
+              +vector_of_metadata[k].label
+              +" "+vector_of_metadata[l].label+")";
+          }
+        }
+      }
+    }
+
+    write_named_series(opt.outfile_npsd,
+                       "power spectral density of incoherent component",
+                       frequencies,
+                       npsd_vector,
+                       opt.verbose,
+                       opt.overwrite);
+
+  } // if (opt.compute_npsd)
+
 } // main()
 
 /* ----- END OF croposp.cc ----- */