mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2025-06-22 16:52:02 +01:00
Compare commits
1 Commits
andrew-pr
...
7ba88c496d
| Author | SHA1 | Date | |
|---|---|---|---|
| 7ba88c496d |
@ -184,7 +184,7 @@ PlotData::PlotData(const DMatSample &x, const DVec &y, const bool abs)
|
||||
}
|
||||
else
|
||||
{
|
||||
setCommand("'" + tmpFileName + "' u 1:(abs($3)):2 w xerr");
|
||||
setCommand("'" + tmpFileName + "' u 1:($3):2 w xerr");
|
||||
}
|
||||
}
|
||||
|
||||
@ -206,60 +206,6 @@ PlotData::PlotData(const XYStatData &data, const Index i, const Index j, const b
|
||||
setCommand("'" + tmpFileName + "' " + usingCmd);
|
||||
}
|
||||
|
||||
// PlotPoint constructor ///////////////////////////////////////////////////////
|
||||
PlotPoint::PlotPoint(const double x, const double y)
|
||||
{
|
||||
DMat d(1, 2);
|
||||
string usingCmd, tmpFileName;
|
||||
|
||||
d(0, 0) = x;
|
||||
d(0, 1) = y;
|
||||
tmpFileName = dumpToTmpFile(d);
|
||||
pushTmpFile(tmpFileName);
|
||||
setCommand("'" + tmpFileName + "' u 1:2");
|
||||
}
|
||||
|
||||
PlotPoint::PlotPoint(const DSample &x, const double y)
|
||||
{
|
||||
DMat d(1, 3);
|
||||
string usingCmd, tmpFileName;
|
||||
|
||||
d(0, 0) = x[central];
|
||||
d(0, 2) = y;
|
||||
d(0, 1) = sqrt(x.variance());
|
||||
tmpFileName = dumpToTmpFile(d);
|
||||
pushTmpFile(tmpFileName);
|
||||
setCommand("'" + tmpFileName + "' u 1:3:2 w xerr");
|
||||
}
|
||||
|
||||
PlotPoint::PlotPoint(const double x, const DSample &y)
|
||||
{
|
||||
DMat d(1, 3);
|
||||
string usingCmd, tmpFileName;
|
||||
|
||||
d(0, 0) = x;
|
||||
d(0, 1) = y[central];
|
||||
d(0, 2) = sqrt(y.variance());
|
||||
tmpFileName = dumpToTmpFile(d);
|
||||
pushTmpFile(tmpFileName);
|
||||
setCommand("'" + tmpFileName + "' u 1:2:3 w yerr");
|
||||
}
|
||||
|
||||
PlotPoint::PlotPoint(const DSample &x, const DSample &y)
|
||||
{
|
||||
DMat d(1, 4);
|
||||
string usingCmd, tmpFileName;
|
||||
|
||||
d(0, 0) = x[central];
|
||||
d(0, 2) = y[central];
|
||||
d(0, 1) = sqrt(x.variance());
|
||||
d(0, 3) = sqrt(y.variance());
|
||||
tmpFileName = dumpToTmpFile(d);
|
||||
pushTmpFile(tmpFileName);
|
||||
setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr");
|
||||
}
|
||||
|
||||
|
||||
// PlotLine constructor ////////////////////////////////////////////////////////
|
||||
PlotLine::PlotLine(const DVec &x, const DVec &y)
|
||||
{
|
||||
|
||||
@ -98,18 +98,6 @@ public:
|
||||
virtual ~PlotData(void) = default;
|
||||
};
|
||||
|
||||
class PlotPoint: public PlotObject
|
||||
{
|
||||
public:
|
||||
// constructor
|
||||
PlotPoint(const double x, const double y);
|
||||
PlotPoint(const DSample &x, const double y);
|
||||
PlotPoint(const double x, const DSample &y);
|
||||
PlotPoint(const DSample &x, const DSample &y);
|
||||
// destructor
|
||||
virtual ~PlotPoint(void) = default;
|
||||
};
|
||||
|
||||
class PlotHLine: public PlotObject
|
||||
{
|
||||
public:
|
||||
|
||||
@ -108,23 +108,6 @@ inline std::string strFrom(const T x)
|
||||
}
|
||||
|
||||
// specialization for vectors
|
||||
template<>
|
||||
inline std::vector<Index> strTo<std::vector<Index>>(const std::string &str)
|
||||
{
|
||||
std::vector<Index> res;
|
||||
std::vector<double> vbuf;
|
||||
double buf;
|
||||
std::istringstream stream(str);
|
||||
|
||||
while (!stream.eof())
|
||||
{
|
||||
stream >> buf;
|
||||
res.push_back(buf);
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
template<>
|
||||
inline DVec strTo<DVec>(const std::string &str)
|
||||
{
|
||||
|
||||
@ -253,39 +253,16 @@ DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
|
||||
}
|
||||
}
|
||||
|
||||
DMatSample CorrelatorUtils::fold(const DMatSample &c, const CorrelatorModels::ModelPar &par)
|
||||
DMatSample CorrelatorUtils::fold(const DMatSample &c)
|
||||
{
|
||||
const Index nt = c[central].rows();
|
||||
DMatSample buf = c;
|
||||
int sign;
|
||||
bool fold = false;
|
||||
|
||||
switch (par.type)
|
||||
{
|
||||
case CorrelatorType::cosh:
|
||||
case CorrelatorType::cst:
|
||||
sign = 1;
|
||||
fold = true;
|
||||
break;
|
||||
case CorrelatorType::sinh:
|
||||
sign = -1;
|
||||
fold = true;
|
||||
break;
|
||||
case CorrelatorType::linear:
|
||||
cout << "Linear model is asymmetric: will not fold." << endl;
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
if (fold)
|
||||
FOR_STAT_ARRAY(buf, s)
|
||||
{
|
||||
FOR_STAT_ARRAY(buf, s)
|
||||
for (Index t = 0; t < nt; ++t)
|
||||
{
|
||||
for (Index t = 0; t < nt; ++t)
|
||||
{
|
||||
buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt));
|
||||
}
|
||||
buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -56,7 +56,7 @@ namespace CorrelatorModels
|
||||
namespace CorrelatorUtils
|
||||
{
|
||||
DMatSample shift(const DMatSample &c, const Index ts);
|
||||
DMatSample fold(const DMatSample &c, const CorrelatorModels::ModelPar &par);
|
||||
DMatSample fold(const DMatSample &c);
|
||||
DMatSample fourierTransform(const DMatSample &c, FFT &fft,
|
||||
const unsigned int dir = FFT::Forward);
|
||||
};
|
||||
|
||||
@ -146,16 +146,6 @@ double Histogram::getX(const Index i) const
|
||||
return x_(i);
|
||||
}
|
||||
|
||||
double Histogram::getXMin(void) const
|
||||
{
|
||||
return xMin_;
|
||||
}
|
||||
|
||||
double Histogram::getXMax(void) const
|
||||
{
|
||||
return xMax_;
|
||||
}
|
||||
|
||||
double Histogram::operator[](const Index i) const
|
||||
{
|
||||
return bin_(i)*(isNormalized() ? norm_ : 1.);
|
||||
|
||||
@ -52,8 +52,6 @@ public:
|
||||
const StatArray<double> & getData(void) const;
|
||||
const StatArray<double> & getWeight(void) const;
|
||||
double getX(const Index i) const;
|
||||
double getXMin(void) const;
|
||||
double getXMax(void) const;
|
||||
double operator[](const Index i) const;
|
||||
double operator()(const double x) const;
|
||||
// percentiles & confidence interval
|
||||
|
||||
@ -300,67 +300,6 @@ const XYStatData & XYSampleData::getData(void)
|
||||
}
|
||||
|
||||
// fit /////////////////////////////////////////////////////////////////////////
|
||||
|
||||
void XYSampleData::fitSample(std::vector<Minimizer *> &minimizer,
|
||||
const std::vector<const DoubleModel *> &v,
|
||||
SampleFitResult &result,
|
||||
DVec &init,
|
||||
Index s)
|
||||
{
|
||||
result.resize(nSample_);
|
||||
result.chi2_.resize(nSample_);
|
||||
result.model_.resize(v.size());
|
||||
FitResult sampleResult;
|
||||
setDataToSample(s);
|
||||
if (s == central)
|
||||
{
|
||||
sampleResult = data_.fit(minimizer, init, v);
|
||||
init = sampleResult.segment(0, init.size());
|
||||
result.nPar_ = sampleResult.getNPar();
|
||||
result.nDof_ = sampleResult.nDof_;
|
||||
result.parName_ = sampleResult.parName_;
|
||||
result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
|
||||
}
|
||||
else
|
||||
{
|
||||
sampleResult = data_.fit(*(minimizer.back()), init, v);
|
||||
}
|
||||
result[s] = sampleResult;
|
||||
result.chi2_[s] = sampleResult.getChi2();
|
||||
for (unsigned int j = 0; j < v.size(); ++j)
|
||||
{
|
||||
result.model_[j].resize(nSample_);
|
||||
result.model_[j][s] = sampleResult.getModel(j);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
|
||||
const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v,
|
||||
Index s)
|
||||
{
|
||||
computeVarMat();
|
||||
|
||||
SampleFitResult result;
|
||||
DVec initCopy = init;
|
||||
|
||||
fitSample(minimizer, v, result, initCopy, s);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
SampleFitResult XYSampleData::fit(Minimizer &minimizer,
|
||||
const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v,
|
||||
Index s)
|
||||
{
|
||||
vector<Minimizer *> mv{&minimizer};
|
||||
|
||||
return fit(mv, init, v, s);
|
||||
}
|
||||
|
||||
SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
|
||||
const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v)
|
||||
@ -368,14 +307,43 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
|
||||
computeVarMat();
|
||||
|
||||
SampleFitResult result;
|
||||
FitResult sampleResult;
|
||||
DVec initCopy = init;
|
||||
Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity();
|
||||
|
||||
result.resize(nSample_);
|
||||
result.chi2_.resize(nSample_);
|
||||
result.model_.resize(v.size());
|
||||
FOR_STAT_ARRAY(result, s)
|
||||
{
|
||||
fitSample(minimizer, v, result, initCopy, s);
|
||||
setDataToSample(s);
|
||||
if (s == central)
|
||||
{
|
||||
sampleResult = data_.fit(minimizer, initCopy, v);
|
||||
initCopy = sampleResult.segment(0, initCopy.size());
|
||||
if (verbCopy != Minimizer::Verbosity::Debug)
|
||||
{
|
||||
minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
|
||||
}
|
||||
result[s] = sampleResult;
|
||||
result.chi2_[s] = sampleResult.getChi2();
|
||||
for (unsigned int j = 0; j < v.size(); ++j)
|
||||
{
|
||||
result.model_[j].resize(nSample_);
|
||||
result.model_[j][s] = sampleResult.getModel(j);
|
||||
}
|
||||
}
|
||||
minimizer.back()->setVerbosity(verbCopy);
|
||||
result.nPar_ = sampleResult.getNPar();
|
||||
result.nDof_ = sampleResult.nDof_;
|
||||
result.parName_ = sampleResult.parName_;
|
||||
result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
@ -103,16 +103,9 @@ public:
|
||||
// get internal XYStatData
|
||||
const XYStatData & getData(void);
|
||||
// fit
|
||||
void fitSample(std::vector<Minimizer *> &minimizer,
|
||||
const std::vector<const DoubleModel *> &v,
|
||||
SampleFitResult &sampleResult, DVec &init, Index s);
|
||||
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v, Index s);
|
||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v, Index s);
|
||||
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v);
|
||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v);
|
||||
template <typename... Ts>
|
||||
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
|
||||
@ -114,7 +114,6 @@ int main(int argc, char *argv[])
|
||||
nt = corr[central].rows();
|
||||
corr = corr.block(0, 0, nt, 1);
|
||||
corr = CorrelatorUtils::shift(corr, shift);
|
||||
|
||||
if (doLaplace)
|
||||
{
|
||||
vector<double> filter = {1., -2., 1.};
|
||||
@ -156,11 +155,6 @@ int main(int argc, char *argv[])
|
||||
}
|
||||
}
|
||||
|
||||
if (fold)
|
||||
{
|
||||
corr = CorrelatorUtils::fold(corr,modelPar);
|
||||
}
|
||||
|
||||
// fit /////////////////////////////////////////////////////////////////////
|
||||
DVec init(nPar);
|
||||
NloptMinimizer globMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
|
||||
|
||||
@ -6,17 +6,18 @@ if CXX_INTEL
|
||||
endif
|
||||
endif
|
||||
|
||||
bin_PROGRAMS = \
|
||||
latan-plot \
|
||||
latan-sample-combine \
|
||||
latan-sample-dwt \
|
||||
latan-sample-element \
|
||||
latan-sample-fake \
|
||||
latan-sample-ft \
|
||||
latan-sample-merge \
|
||||
latan-sample-plot \
|
||||
latan-sample-plot-corr\
|
||||
latan-sample-read \
|
||||
bin_PROGRAMS = \
|
||||
latan-plot \
|
||||
latan-sample-combine \
|
||||
latan-sample-dwt \
|
||||
latan-sample-element \
|
||||
latan-sample-fake \
|
||||
latan-sample-ft \
|
||||
latan-sample-merge \
|
||||
latan-sample-noise-analysis\
|
||||
latan-sample-plot \
|
||||
latan-sample-plot-corr \
|
||||
latan-sample-read \
|
||||
latan-resample
|
||||
|
||||
latan_plot_SOURCES = plot.cpp
|
||||
@ -47,6 +48,10 @@ latan_sample_merge_SOURCES = sample-merge.cpp
|
||||
latan_sample_merge_CXXFLAGS = $(COM_CXXFLAGS)
|
||||
latan_sample_merge_LDFLAGS = -L../lib/.libs -lLatAnalyze
|
||||
|
||||
latan_sample_noise_analysis_SOURCES = sample-noise-analysis.cpp
|
||||
latan_sample_noise_analysis_CXXFLAGS = $(COM_CXXFLAGS)
|
||||
latan_sample_noise_analysis_LDFLAGS = -L../lib/.libs -lLatAnalyze
|
||||
|
||||
latan_sample_plot_corr_SOURCES = sample-plot-corr.cpp
|
||||
latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS)
|
||||
latan_sample_plot_corr_LDFLAGS = -L../lib/.libs -lLatAnalyze
|
||||
|
||||
@ -18,42 +18,30 @@
|
||||
*/
|
||||
|
||||
#include <LatAnalyze/Io/Io.hpp>
|
||||
#include <LatAnalyze/Core/OptParser.hpp>
|
||||
|
||||
|
||||
using namespace std;
|
||||
using namespace Latan;
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
OptParser opt;
|
||||
Index nSample;
|
||||
double val, err;
|
||||
string outFileName;
|
||||
|
||||
opt.addOption("r", "seed" , OptParser::OptType::value, true,
|
||||
"random generator seed (default: random)");
|
||||
opt.addOption("", "help" , OptParser::OptType::trigger, true,
|
||||
"show this help message and exit");
|
||||
|
||||
bool parsed = opt.parse(argc, argv);
|
||||
if (!parsed or (opt.getArgs().size() != 4) or opt.gotOption("help"))
|
||||
if (argc != 5)
|
||||
{
|
||||
cerr << "usage: " << argv[0];
|
||||
cerr << " <central value> <error> <nSample> <output file>" << endl;
|
||||
cerr << endl << "Possible options:" << endl << opt << endl;
|
||||
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
|
||||
val = strTo<double>(argv[1]);
|
||||
err = strTo<double>(argv[2]);
|
||||
nSample = strTo<Index>(argv[3]);
|
||||
outFileName = argv[4];
|
||||
|
||||
random_device rd;
|
||||
SeedType seed = (opt.gotOption("r")) ? opt.optionValue<SeedType>("r") : rd();
|
||||
mt19937 gen(seed);
|
||||
mt19937 gen(rd());
|
||||
normal_distribution<> dis(val, err);
|
||||
DSample res(nSample);
|
||||
|
||||
@ -71,4 +59,4 @@ int main(int argc, char *argv[])
|
||||
Io::save<DSample>(res, outFileName);
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
}
|
||||
|
||||
193
utils/sample-noise-analysis.cpp
Normal file
193
utils/sample-noise-analysis.cpp
Normal file
@ -0,0 +1,193 @@
|
||||
/*
|
||||
* sample-noise-analysis.cpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2020 Antonin Portelli
|
||||
*
|
||||
* LatAnalyze 3 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 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* LatAnalyze 3 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 LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <LatAnalyze/Core/OptParser.hpp>
|
||||
#include <LatAnalyze/Io/Io.hpp>
|
||||
#include <LatAnalyze/Core/Math.hpp>
|
||||
#include <LatAnalyze/Core/Plot.hpp>
|
||||
#include <LatAnalyze/Numerical/GslFFT.hpp>
|
||||
#include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
|
||||
#include <LatAnalyze/Statistics/XYSampleData.hpp>
|
||||
|
||||
using namespace std;
|
||||
using namespace Latan;
|
||||
using namespace Math;
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// argument parsing ////////////////////////////////////////////////////////
|
||||
OptParser opt;
|
||||
bool parsed;
|
||||
string filename;
|
||||
|
||||
opt.addOption("" , "help" , OptParser::OptType::trigger, true,
|
||||
"show this help message and exit");
|
||||
parsed = opt.parse(argc, argv);
|
||||
if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help"))
|
||||
{
|
||||
cerr << "usage: " << argv[0];
|
||||
cerr << "<options> <sample file>" << endl;
|
||||
cerr << endl << "Possible options:" << endl << opt << endl;
|
||||
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
filename = opt.getArgs()[0];
|
||||
|
||||
// load data ///////////////////////////////////////////////////////////////
|
||||
DMatSample sample;
|
||||
|
||||
cout << "-- load data" << endl;
|
||||
sample = Io::load<DMatSample>(filename);
|
||||
|
||||
// compute power spectrum //////////////////////////////////////////////////
|
||||
DMat av, err;
|
||||
double l0;
|
||||
Index nSample = sample.size(), n = sample[central].rows();
|
||||
DMatSample noise(nSample), pow(nSample, n, 1);
|
||||
CMatSample ftBuf(nSample, n, 1);
|
||||
GslFFT fft(n);
|
||||
|
||||
cout << "-- compute power spectrum" << endl;
|
||||
FOR_STAT_ARRAY(sample, s)
|
||||
{
|
||||
sample[s].conservativeResize(n, 1);
|
||||
}
|
||||
av = sample.mean();
|
||||
err = sample.variance().cwiseSqrt();
|
||||
FOR_STAT_ARRAY(sample, s)
|
||||
{
|
||||
noise[s] = sample[s] - av;
|
||||
ftBuf[s].real() = noise[s];
|
||||
ftBuf[s].imag().fill(0.);
|
||||
fft(ftBuf[s]);
|
||||
pow[s] = ftBuf[s].cwiseAbs2().unaryExpr([](const double x){return 10.*log10(x);});
|
||||
//pow[s] = ftBuf[s].cwiseAbs2();
|
||||
pow[s].conservativeResize(n/2, 1);
|
||||
}
|
||||
pow[central] = pow.mean();
|
||||
// {
|
||||
// Plot p;
|
||||
// DVec x;
|
||||
|
||||
// x.setLinSpaced(n/2, 0., n/2 - 1.);
|
||||
// p << LogScale(Axis::x);
|
||||
// p << PlotData(x, pow);
|
||||
// p.display();
|
||||
// }
|
||||
// l0 = pow.mean()(1);
|
||||
// FOR_STAT_ARRAY(sample, s)
|
||||
// {
|
||||
// pow[s] = pow[s].unaryExpr([l0](const double x){return x - l0;});
|
||||
// }
|
||||
|
||||
// fit decay ///////////////////////////////////////////////////////////////
|
||||
DVec x, init(2);
|
||||
DMat fitErr;
|
||||
DMatSample xs(nSample, n/2, 1);
|
||||
DSample beta(nSample);
|
||||
XYSampleData data(nSample);
|
||||
MinuitMinimizer min;
|
||||
DoubleModel lin([](const double *x, const double *p){return x[0]*p[0] + p[1];}, 1, 2);
|
||||
|
||||
cout << "-- fit decay" << endl;
|
||||
x.setLinSpaced(n/2, 0., n/2 - 1.);
|
||||
FOR_VEC(x, i)
|
||||
{
|
||||
x(i) = log2(x(i));
|
||||
}
|
||||
xs.fill(x);
|
||||
data.addXDim(n/2, "f", true);
|
||||
data.addYDim("pow");
|
||||
data.setUnidimData(xs, pow);
|
||||
data.assumeYYCorrelated(true, 0, 0);
|
||||
for (unsigned int i = 0; i < n/2; ++i)
|
||||
{
|
||||
data.fitPoint((x(i) >= 2.) and (x(i) <= log2(n/2) - 0.5), i);
|
||||
}
|
||||
init(0) = -0.1; init(1) = -0.1;
|
||||
auto fit = data.fit(min, init, lin);
|
||||
fitErr = fit.variance().cwiseSqrt();
|
||||
FOR_STAT_ARRAY(beta, s)
|
||||
{
|
||||
beta[s] = fit[s](0)/(10.*log10(2.));
|
||||
}
|
||||
printf("chi^2/dof = %.1e/%d= %.2e -- chi^2 CCDF = %.2e -- p-value = %.2e -- CDR = %.1f dB\n",
|
||||
fit.getChi2(), static_cast<int>(fit.getNDof()), fit.getChi2PerDof(),
|
||||
fit.getCcdf(), fit.getPValue(), fit.getCorrRangeDb());
|
||||
printf(" decay = %.2f +/- %.2f dB/oct\n", fit[central](0), fitErr(0));
|
||||
printf(" exponent = %.2f +/- %.2f\n", beta[central], sqrt(beta.variance()));
|
||||
|
||||
Plot p;
|
||||
|
||||
p << Caption("noise power spectrum");
|
||||
p << PlotRange(Axis::x, -0.5, log2(n/2) + 0.5)
|
||||
<< Label("frequency (oct)", Axis::x) << Label("power (dB)", Axis::y);
|
||||
p << Color("1") << PlotPredBand(fit.getModel(_), 0., log2(n/2) + 0.5);
|
||||
p << Color("1") << PlotFunction(fit.getModel(), 0., log2(n/2) + 0.5);
|
||||
p << Color("2") << PlotData(x, pow);
|
||||
|
||||
p.display();
|
||||
|
||||
// p.reset();
|
||||
// p << PlotCorrMatrix(data.getFitCorrMat());
|
||||
// p.display();
|
||||
|
||||
// filter correlator ///////////////////////////////////////////////////////
|
||||
DVec filter(n);
|
||||
DMatSample fsample(nSample, n, 1);
|
||||
|
||||
FOR_VEC(filter, i)
|
||||
{
|
||||
filter(i) = -std::pow(2.*sin(pi/n*i), 2);//-beta[central]*.5);
|
||||
}
|
||||
FOR_STAT_ARRAY(sample, s)
|
||||
{
|
||||
ftBuf[s].real() = sample[s].col(0);
|
||||
ftBuf[s].imag().fill(0.);
|
||||
fft(ftBuf[s], FFT::Forward);
|
||||
ftBuf[s] = ftBuf[s].cwiseProduct(filter);
|
||||
fft(ftBuf[s], FFT::Backward);
|
||||
fsample[s] = ftBuf[s].real();
|
||||
}
|
||||
|
||||
// p.reset();
|
||||
x.setLinSpaced(n, 0., n - 1.);
|
||||
// p << PlotData(x, sample);
|
||||
// p << PlotData(x, fsample);
|
||||
// p.display();
|
||||
|
||||
p.reset();
|
||||
FOR_VEC(x, i)
|
||||
{
|
||||
x(i) = log2(x(i));
|
||||
}
|
||||
p << PlotRange(Axis::x, -0.5, log2(n/2) + 0.5);
|
||||
p << PlotPoints(x, -filter);
|
||||
p.display();
|
||||
p.reset();
|
||||
p << PlotCorrMatrix(sample.correlationMatrix());
|
||||
p.display();
|
||||
p.reset();
|
||||
p << PlotCorrMatrix(fsample.correlationMatrix());
|
||||
p.display();
|
||||
Io::save(fsample, "test.h5");
|
||||
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
@ -38,23 +38,9 @@ int main(int argc, char *argv[])
|
||||
{
|
||||
DMatSample s = Io::load<DMatSample>(fileName);
|
||||
string name = Io::getFirstName(fileName);
|
||||
Index nRows = s[central].rows();
|
||||
Index nCols = s[central].cols();
|
||||
cout << scientific;
|
||||
cout << "central value +/- standard deviation\n" << endl;
|
||||
cout << "Re:" << endl;
|
||||
for(Index i = 0; i < nRows; i++)
|
||||
{
|
||||
cout << s[central](i,0) << " +/- " << s.variance().cwiseSqrt()(i,0) << endl;
|
||||
}
|
||||
if(nCols == 2)
|
||||
{
|
||||
cout << "\nIm:" << endl;
|
||||
for(Index i = 0; i < nRows; i++)
|
||||
{
|
||||
cout << s[central](i,1) << " +/- " << s.variance().cwiseSqrt()(i,1) << endl;
|
||||
}
|
||||
}
|
||||
cout << "central value:\n" << s[central] << endl;
|
||||
cout << "standard deviation:\n" << s.variance().cwiseSqrt() << endl;
|
||||
if (!copy.empty())
|
||||
{
|
||||
Io::save(s, copy, File::Mode::write, name);
|
||||
@ -65,8 +51,8 @@ int main(int argc, char *argv[])
|
||||
DSample s = Io::load<DSample>(fileName);
|
||||
string name = Io::getFirstName(fileName);
|
||||
cout << scientific;
|
||||
cout << "central value +/- standard deviation\n" << endl;
|
||||
cout << s[central] << " +/- " << sqrt(s.variance()) << endl;
|
||||
cout << "central value:\n" << s[central] << endl;
|
||||
cout << "standard deviation:\n" << sqrt(s.variance()) << endl;
|
||||
if (!copy.empty())
|
||||
{
|
||||
Io::save(s, copy, File::Mode::write, name);
|
||||
|
||||
Reference in New Issue
Block a user