croposp [WP][FEATURE]: compute cross-psd

src/ts/croposp/croposp.cc

@@ -396,6 +396,13 @@ int main(int iargc, char* argv[])
   /* compute power spectral density
    * ------------------------------
    */
+  psd::DPSDComputer psd_computer;
+  psd_computer.set_debug(opt.debug);
+  psd_computer.set_verbose(opt.verbose);
+  psd_computer.set_nsegments(opt.nsegments);
+  psd_computer.set_divisor(opt.divisor);
+  psd_computer.set_padfactor(opt.padfactor);
+  psd_computer.set_overlap(opt.overlap);
 
   // provide container for PSD values
   croposp::TNamedSeriesVector PSD_vector;
@@ -403,14 +410,6 @@ int main(int iargc, char* argv[])
   if (opt.compute_psd || opt.compute_npsd ||
       opt.compute_coherency)
   {
-    psd::DPSDComputer psd_computer;
-    psd_computer.set_debug(opt.debug);
-    psd_computer.set_verbose(opt.verbose);
-    psd_computer.set_nsegments(opt.nsegments);
-    psd_computer.set_divisor(opt.divisor);
-    psd_computer.set_padfactor(opt.padfactor);
-    psd_computer.set_overlap(opt.overlap);
-
     croposp::TCollection::const_iterator i_series=collection_of_series.begin();
     croposp::TMetaVector::const_iterator i_meta=vector_of_metadata.begin();
     while (i_series != collection_of_series.end() &&
@@ -464,14 +463,68 @@ int main(int iargc, char* argv[])
                          opt.verbose,
                          opt.overwrite);
     } // if (opt.compute_psd)
-  } // if (opt.compute_psd || opt.compute_npsd || opt.compute_transfer ||
-    //  //      opt.compute_coherency)
+  } // if (opt.compute_psd || opt.compute_npsd || 
+    //     opt.compute_coherency)
 
   /* ---------------------------------------------------------------------- */
 
   // provide container for CPSD values
   croposp::TNamedCPSDVector CPSD_vector;
 
+  if (opt.compute_transfer || opt.compute_phase || opt.compute_npsd ||
+      opt.compute_coherency)
+  {
+    TFXX_assert(collection_of_series.size()>1,
+           "requested computation requires at least two input time series");
+
+    croposp::Pairs pairs(collection_of_series.size());
+
+    for (unsigned int k=1; k<collection_of_series.size(); ++k)
+    {
+      for (unsigned int l=0; l<k; ++l)
+      {
+        psd::TDISeries::Tcoc isk, isl;
+        isk.data=collection_of_series[k];
+        isk.interval=collection_of_series[k].header.dt;
+        isl.data=collection_of_series[l];
+        isl.interval=collection_of_series[l].header.dt;
+
+        psd::TDCISeries cpsd=psd_computer.cross_psd(isk, isl);
+
+        // compute array of frequency values, if not yet done
+        if (frequencies.size()<1)
+        {
+          if (opt.logscale)
+          {
+            frequencies=psd::log_frequency(cpsd.interval, cpsd.data.size(),
+                                           opt.n_per_decade);
+          }
+          else
+          {
+            frequencies=psd::lin_frequency(cpsd.interval, cpsd.data.size());
+          } // if (opt.logscale) ... else
+        } // if (frequencies.size()<1)
+
+        croposp::NamedCPSD named_cpsd;
+        if (opt.logscale)
+        {
+          named_cpsd.series=psd::log_sampling(cpsd, frequencies);
+        } // if (opt.logscale)
+        else
+        {
+          named_cpsd.series=cpsd.data;
+        } // if (opt.logscale) ... else
+
+        named_cpsd.label=vector_of_metadata[k].label+" "
+          +vector_of_metadata[l].label;
+        CPSD_vector[pairs(k,l)]=named_cpsd;
+
+      } // for (unsigned int l=0; l<k; ++l)
+    } // for (unsigned int l=0; l<k; ++l)
+
+  } // if (opt.compute_transfer || opt.compute_phase || opt.compute_npsd ||
+    //     opt.compute_coherency)
+
 } // main()
 
 /* ----- END OF croposp.cc ----- */