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

libpsdxx [FEATURE]: provide logarithmic mapping of raw series

parent 7083166f
......@@ -87,6 +87,7 @@ LIBCCSRC= \
function_log_frequency.cc \
function_log_sampling_c.cc \
function_log_sampling_r.cc \
function_lin_log_sampling.cc \
function_sqrt.cc \
log_index.cc \
log_taper.cc \
......
/*! \file function_lin_log_sampling.cc
* \brief take a series with linear sampling and resample to logarithmic sampling (implementation)
*
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 14/02/2019
*
* take a series with linear sampling and resample to logarithmic sampling
* (implementation)
*
* Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
* ----
*
* REVISIONS and CHANGES
* - 14/02/2019 V1.0 Thomas Forbriger
*
* ============================================================================
*/
#define PSDXX_FUNCTION_LIN_LOG_SAMPLING_CC_VERSION \
"PSDXX_FUNCTION_LIN_LOG_SAMPLING_CC V1.0"
#include <psdxx/psd.h>
#include <psdxx/error.h>
namespace psd {
TDseries log_sampling(const TDseries::Tcoc& s,
const TDseries::Tcoc& flin,
const TDseries::Tcoc& flog)
{
// consistency checks
PSDXX_assert(s.size()>1, "input series must contain at least 2 values");
PSDXX_assert(s.size()==flin.size(),
"series of frequencies does not match data series");
double dt1=flin(flin.f()+1)-flin(flin.f());
double dt2=flin(flin.l())-flin(flin.l()-1);
PSDXX_assert(std::abs((dt2/dt1)-1.)<1.e-8,
"input frequency scale is not linear");
TDISeries::Tcoc data;
data.interval=dt1;
data.data=s;
TDseries retval=log_sampling(data, flog);
return(retval);
} // TDseries log_sampling(const TDseries::Tcoc& s,
// const TDseries::Tcoc& flin,
// const TDseries::Tcoc& flog)
} // namespace psd
/* ----- END OF function_lin_log_sampling.cc ----- */
......@@ -292,6 +292,21 @@ namespace psd {
TDISeries::Tseries log_sampling(const TDISeries::Tcoc& s,
const TDseries::Tcoc& f);
/*! sample real values on a logarithmic frequency axis
*
* \param s spectral values with uniform sampling along frequency axis
* \param flin frequency values on a linear scale, matching the input
* series \p s
* \param flog frequency values on a logarithmic scale like returned by
* function log_frequency()
* \return series of sampled spectral values, where each sample of the
* returned series corresponds to the frequency with the
* corresponding index in \p flog
*/
TDseries log_sampling(const TDseries::Tcoc& s,
const TDseries::Tcoc& flin,
const TDseries::Tcoc& flog);
} // namespace psd
#endif // PSDXX_PSD_H_VERSION (includeguard)
......
......@@ -70,6 +70,7 @@ int main(int iargc, char* argv[])
CROPOSP_VERSION "\n"
"usage: croposp [-verbose] [-itype f] [-trim] [-datetolerance t]" "\n"
" [-pattern p] [-overwrite] [-DEBUG]" "\n"
" [-log n] [-avgout]\n"
" [-nsegments n] [-divisor n] [-padfactor n] [-overlap f]\n"
" [-psd f] [-npsd f] [-transfer f] [-phase f] [-coherency f]" "\n"
" [-prefixpsd s] [-prefixnpsd s] [-prefixtransfer s]\n"
......@@ -105,6 +106,13 @@ int main(int iargc, char* argv[])
"computation options:\n"
"-log n map averages to logarithmic scale with \"n\"\n"
" samples pre decade\n"
"-avgout do mapping on logarithmic scale not for values\n"
" of raw PSD and cross-PSD but for final results of\n"
" computation; in cases where computation is based on\n"
" cross power spectra density and where the differential\n"
" phase response of input channels strongly varies with\n"
" frequency, results may differ; using this option\n"
" increases computation time\n"
"-nsegments n set number of segments to be used for each input signal\n"
"-divisor n set number of samples to an integer multiple of \"n\"\n"
"-padfactor n pad with zeroes to increase length of series by a factor\n"
......@@ -182,6 +190,8 @@ int main(int iargc, char* argv[])
{"prefixphase",arg_yes,"PHA"},
// 22: set coherency prefix
{"prefixcoherency",arg_yes,"COH"},
// 23: average upon output
{"avgout",arg_no,"-"},
{NULL}
};
......@@ -254,6 +264,8 @@ int main(int iargc, char* argv[])
opt.prefix_phase=cmdline.string_arg(21);
opt.prefix_coherency=cmdline.string_arg(22);
opt.avg_on_output=cmdline.optset(23);
TFXX_assert(opt.n_per_decade>0,
"number of samples per decade must be finite and positive");
TFXX_assert(opt.nsegments>0,
......@@ -366,6 +378,17 @@ int main(int iargc, char* argv[])
TFXX_assert(collection_of_series.are_consistent(comparer),
"time series are inconsistent");
if (opt.compute_npsd || opt.compute_coherency)
{
if ((opt.nsegments < 2) && (!opt.logscale))
{
cout << "WARNING! The requested result requires averaging." << endl;
cout << " Current parameters avoid averaging of" << endl;
cout << " cross power spectral density. Results" << endl;
cout << " might be meaningless." << endl;
}
}
/* ====================================================================== */
// provide container for frequency values (will be filled upon computation
......@@ -396,7 +419,12 @@ int main(int iargc, char* argv[])
<< "Compute power spectral density for all signals." << endl;
if (opt.logscale)
{
cout << "Map to logarithmic frequency scale." << endl;
cout << "Map to logarithmic frequency scale";
if (opt.avg_on_output)
{
cout << " upon output to file";
}
cout << "." << endl;
}
}
......@@ -416,7 +444,7 @@ int main(int iargc, char* argv[])
if (i_series==collection_of_series.begin())
{
if (opt.logscale)
if (opt.logscale && !opt.avg_on_output)
{
frequencies=psd::log_frequency(psd.interval, psd.data.size(),
opt.n_per_decade);
......@@ -424,18 +452,18 @@ int main(int iargc, char* argv[])
else
{
frequencies=psd::lin_frequency(psd.interval, psd.data.size());
} // if (opt.logscale) ... else
} // if (opt.logscale && !opt.avg_on_output) ... else
} // if (i_series==collection_of_series.begin())
croposp::NamedSeries named_psd;
if (opt.logscale)
if (opt.logscale && !opt.avg_on_output)
{
named_psd.series=psd::log_sampling(psd, frequencies);
} // if (opt.logscale)
} // if (opt.logscale && !opt.avg_on_output)
else
{
named_psd.series=psd.data;
} // if (opt.logscale) ... else
} // if (opt.logscale && !opt.avg_on_output) ... else
named_psd.label=opt.prefix_psd+" "+i_meta->label;
PSD_vector.push_back(named_psd);
......@@ -475,7 +503,12 @@ int main(int iargc, char* argv[])
<< endl;
if (opt.logscale)
{
cout << "Map to logarithmic frequency scale." << endl;
cout << "Map to logarithmic frequency scale";
if (opt.avg_on_output)
{
cout << " upon output to file";
}
cout << "." << endl;
}
}
......@@ -503,7 +536,7 @@ int main(int iargc, char* argv[])
// compute array of frequency values, if not yet done
if (frequencies.size()<1)
{
if (opt.logscale)
if (opt.logscale && !opt.avg_on_output)
{
frequencies=psd::log_frequency(cpsd.interval, cpsd.data.size(),
opt.n_per_decade);
......@@ -511,18 +544,18 @@ int main(int iargc, char* argv[])
else
{
frequencies=psd::lin_frequency(cpsd.interval, cpsd.data.size());
} // if (opt.logscale) ... else
} // if (opt.logscale && !opt.avg_on_output) ... else
} // if (frequencies.size()<1)
croposp::NamedCPSD named_cpsd;
if (opt.logscale)
if (opt.logscale && !opt.avg_on_output)
{
named_cpsd.series=psd::log_sampling(cpsd, frequencies);
} // if (opt.logscale)
} // if (opt.logscale && !opt.avg_on_output)
else
{
named_cpsd.series=cpsd.data;
} // if (opt.logscale) ... else
} // if (opt.logscale && !opt.avg_on_output) ... else
TFXX_debug(opt.debug, "croposp",
"CPSD " << TFXX_value(named_cpsd.series.size()));
......
......@@ -54,7 +54,7 @@ namespace croposp {
bool verbose, trim, debug, overwrite;
std::string inputformat, labelpattern;
double datetolerance;
bool logscale;
bool logscale, avg_on_output;
unsigned int n_per_decade;
unsigned int nsegments, divisor, padfactor;
double overlap;
......
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