croposp [FIX][WP]: fix initialization of vectors

src/ts/croposp/Makefile

@@ -73,7 +73,7 @@ LDFLAGS=$(addprefix -L,$(LOCLIBDIR) $(subst :, ,$(SERVERLIBDIR))) $(STATIC)
 clean: ; 
 	-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
 	-find . -name \*.d | xargs --no-run-if-empty /bin/rm -v
-	-/bin/rm -vf flist *.o *.xxx *.ps *~ $(PROGRAMS)
+	-/bin/rm -vf flist *.o *.xxx *.ps *.pdf *~ $(PROGRAMS)
 
 TESTCASEMAKE=$(filter-out %.bak,$(wildcard testcases/Makefile*))
 
@@ -229,18 +229,24 @@ CPSD_OVERLAP=0.5
 CPSD_DIVISOR=100
 CPSD_PADFACTOR=5
 CPSD_LOG=--log 10
-psd.xxx: cpsd1.bin cpsd2.bin cpsd3.bin
+coherency.xxx psd.xxx: cpsd1.bin cpsd2.bin cpsd3.bin
 	croposp --DEBUG --verbose --overwrite --itype bin --trim \
 	  $(CPSD_LOG) \
 	  --nsegments=$(CPSD_NSEGMENTS) \
 	  --divisor=$(CPSD_DIVISOR) \
 	  --padfactor=$(CPSD_PADFACTOR) \
 	  --overlap=$(CPSD_OVERLAP) \
-	  --datetolerance 0.1 -psd $@ $^
+	  --datetolerance 0.1 \
+	  -psd psd.xxx \
+	  -coherency coherency.xxx \
+	  $^
 
 psd.pdf: psd.xxx
 	$(CROPOSPLOT) -v --grid -o $@ $<
 
+coherency.pdf: coherency.xxx
+	$(CROPOSPLOT) -v --grid --nology -o $@ $<
+
 %.pdp: %.pdf; evince $<; /bin/rm -fv $<
 #
 #======================================================================

src/ts/croposp/croposp.cc

@@ -387,7 +387,6 @@ int main(int iargc, char* argv[])
   TFXX_assert(!opt.compute_npsd, "-npsd is not yet implemented");
   TFXX_assert(!opt.compute_transfer, "-transfer is not yet implemented");
   TFXX_assert(!opt.compute_phase, "-phase is not yet implemented");
-  TFXX_assert(!opt.compute_coherency, "-coherency is not yet implemented");
 
   // provide container for frequency values (will be filled upon computation
   // of PSD values)
@@ -411,6 +410,16 @@ int main(int iargc, char* argv[])
   if (opt.compute_psd || opt.compute_npsd ||
       opt.compute_coherency)
   {
+    if (opt.verbose)
+    {
+      cout << endl 
+        << "Compute power spectral density for all signals." << endl;
+      if (opt.logscale)
+      {
+        cout << "Map to logarithmic frequency scale." << endl;
+      }
+    }
+
     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() &&
@@ -451,6 +460,11 @@ int main(int iargc, char* argv[])
       named_psd.label=opt.prefix_psd+" "+i_meta->label;
       PSD_vector.push_back(named_psd);
 
+      if (opt.verbose)
+      {
+        cout << "  computed: " << named_psd.label << endl;
+      }
+
       ++i_series;
       ++i_meta;
     } // while (i_series != collection_of_series.end() &&
@@ -469,21 +483,35 @@ int main(int iargc, char* argv[])
 
   /* ---------------------------------------------------------------------- */
 
+  croposp::Pairs pairs(collection_of_series.size());
+
   // provide container for CPSD values
-  croposp::TNamedCPSDVector CPSD_vector;
+  croposp::TNamedCPSDVector CPSD_vector(pairs.pairs());
 
   if (opt.compute_transfer || opt.compute_phase || opt.compute_npsd ||
       opt.compute_coherency)
   {
+    if (opt.verbose)
+    {
+      cout << endl 
+        << "Compute cross power spectral density for all pairs of signals." 
+        << endl;
+      if (opt.logscale)
+      {
+        cout << "Map to logarithmic frequency scale." << endl;
+      }
+    }
+
     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)
       {
+        TFXX_debug(opt.debug, "croposp",
+                   "next pair " << TFXX_value(k) << " " << TFXX_value(l));
+
         psd::TDISeries::Tcoc isk, isl;
         isk.data=collection_of_series[k];
         isk.interval=collection_of_series[k].header.dt;
@@ -492,6 +520,9 @@ int main(int iargc, char* argv[])
 
         psd::TDCISeries cpsd=psd_computer.cross_psd(isk, isl);
 
+        TFXX_debug(opt.debug, "croposp",
+                   "CPSD " << TFXX_value(cpsd.data.size()));
+
         // compute array of frequency values, if not yet done
         if (frequencies.size()<1)
         {
@@ -516,10 +547,21 @@ int main(int iargc, char* argv[])
           named_cpsd.series=cpsd.data;
         } // if (opt.logscale) ... else
 
+        TFXX_debug(opt.debug, "croposp",
+                   "CPSD " << TFXX_value(named_cpsd.series.size()));
+
         named_cpsd.label=vector_of_metadata[k].label+" "
           +vector_of_metadata[l].label;
         CPSD_vector[pairs(k,l)]=named_cpsd;
 
+        TFXX_debug(opt.debug, "croposp",
+                   "CPSD " << TFXX_value(CPSD_vector[pairs(k,l)].label));
+
+        if (opt.verbose)
+        {
+          cout << "  computed: " << named_cpsd.label << endl;
+        }
+
       } // for (unsigned int l=0; l<k; ++l)
     } // for (unsigned int l=0; l<k; ++l)
 
@@ -533,15 +575,21 @@ int main(int iargc, char* argv[])
 
   if (opt.compute_coherency)
   {
-    croposp::TNamedSeriesVector coherency_vector;
-    croposp::Pairs pairs(collection_of_series.size());
+    if (opt.verbose)
+    {
+      cout << endl 
+        << "Compute coherency for all pairs of signals." 
+        << endl;
+    }
+
+    croposp::TNamedSeriesVector coherency_vector(pairs.pairs());
 
     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;
+    croposp::TNamedSeriesVector sqrtpsd_vector(pairs.size());
 
     for (unsigned int k=0; k<PSD_vector.size(); ++k)
     {
@@ -563,6 +611,11 @@ int main(int iargc, char* argv[])
         coherency_vector[pairs(k,l)].label=
           opt.prefix_coherency+" "+CPSD_vector[pairs(k,l)].label;
 
+        if (opt.verbose)
+        {
+          cout << "  computed: " << coherency_vector[pairs(k,l)].label << endl;
+        }
+
       } // for (unsigned int l=0; l<k; ++l)
     } // for (unsigned int l=0; l<k; ++l)