libpsdxx [WP]: implement PSD computation

Code is copied from foutra.cc

src/libs/libpsdxx/dpsdcomputer_psd.cc

@@ -33,13 +33,232 @@
 #define PSDXX_DPSDCOMPUTER_PSD_CC_VERSION \
   "PSDXX_DPSDCOMPUTER_PSD_CC   V1.0"
 
-#include <psd.h>
+#include <iostream>
+#include <psdxx/psd.h>
+#include <psdxx/error.h>
+#include <aff/iterator.h>
+#include <aff/dump.h>
+#include <aff/seriesoperators.h>
+#include <aff/functions/avg.h>
+#include <aff/functions/rms.h>
+#include <aff/functions/sqrsum.h>
+#include <aff/subarray.h>
+#include <tsxx/tsxx.h>
+#include <tsxx/tapers.h>
+#include <tsxx/filter.h>
+#include <tsxx/wid2timeseries.h>
+#include <fourier/fftwaff.h>
+
+using std::cout;
+using std::endl;
+
+// time series
+typedef aff::Series<double> Tseries;
+
+typedef fourier::fft::DRFFTWAFF Tfft;
 
 namespace psd {
 
   Tdouble_series DPSDComputer::psd(const Tdouble_series& s) const
   {
     Tdouble_series retval;
+
+    /*======================================================================*/
+    // start processing
+
+    // create taper
+    ts::tapers::Hanning hanning_taper;
+    // create FFTW processor
+    Tfft fft;
+
+    const Tseries& inputseries=s.data;
+    const double& dt=s.interval;
+
+    // process data
+    //
+    // inputseries: full series passed to function
+    // series: series for one segment to be analysed 
+    //         the spectral representation of a single segment is also stored
+    //         in this container
+    // spectrum: resulting spectral representation
+    //
+    // in the previous (pre segmentation) version of foutra only "series" was
+    // used; we will make use of the aff::Series reference semantics to use the
+    // old code with no more modifications than necessary; such "series" will
+    // be used for different things
+    //
+    // wid2.dt remains valid throughout the processing
+    // wid2.nsamples is only correct for the inputseries
+    //   use series.size() for other request for number of samples
+
+    // prepare segmentation
+    // --------------------
+    int nsegsamples=inputseries.size();
+    if (this->Mnsegments>1)
+    {
+      nsegsamples=2*inputseries.size()/(this->Mnsegments+1); 
+    }
+
+    // adjust length of time series
+    if (this->Mdivisor>1)
+    {
+      if (this->Mverbose)
+      {
+        std::cout << "    adjust divisor to " << this->Mdivisor << std::endl;
+        std::cout << "    number of input samples:      " <<
+          inputseries.size() <<std::endl;
+        std::cout << "    number of segment samples:    " << nsegsamples
+          <<std::endl;
+      }
+      int rest=nsegsamples % this->Mdivisor;
+      nsegsamples=nsegsamples-rest;
+      if (this->Mverbose)
+      {
+        std::cout << "    number of processing samples: " << nsegsamples
+          <<std::endl;
+      }
+    }
+
+    // segment stride
+    int segstride=inputseries.size()/(this->Mnsegments+1);
+
+    if (this->Mverbose && (this->Mnsegments>1))
+    {
+      cout << "    number of input samples:           " 
+        << inputseries.size() << endl;
+      cout << "    number of segments:                " 
+        << this->Mnsegments << endl;
+      cout << "    number of samples in each segment: " 
+        << nsegsamples << endl;
+      cout << "    stride between segments:           "
+        << segstride << endl;
+      cout << "    overlap of proximate segments:     "
+        << 100.*(nsegsamples-segstride)/nsegsamples << "%" << endl;
+    }
+
+    // length of time series segment
+    double T=dt*nsegsamples;
+    // length of taper window
+    double Tw=T;
+    // frequency sampling interval
+    double df=1./T;
+
+    if (this->Mpadfactor>1)
+    {
+      T *= this->Mpadfactor;
+      df /= this->Mpadfactor;
+    }
+
+    if (this->Mverbose)
+    {
+      cout << "    duration T of each segment:     " 
+        << T << "s" << endl;
+      cout << "    duration Tw of signal window:   "
+        << Tw << "s" << endl;
+      cout << "    frequency sampling interval df: "
+        << df << "Hz" << endl;
+    }
+
+    // the result will be collected in
+    Tseries spectrum;
+
+    // segment counter
+    unsigned int isegment=0;
+
+    // start actual processing
+    // -----------------------
+    while (isegment < this->Mnsegments)
+    {
+      int firstindex=inputseries.first()+isegment*segstride;
+      int lastindex=firstindex+nsegsamples-1;
+      PSDXX_assert((firstindex >= inputseries.first()) &&
+                   (firstindex <= inputseries.last()) &&
+                   (lastindex >= inputseries.first()) &&
+                   (lastindex <= inputseries.last()),
+                   "ERROR: index out of range; program design error");
+      Tseries segseries=aff::subarray(inputseries)(firstindex,lastindex);
+      Tseries series=segseries;
+      if (this->Mnsegments>1) { series=segseries.copyout(); }
+      series.shift(-series.first());
+
+      // demean
+      if (this->Mdemean) { 
+        ts::filter::RemoveAverage demean(0);
+        demean(ts::filter::Ttimeseries(series, dt));
+      }
+
+      // detrend
+      if (this->Mdetrend) { 
+        ts::filter::RemoveTrend detrend(0);
+        detrend(ts::filter::Ttimeseries(series, dt));
+      }
+
+      // apply taper
+      hanning_taper.apply(series); 
+
+      // pad time series
+      // Notice: padding the signals requires correct scaling of amplitudes
+      // this is only implemented for harmonic signals up to now
+      if (this->Mpadfactor>0)
+      {
+        Tseries newseries(series.size()*this->Mpadfactor);
+        newseries=0.;
+        Tseries subseries(aff::subarray(newseries)(series.f(),series.l()));
+        subseries.copyin(series);
+        series=newseries;
+      }
+
+      Tfft::Tspectrum coeff=fft(series, dt);
+      series=Tseries(coeff.size());
+      aff::Iterator<Tseries> S(series);
+      aff::Browser<Tfft::Tspectrum> C(coeff);
+
+      while(S.valid() && C.valid())
+      {
+        *S = C->real()*C->real()+C->imag()*C->imag();
+        ++S; ++C;
+      }
+
+      // scaling factor to adjust for taper effect
+      double tapscaling=1.;
+      // scaling factor for Hanning
+      // see below for derivation of this factor
+      tapscaling=8./3.;
+
+      // we have an energy spectrum so far
+      // adjust scaling factor to obtain signal power
+      double scalingfactor=2.*tapscaling/Tw;
+
+      series *= scalingfactor;
+
+      // collect result
+      if (isegment==0)
+      {
+        spectrum=series.copyout();
+      }
+      else
+      {
+        // spectrum += series;
+        for (int i=spectrum.f(),
+             j=series.f();
+             i<=spectrum.l(); ++i, ++j)
+        {
+          PSDXX_assert(j<=series.l(),
+                       "ERROR: program design error");
+          spectrum(i) += series(j);
+        }
+      }
+
+      ++isegment;
+    } // while (isegment<this->Mnsegments)
+
+    if (this->Mnsegments>1) { spectrum /= this->Mnsegments; }
+
+    // use reference semantics to make results available under previous
+    // name
+    retval.data=spectrum;
+    retval.interval=df;
+
     return(retval);
   } // Tdouble_series DPSDComputer::psd(const Tdouble_series& s) const