Commit 91b648f4 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

libpsdxx [WP][FEATURE]: provide reverse of logarithmic mapping

parent 9f63f413
......@@ -43,19 +43,77 @@ namespace psd {
namespace helper {
/*! Class function to return index on logarithmic scale
/*! \brief Function class to return index on logarithmic scale
*
* Support inverse operation to map index to lograrithmic frequency axis.
* The returned index value may be fractional, it is not an integer.
*
* Consider a series
* \f[
* f_l = a\, c^{l} + b
* \f]
* of frequency values on a logarithmic scale.
* Then
* \f[
* f_0 = a + b
* \f]
* \f[
* f_1 = a\, c + b,
* \f]
* and
* \f[
* f_2 = a\, c^{2} + b.
* \f]
* Hence
* \f[
* f_1 - f_0 = a\,(c-1)
* \f]
* and
* \f[
* f_2 - f_1 = a\,c\,(c-1)
* \f]
* such that
* \f[
* c=\frac{f_2-f_1}{f_1-f_0},
* \f]
* \f[
* a=\frac{f_1-f_0}{c-1},
* \f]
* and
* \f[
* b=f_0-a.
* \f]
* The index value of a given frequency \f$ f \f$ then is
* \f[
* l=\log_{c}\frac{f-b}{a}=\frac{\ln\left(\frac{f-b}{a}\right)}{\ln(c)}
* \f]
*/
class LogIndex
{
public:
/*! \brief Constructor.
* Initialize class with any three consecutive values of frequency on
* a logaritmic scale. The LogIndex::operator()() function will return
* 0 for \p f0.
*/
LogIndex(const double& f0,
const double& f1,
const double& f2);
/*! \brief Return index value for given frequency.
* Returns 0 for \p f0, 1 for \p f1, and 2 for \p f2, where \p f0, \p
* f1, and \p f2 are the arguments passed to the constructor.
*
* \param f frequency
* \return fractional index of frequency \p f
*/
double operator()(const double& f);
private:
double Mn_per_decade;
double Mc;
double Mdf;
//! factor in front of exponential function in mapping function
double Ma;
//! constant in mapping function
double Mb;
//! logarithm of base of exponential function
double Mlnc;
}; // class LogIndex
} // namespace helper
......
......@@ -46,10 +46,10 @@ namespace psd {
{
PSDXX_assert(f1>f0, "frequency values must increase");
PSDXX_assert(f2>f1, "frequency values must increase");
Mdf=f0;
double logfac=(f2-f1)/(f1-f0);
Mn_per_decade=std::log(10.)/std::log(logfac);
Mc=(f1-f0)/(f0*(logfac-1.));
double c=(f2-f1)/(f1-f0);
Ma=(f1-f0)/(c-1.);
Mb=f0-Ma;
Mlnc=std::log(c);
} // LogIndex::LogIndex(const double& f0,
// const double& f1,
// const double& f2)
......@@ -59,7 +59,8 @@ namespace psd {
double LogIndex::operator()(const double& f)
{
double l =Mn_per_decade*std::log(1.+((f/Mdf)-1.)/Mc)/std::log(10.);
double l = std::log((f-Mb)/Ma)/Mlnc;
return l;
} // double LogIndex::operator()(const double& f)
} // namespace helper
......
......@@ -21,6 +21,7 @@
#include <iostream>
#include <tfxx/commandline.h>
#include <psdxx/psd.h>
#include <psdxx/debug.h>
#include <psdxx/helper.h>
#include <aff/dump.h>
......@@ -37,10 +38,11 @@ struct Options
void testlogindex(const psd::TDseries::Tcoc& f, const int& i)
{
CODE( psd::helper::LogIndex li(f(i), f(i+1), f(i+2)) );
CODE( psd::helper::LogIndex index(f(i), f(i+1), f(i+2)) );
cout << "with " << PSDXX_value(i) << endl;
for (int i=f.f(); i<=f.l(); ++i)
{
cout << "i= " << i << " index(f(i))=" << li(f(i)) << endl;
cout << PSDXX_value(i) << " " << PSDXX_value(index(f(i))) << endl;
}
}
......
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