libpsdxx [WP]: implement actual processing

src/libs/libpsdxx/dpsdcomputer_processor.cc

@@ -33,7 +33,25 @@
 #define PSDXX_DPSDCOMPUTER_PROCESSOR_CC_VERSION \
   "PSDXX_DPSDCOMPUTER_PROCESSOR_CC   V1.0"
 
+#include <iostream>
 #include <psdxx/psd.h>
+#include <psdxx/error.h>
+#include <aff/iterator.h>
+#include <aff/seriesoperators.h>
+#include <aff/functions/min.h>
+#include <aff/subarray.h>
+#include <tsxx/tapers.h>
+#include <tsxx/filter.h>
+#include <fourier/fftwaff.h>
+
+using std::cout;
+using std::endl;
+
+// time series
+typedef ::psd::TDISeries::Tseries Tseries;
+typedef ::psd::TDCISeries::Tseries Tcseries;
+
+typedef fourier::fft::DRFFTWAFF Tfft;
 
 namespace psd {
 
@@ -81,16 +99,301 @@ namespace psd {
    * 
    */
   DPSDComputer::SpectralValues 
-    DPSDComputer::processor(const Tdouble_series& s1,
-                            const Tdouble_series& s2,
+    DPSDComputer::processor(const TDISeries::Tcoc& s1,
+                            const TDISeries::Tcoc& s2,
                             const bool& cpsd_flag,
                             const bool& psd_flag) const
     {
       DPSDComputer::SpectralValues retval;
+
+      PSDXX_assert(cpsd_flag || psd_flag,
+                   "invalid flag settings; this is a code bug");
+
+      /*======================================================================*/
+      // start processing
+
+      // create taper
+      ts::tapers::Hanning hanning_taper;
+      // create FFTW processor
+      Tfft fft;
+
+      const Tseries::Tcoc& inputseries1=s1.data;
+      const Tseries::Tcoc& inputseries2=s2.data;
+
+      const double& dt=s1.interval;
+
+      if (cpsd_flag)
+      {
+        PSDXX_assert(std::abs((s1.interval/s2.interval)-1.)<1.e-8,
+                     "input series have inconsistent interval");
+        PSDXX_assert(inputseries1.size() == inputseries2.size(),
+                     "input series have inconsistent size");
+      }
+
+      // process data
+      //
+      // inputseries1: full series passed to function
+      // inputseries2: 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
+      //
+      // dt remains valid throughout the processing
+
+      // prepare segmentation
+      // --------------------
+      int nsegsamples=inputseries1.size();
+      if (this->Mnsegments>1)
+      {
+        nsegsamples=2*inputseries1.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:      " <<
+            inputseries1.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=inputseries1.size()/(this->Mnsegments+1);
+
+      if (this->Mverbose && (this->Mnsegments>1))
+      {
+        cout << "    number of input samples:           " 
+          << inputseries1.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& psd1(retval.psd1.data);
+      Tseries& psd2(retval.psd2.data);
+      Tcseries& cpsd(retval.cpsd.data);
+      retval.psd1.interval=df;
+      retval.psd2.interval=df;
+      retval.cpsd.interval=df;
+
+      // prepare processing containers
+      // padding is done by filling paddedseries through a reference subarray
+      Tseries paddedseries1(nsegsamples*this->Mpadfactor);
+      Tseries paddedseries2(nsegsamples*this->Mpadfactor);
+      paddedseries1=0.;
+      paddedseries2=0.;
+      Tseries series1(aff::subarray(paddedseries1)(nsegsamples-1));
+      Tseries series2(aff::subarray(paddedseries2)(nsegsamples-1));
+
+      // segment counter
+      unsigned int isegment=0;
+
+      // start actual processing
+      // -----------------------
+      while (isegment < this->Mnsegments)
+      {
+        // processing containers
+
+        // extract segment of inputseries1
+        int firstindex=inputseries1.first()+isegment*segstride;
+        int lastindex=firstindex+nsegsamples-1;
+        PSDXX_assert((firstindex >= inputseries1.first()) &&
+                     (firstindex <= inputseries1.last()) &&
+                     (lastindex >= inputseries1.first()) &&
+                     (lastindex <= inputseries1.last()),
+                     "ERROR: index out of range; program design error");
+        Tseries::Tcoc segseries
+          =aff::subarray(inputseries1)(firstindex,lastindex);
+        series1.copyin(segseries);
+
+        // extract segment of inputseries2
+        if (cpsd_flag)
+        {
+          firstindex=inputseries2.first()+isegment*segstride;
+          lastindex=firstindex+nsegsamples-1;
+          PSDXX_assert((firstindex >= inputseries2.first()) &&
+                       (firstindex <= inputseries2.last()) &&
+                       (lastindex >= inputseries2.first()) &&
+                       (lastindex <= inputseries2.last()),
+                       "ERROR: index out of range; program design error");
+          segseries=aff::subarray(inputseries2)(firstindex,lastindex);
+          series2.copyin(segseries);
+        }
+
+        // demean
+        if (this->Mdemean) { 
+          ts::filter::RemoveAverage demean(0);
+          demean(ts::filter::Ttimeseries(series1, dt));
+          if (cpsd_flag) { demean(ts::filter::Ttimeseries(series2, dt)); }
+        }
+
+        // detrend
+        if (this->Mdetrend) { 
+          ts::filter::RemoveTrend detrend(0);
+          detrend(ts::filter::Ttimeseries(series1, dt));
+          if (cpsd_flag) { detrend(ts::filter::Ttimeseries(series2, dt)); }
+        }
+
+        // apply taper
+        hanning_taper.apply(series1); 
+        if (cpsd_flag) { hanning_taper.apply(series2); }
+
+        // compute FFT
+        Tfft::Tspectrum coeff1=fft(series1, dt);
+        Tfft::Tspectrum coeff2;
+        if (cpsd_flag) { coeff2=fft(series2, dt); }
+
+        // allocate space for results
+        if (isegment==0)
+        {
+          if (psd_flag) 
+          {
+            psd1=Tseries(coeff1.size());
+            psd1=0.;
+            if (cpsd_flag)
+            {
+              psd2=Tseries(coeff2.size());
+              psd2=0.;
+            }
+          }
+          if (cpsd_flag)
+          {
+            cpsd=Tcseries(coeff1.size());
+            cpsd=Tcseries::Tvalue(0.);
+          }
+        }
+        // end of allocation
+
+        // compute results
+        if (psd_flag) 
+        {
+          { // make iterators local
+            aff::Iterator<Tseries> S(psd1);
+            aff::Browser<Tfft::Tspectrum> C(coeff1);
+
+            while(S.valid() && C.valid())
+            {
+              *S = C->real()*C->real()+C->imag()*C->imag();
+              ++S; ++C;
+            }
+          }
+
+          if (cpsd_flag)
+          {
+            aff::Iterator<Tseries> S(psd2);
+            aff::Browser<Tfft::Tspectrum> C(coeff2);
+
+            while(S.valid() && C.valid())
+            {
+              *S = C->real()*C->real()+C->imag()*C->imag();
+              ++S; ++C;
+            }
+          }
+        }
+        if (cpsd_flag)
+        {
+            aff::Iterator<Tcseries> S(cpsd);
+            aff::Browser<Tfft::Tspectrum> C1(coeff1);
+            aff::Browser<Tfft::Tspectrum> C2(coeff2);
+
+            while(S.valid() && C1.valid() && C2.valid())
+            {
+              *S = *C1 * conj(*C2);
+              ++S; ++C1; ++C2;
+            }
+        }
+        // end of computation
+
+        ++isegment;
+      } // while (isegment<this->Mnsegments)
+
+      // scaling of results
+      // ------------------
+      // 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;
+
+      if (this->Mnsegments>1) { scalingfactor /= this->Mnsegments; }
+        
+      // apply scaling to results
+      if (psd_flag) 
+      {
+        psd1*=scalingfactor;
+        PSDXX_assert(aff::func::min(psd1)>=0.,
+                     "must be positive; this error indicates a bug");
+        if (cpsd_flag)
+        {
+          psd2*=scalingfactor;
+          PSDXX_assert(aff::func::min(psd2)>=0.,
+                       "must be positive; this error indicates a bug");
+        }
+      }
+      if (cpsd_flag)
+      {
+        cpsd*=scalingfactor;
+      }
+      // end of scaling
+
+
+      // use reference semantics to make results available under previous
+      // name
+
       return(retval);
     } /* DPSDComputer::SpectralValues 
-           DPSDComputer::processor(const Tdouble_series& s1,
-                                   const Tdouble_series& s2;
+           DPSDComputer::processor(const TDISeries::Tcoc& s1,
+                                   const TDISeries::Tcoc& s2;
                                    const bool& cpsd_flag
                                    const bool& psd_flag) const 
                                    */