croposp [WP][FEATURE]: computer phase transfer function

src/ts/croposp/croposp.cc

@@ -628,6 +628,62 @@ int main(int iargc, char* argv[])
 
   } // if (opt.compute_coherency)
 
+  /* ---------------------------------------------------------------------- */
+  /* compute amplitude transfer function
+   * -----------------------------------
+   */
+
+  if (opt.compute_phase)
+  {
+    if (opt.verbose)
+    {
+      cout << endl 
+        << "Compute amplitude phase transfer function for all pairs of signals." 
+        << endl;
+    }
+
+    TFXX_assert(pairs.size()>2, "requires at least three input time series");
+
+    croposp::Triples triples(pairs, true);
+    croposp::TNamedSeriesVector phase_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 k=1; k<pairs.size(); ++k)
+    {
+      for (unsigned int l=0; l<k; ++l)
+      {
+        for (unsigned int i=0; i<pairs.size(); ++i)
+        {
+          if ((k!=i) && (l!=i))
+          {
+            psd::TDseries::Tcoc cpsd1, cpsd2;
+            cpsd1=psd::arg(CPSD_vector[pairs(k,i)].series);
+            if (pairs.swap(k,i)) { cpsd1 = -1.*cpsd1; }
+            cpsd2=psd::arg(CPSD_vector[pairs(l,i)].series);
+            if (pairs.swap(l,i)) { cpsd2 = -1.*cpsd1; }
+            phase_vector[triples(i,k,l)].series=cpsd1-cpsd2;
+            phase_vector[triples(i,k,l)].label
+              =opt.prefix_phase+" "+vector_of_metadata[k].label
+              +"-"+vector_of_metadata[l].label+" (ref: "
+              +vector_of_metadata[i].label+")";
+          }
+        }
+      }
+    }
+
+    write_named_series(opt.outfile_phase,
+                       "phase transfer function of pairs of signals",
+                       frequencies,
+                       phase_vector,
+                       opt.verbose,
+                       opt.overwrite);
+
+  } // if (opt.compute_phase)
+
   /* ---------------------------------------------------------------------- */
   /* compute amplitude transfer function
    * -----------------------------------
@@ -673,7 +729,14 @@ int main(int iargc, char* argv[])
       }
     }
 
-  } // if (opt.compute_coherency)
+    write_named_series(opt.outfile_transfer,
+                       "amplitude transfer ratio of pairs of signals",
+                       frequencies,
+                       transfer_vector,
+                       opt.verbose,
+                       opt.overwrite);
+
+  } // if (opt.compute_transfer)
 
 } // main()