[FEATURE] (foutra): add option -derivative

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.

This option provides the caluclation of values representing the values for a
time series of which the derivative with time was taken. Support the
calculation of PSD for acceleration where input time series is velocity.

src/ts/wf/foutra.cc

@@ -51,6 +51,8 @@
  *  - 03/12/2014   V1.9   provide output file format selector
  *  - 08/01/2015   V1.10  FIX: when using -scalerbw take bandwidth from
  *                        opt.scaledecades rather than opt.decades
+ *  - 29/01/2015   V1.11  FEATURE: provide calculation of n-th derivative
+ *                        of time series
  *
  *  \note 08/01/2010:
  *  Scaling for foutra power spectrum was tested against theory.
@@ -107,8 +109,9 @@ struct Options {
   bool amplitudespectrum, powerspectrum, boxcartaper;
   bool avgconstbw, avgrelbw, avgasciionly;
   bool demean, detrend, scalerbw, adjustdivisor;
+  bool derivative;
   int ndemean, ndetrend, divisor;
-  double decades, scaledecades, asciidecades;
+  double decades, scaledecades, asciidecades, nderivative;
   int avgsamples;
   bool reportrms, harmonicsignal, padzeroes;
   int padfactor, nsegments;
@@ -142,6 +145,7 @@ int main(int iargc, char* argv[])
     "              [-boxcar] [-avg[=n]] [-rbw[=n]] [-avgascii]" "\n"
     "              [-demean[=n]] [-dtrend[=n]] [-scalerbw[=n]]" "\n"
     "              [-divisor[=n]] [-rms] [-harmonic] [-pad n]" "\n"
+    "              [-derivative n]" "\n"
     "              [-nsegments n]" "\n"
     "              outfile infile [t:T] [infile [t:T] ...]" "\n"
     "   or: foutra --help|-h" "\n"
@@ -173,6 +177,7 @@ int main(int iargc, char* argv[])
     "-rbw[=n]     smooth power spectrum by averaging over n decades" "\n"
     "-demean[=n]  remove average (determined from n samples)" "\n"
     "-detrend[=n] remove trend (determined from n samples)" "\n"
+    "-derivative[=n] take n-th derivative of time series" "\n"
     "-scalerbw[=n] scale to mean value in n decades" "\n"
     "-divisor[=n] FFT becomes very inefficient if the factorization" "\n"
     "             of the number of samples includes large prime numbers." "\n"
@@ -290,6 +295,8 @@ int main(int iargc, char* argv[])
     {"xhelp",arg_no,"-"},
     // 22: input format
     {"Type",arg_yes,"sff"},
+    // 23: input format
+    {"derivative",arg_opt,"1"},
     {NULL}
   };
 
@@ -350,6 +357,8 @@ int main(int iargc, char* argv[])
   opt.padfactor=cmdline.int_arg(19);
   opt.nsegments=cmdline.int_arg(20);
   opt.outputformat=cmdline.string_arg(22);
+  opt.derivative=cmdline.optset(23);
+  opt.nderivative=cmdline.double_arg(23);
 
   TFXX_assert((opt.divisor > 0), "illegal value for argument divisor");
   TFXX_assert((opt.nsegments > 0), "illegal value for argument nsegments");
@@ -434,6 +443,15 @@ int main(int iargc, char* argv[])
   if (opt.harmonicsignal || opt.amplitudespectrum || opt.powerspectrum)
   {
     processfree.append("An appropriately scaled FFT (libdrfftw) is applied");
+    if (opt.derivative)
+    {
+      processfree.append("Calculate values for derivative with respect to");
+      processfree.append("  time by multiplication of the Fourier coefficients");
+      std::ostringstream freeline;
+      freeline << "  with the angular frequency to the power of " 
+        << opt.nderivative;
+      processfree.append(freeline.str());
+    }
   }
   if (opt.amplitudespectrum)
   {
@@ -745,6 +763,19 @@ int main(int iargc, char* argv[])
             TFXX_debug(opt.debug, "main", "create iterators");
             aff::Iterator<Tseries> S(series);
             aff::Browser<Tfft::Tspectrum> C(coeff);
+
+            // take derivative if selected
+            if (opt.derivative)
+            {
+              for (int i=coeff.f()+1; i<=coeff.l(); ++i)
+              {
+                double frequency=df*(i-coeff.f())*2*M_PI;
+                double thepower=opt.nderivative;
+                double factor=std::pow(frequency,thepower);
+                coeff(i) *= factor;
+              }
+            }
+
             TFXX_debug(opt.debug, "main", "calculate square of modulus and copy");
             while(S.valid() && C.valid())
             {