diff --git a/utils/Makefile.am b/utils/Makefile.am
index de68500..c39895f 100644
--- a/utils/Makefile.am
+++ b/utils/Makefile.am
@@ -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
diff --git a/utils/sample-noise-analysis.cpp b/utils/sample-noise-analysis.cpp
new file mode 100644
index 0000000..8cf9a4e
--- /dev/null
+++ b/utils/sample-noise-analysis.cpp
@@ -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 .
+ */
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+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 << " " << 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(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(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;
+}