Skip to content

GitLab

Commit ad8d892a authored by thomas.forbriger's avatar thomas.forbriger
Browse files

libpsdxx [WP][FEATURE]: use IntervalSeries struct with Tcoc

parent 5b70da52
Hide whitespace changes
Inline Side-by-side
src/libs/libpsdxx/dpsdcomputer_coherency.cc
View file @ ad8d892a
......@@ -37,12 +37,12 @@
namespace psd {
Tdouble_series DPSDComputer::coherency(const Tdouble_series& s1,
const Tdouble_series& s2) const
TDISeries DPSDComputer::coherency(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const
{
Tdouble_series retval;
TDISeries retval;
return(retval);
} // Tdouble_series DPSDComputer::coherency(...) const
} // TDISeries DPSDComputer::coherency(...) const
} // namespace psd
......
src/libs/libpsdxx/dpsdcomputer_cross_psd.cc
View file @ ad8d892a
......@@ -37,13 +37,13 @@
namespace psd {
Tcomplex_double_series
DPSDComputer::cross_psd(const Tdouble_series& s1,
const Tdouble_series& s2) const
TDCISeries
DPSDComputer::cross_psd(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const
{
Tcomplex_double_series retval;
TDCISeries retval;
return(retval);
} // Tcomplex_double_series DPSDComputer::cross_psd(...) const
} // TDCISeries DPSDComputer::cross_psd(...) const
} // namespace psd
......
src/libs/libpsdxx/dpsdcomputer_psd.cc
View file @ ad8d892a
......@@ -33,234 +33,15 @@
#define PSDXX_DPSDCOMPUTER_PSD_CC_VERSION \
"PSDXX_DPSDCOMPUTER_PSD_CC V1.0"
#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
TDISeries DPSDComputer::psd(const TDISeries::Tcoc& 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;
TDISeries retval;
return(retval);
} // Tdouble_series DPSDComputer::psd(const Tdouble_series& s) const
} // TDISeries DPSDComputer::psd(const TDISeries::Tcoc& s) const
} // namespace psd
......
src/libs/libpsdxx/function_abs.cc
View file @ ad8d892a
......@@ -38,12 +38,14 @@
namespace psd {
TDIseries abs(const TDCIseries& s)
TDISeries abs(const TDCISeries::Tcoc& s)
{
TDIseries retval(s);
TDISeries retval;
retval.interval=s.interval;
retval.data=TDISeries::Tseries(s.data.shape());
aff::Iterator<TDIseries::Tseries> S(retval.data);
aff::Browser<TDCIseries::Tseries> C(s.data);
aff::Iterator<TDISeries::Tseries> S(retval.data);
aff::Browser<TDCISeries::Tcoc::Tseries> C(s.data);
while(S.valid() && C.valid())
{
......@@ -51,7 +53,7 @@ namespace psd {
++S; ++C;
}
return(retval);
} // TDIseries abs(const TDCIseries& s)
} // TDISeries abs(const TDCISeries::Tcoc& s)
} // namespace psd
......
src/libs/libpsdxx/function_arg.cc
View file @ ad8d892a
......@@ -38,12 +38,14 @@
namespace psd {
TDIseries arg(const TDCIseries& s)
TDISeries arg(const TDCISeries::Tcoc& s)
{
TDIseries retval(s);
TDISeries retval;
retval.interval=s.interval;
retval.data=TDISeries::Tseries(s.data.shape());
aff::Iterator<TDIseries::Tseries> S(retval.data);
aff::Browser<TDCIseries::Tseries> C(s.data);
aff::Iterator<TDISeries::Tseries> S(retval.data);
aff::Browser<TDCISeries::Tcoc::Tseries> C(s.data);
while(S.valid() && C.valid())
{
......@@ -51,7 +53,7 @@ namespace psd {
++S; ++C;
}
return(retval);
} // TDIseries arg(const TDCIseries& s)
} // TDISeries arg(const TDCISeries::Tcoc& s)
} // namespace psd
......
src/libs/libpsdxx/psd.h
View file @ ad8d892a
......@@ -120,15 +120,15 @@ namespace psd {
//! time series and real spectral values (PSD and coherency)
typedef IntervalSeries<double> TDISeries;
//! complex coefficients
typedef IntervalSeries<std::complex<double> > TDCIseries;
typedef IntervalSeries<std::complex<double> > TDCISeries;
/* ====================================================================== */
//! return absolute value (magnitude) of complex values
TDIseries abs(const TDCIseries& s);
TDISeries abs(const TDCISeries::Tcoc& s);
//! return argument of complex values
TDIseries arg(const TDCIseries& s);
TDISeries arg(const TDCISeries::Tcoc& s);
/* ====================================================================== */
......@@ -164,13 +164,13 @@ namespace psd {
bool verbose() const { return(Mverbose); }
//! compute power spectral density
TDIseries psd(const TDIseries& s) const;
TDISeries psd(const TDISeries::Tcoc& s) const;
//! compute cross power spectrum
TDCIseries cross_psd(const TDIseries& s1,
const TDIseries& s2) const;
TDCISeries cross_psd(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const;
//! compute coherency
TDIseries coherency(const TDIseries& s1,
const TDIseries& s2) const;
TDISeries coherency(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2) const;
private:
//! number of segments split input time series
......@@ -191,13 +191,13 @@ namespace psd {
//! container for return value of general processor (internal use only)
struct SpectralValues
{
TDIseries psd1, psd2;
TDCIseries cpsd;
TDISeries psd1, psd2;
TDCISeries cpsd;
}; // struct SpectralValues
//! general processor (internal use only)
SpectralValues processor(const TDIseries& s1,
const TDIseries& s2,
SpectralValues processor(const TDISeries::Tcoc& s1,
const TDISeries::Tcoc& s2,
const bool& cpsd_flag,
const bool& psd_flag) const;
}; // class DPSDComputer
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Imprint