diff --git a/Readme.md b/Readme.md
index 9137da0..c7079c3 100644
--- a/Readme.md
+++ b/Readme.md
@@ -1,6 +1,6 @@
# LatAnalyze
-[](https://www.gnu.org/licenses/gpl-3.0) [](https://zenodo.org/badge/latestdoi/10201777)
+[](https://www.gnu.org/licenses/gpl-3.0) [](https://zenodo.org/badge/latestdoi/10201777) [](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml) [](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml)
## Description
LatAnalyze is a C++11 library for statistical data analysis based on bootstrap
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 123bafb..bf23b0d 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -9,6 +9,7 @@ endif
noinst_PROGRAMS = \
exCompiledDoubleFunction\
exDerivative \
+ exDWT \
exFit \
exFitSample \
exIntegrator \
@@ -30,6 +31,10 @@ exDerivative_SOURCES = exDerivative.cpp
exDerivative_CXXFLAGS = $(COM_CXXFLAGS)
exDerivative_LDFLAGS = -L../lib/.libs -lLatAnalyze
+exDWT_SOURCES = exDWT.cpp
+exDWT_CXXFLAGS = $(COM_CXXFLAGS)
+exDWT_LDFLAGS = -L../lib/.libs -lLatAnalyze
+
exFit_SOURCES = exFit.cpp
exFit_CXXFLAGS = $(COM_CXXFLAGS)
exFit_LDFLAGS = -L../lib/.libs -lLatAnalyze
diff --git a/examples/exDWT.cpp b/examples/exDWT.cpp
new file mode 100644
index 0000000..082d8b5
--- /dev/null
+++ b/examples/exDWT.cpp
@@ -0,0 +1,28 @@
+#include <LatAnalyze/Numerical/DWT.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(void)
+{
+ DVec data, dataRec;
+ vector<DWT::DWTLevel> dataDWT;
+ DWT dwt(DWTFilters::db3);
+
+ cout << "-- random data" << endl;
+ data.setRandom(16);
+ cout << data.transpose() << endl;
+ cout << "-- compute Daubechies 3 DWT" << endl;
+ dataDWT = dwt.forward(data, 4);
+ for (unsigned int l = 0; l < dataDWT.size(); ++l)
+ {
+ cout << "* level " << l << endl;
+ cout << "L= " << dataDWT[l].first.transpose() << endl;
+ cout << "H= " << dataDWT[l].second.transpose() << endl;
+ }
+ cout << "-- check inverse DWT" << endl;
+ dataRec = dwt.backward(dataDWT);
+ cout << "rel diff = " << 2.*(data - dataRec).norm()/(data + dataRec).norm() << endl;
+
+ return EXIT_SUCCESS;
+}
diff --git a/lib/Core/Plot.cpp b/lib/Core/Plot.cpp
index 42848ed..3a236cd 100644
--- a/lib/Core/Plot.cpp
+++ b/lib/Core/Plot.cpp
@@ -122,10 +122,10 @@ PlotData::PlotData(const DMatSample &x, const DMatSample &y, const bool abs)
DMat d(x[central].rows(), 4);
string usingCmd, tmpFileName;
- d.col(0) = x[central];
- d.col(2) = y[central];
- d.col(1) = x.variance().cwiseSqrt();
- d.col(3) = y.variance().cwiseSqrt();
+ d.col(0) = x[central].col(0);
+ d.col(2) = y[central].col(0);
+ d.col(1) = x.variance().cwiseSqrt().col(0);
+ d.col(3) = y.variance().cwiseSqrt().col(0);
tmpFileName = dumpToTmpFile(d);
pushTmpFile(tmpFileName);
if (!abs)
@@ -149,8 +149,8 @@ PlotData::PlotData(const DVec &x, const DMatSample &y, const bool abs)
string usingCmd, tmpFileName;
d.col(0) = x;
- d.col(1) = y[central];
- d.col(2) = y.variance().cwiseSqrt();
+ d.col(1) = y[central].col(0);
+ d.col(2) = y.variance().cwiseSqrt().col(0);
tmpFileName = dumpToTmpFile(d);
pushTmpFile(tmpFileName);
if (!abs)
@@ -173,9 +173,9 @@ PlotData::PlotData(const DMatSample &x, const DVec &y, const bool abs)
DMat d(x[central].rows(), 3), xerr, yerr;
string usingCmd, tmpFileName;
- d.col(0) = x[central];
+ d.col(0) = x[central].col(0);
d.col(2) = y;
- d.col(1) = x.variance().cwiseSqrt();
+ d.col(1) = x.variance().cwiseSqrt().col(0);
tmpFileName = dumpToTmpFile(d);
pushTmpFile(tmpFileName);
if (!abs)
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 9fb5e8c..8f33d71 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -48,6 +48,8 @@ libLatAnalyze_la_SOURCES = \
Io/XmlReader.cpp \
Io/Xml/tinyxml2.cpp \
Numerical/Derivative.cpp \
+ Numerical/DWT.cpp \
+ Numerical/DWTFilters.cpp \
Numerical/GslFFT.cpp \
Numerical/GslHybridRootFinder.cpp\
Numerical/GslMinimizer.cpp \
@@ -92,6 +94,8 @@ HPPFILES = \
Io/IoObject.hpp \
Io/XmlReader.hpp \
Numerical/Derivative.hpp \
+ Numerical/DWT.hpp \
+ Numerical/DWTFilters.hpp \
Numerical/FFT.hpp \
Numerical/GslFFT.hpp \
Numerical/GslHybridRootFinder.hpp\
diff --git a/lib/Numerical/DWT.cpp b/lib/Numerical/DWT.cpp
new file mode 100644
index 0000000..bc0e03a
--- /dev/null
+++ b/lib/Numerical/DWT.cpp
@@ -0,0 +1,137 @@
+/*
+ * DWT.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 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 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. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/DWT.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ * DWT implementation *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+DWT::DWT(const DWTFilter &filter)
+: filter_(filter)
+{}
+
+// convolution primitive ///////////////////////////////////////////////////////
+void DWT::filterConvolution(DVec &out, const DVec &data,
+ const std::vector<double> &filter, const Index offset)
+{
+ Index n = data.size(), nf = n*filter.size();
+
+ out.resize(n);
+ out.fill(0.);
+ for (unsigned int i = 0; i < filter.size(); ++i)
+ {
+ FOR_VEC(out, j)
+ {
+ out(j) += filter[i]*data((j + i + nf - offset) % n);
+ }
+ }
+}
+
+// downsampling/upsampling primitives //////////////////////////////////////////
+void DWT::downsample(DVec &out, const DVec &in)
+{
+ if (out.size() < in.size()/2)
+ {
+ LATAN_ERROR(Size, "output vector smaller than half the input vector size");
+ }
+ for (Index i = 0; i < in.size(); i += 2)
+ {
+ out(i/2) = in(i);
+ }
+}
+
+void DWT::upsample(DVec &out, const DVec &in)
+{
+ if (out.size() < 2*in.size())
+ {
+ LATAN_ERROR(Size, "output vector smaller than twice the input vector size");
+ }
+ out.segment(0, 2*in.size()).fill(0.);
+ for (Index i = 0; i < in.size(); i ++)
+ {
+ out(2*i) = in(i);
+ }
+}
+
+// DWT /////////////////////////////////////////////////////////////////////////
+std::vector<DWT::DWTLevel>
+DWT::forward(const DVec &data, const unsigned int level) const
+{
+ std::vector<DWTLevel> dwt(level);
+ DVec *finePt = const_cast<DVec *>(&data);
+ DVec tmp;
+ Index n = data.size(), o = filter_.fwdL.size()/2, minSize;
+
+ minSize = 1;
+ for (unsigned int l = 0; l < level; ++l) minSize *= 2;
+ if (n < minSize)
+ {
+ LATAN_ERROR(Size, "data vector too small for a " + strFrom(level)
+ + "-level DWT (data size is " + strFrom(n) + ")");
+ }
+ for (unsigned int l = 0; l < level; ++l)
+ {
+ n /= 2;
+ dwt[l].first.resize(n);
+ dwt[l].second.resize(n);
+ filterConvolution(tmp, *finePt, filter_.fwdL, o);
+ downsample(dwt[l].first, tmp);
+ filterConvolution(tmp, *finePt, filter_.fwdH, o);
+ downsample(dwt[l].second, tmp);
+ finePt = &dwt[l].first;
+ }
+
+ return dwt;
+}
+
+DVec DWT::backward(const std::vector<DWTLevel>& dwt) const
+{
+ unsigned int level = dwt.size();
+ Index n = dwt.back().second.size(), o = filter_.bwdL.size()/2 - 1;
+ DVec res, tmp, conv;
+
+ res = dwt.back().first;
+ for (int l = level - 2; l >= 0; --l)
+ {
+ n *= 2;
+ if (dwt[l].second.size() != n)
+ {
+ LATAN_ERROR(Size, "DWT result size mismatch");
+ }
+ }
+ n = dwt.back().second.size();
+ for (int l = level - 1; l >= 0; --l)
+ {
+ n *= 2;
+ tmp.resize(n);
+ upsample(tmp, res);
+ filterConvolution(conv, tmp, filter_.bwdL, o);
+ res = conv;
+ upsample(tmp, dwt[l].second);
+ filterConvolution(conv, tmp, filter_.bwdH, o);
+ res += conv;
+ }
+
+ return res;
+}
diff --git a/lib/Numerical/DWT.hpp b/lib/Numerical/DWT.hpp
new file mode 100644
index 0000000..5745245
--- /dev/null
+++ b/lib/Numerical/DWT.hpp
@@ -0,0 +1,55 @@
+/*
+ * DWT.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 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 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. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DWT_hpp_
+#define Latan_DWT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Numerical/DWTFilters.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ * Discrete wavelet transform class *
+ ******************************************************************************/
+class DWT
+{
+public:
+ typedef std::pair<DVec, DVec> DWTLevel;
+public:
+ // constructor
+ DWT(const DWTFilter &filter);
+ // destructor
+ virtual ~DWT(void) = default;
+ // convolution primitive
+ static void filterConvolution(DVec &out, const DVec &data,
+ const std::vector<double> &filter, const Index offset);
+ // downsampling/upsampling primitives
+ static void downsample(DVec &out, const DVec &in);
+ static void upsample(DVec &out, const DVec &in);
+ // DWT
+ std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
+ DVec backward(const std::vector<DWTLevel>& dwt) const;
+private:
+ DWTFilter filter_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DWT_hpp_
diff --git a/lib/Numerical/DWTFilters.cpp b/lib/Numerical/DWTFilters.cpp
new file mode 100644
index 0000000..be3accf
--- /dev/null
+++ b/lib/Numerical/DWTFilters.cpp
@@ -0,0 +1,528 @@
+/*
+ * DWTFilters.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 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 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. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/DWTFilters.hpp>
+#include <LatAnalyze/includes.hpp>
+
+// cf. http://wavelets.pybytes.com
+// *here we implement the reverse filters more convenient for convolutions*
+
+using namespace std;
+using namespace Latan;
+
+#define FILTDICT(x) {#x, &DWTFilters::x}
+
+std::map<std::string, const DWTFilter *> DWTFilters::fromName = {
+ FILTDICT(haar),
+ FILTDICT(db2),
+ FILTDICT(db3),
+ FILTDICT(db4),
+ FILTDICT(db5),
+ FILTDICT(db6),
+ FILTDICT(bior13),
+ FILTDICT(bior15),
+ FILTDICT(bior22),
+ FILTDICT(bior24),
+ FILTDICT(bior31),
+ FILTDICT(bior33),
+ FILTDICT(bior35)
+};
+
+DWTFilter DWTFilters::haar = {
+ // fwdL
+ {0.7071067811865476,
+ 0.7071067811865476},
+ // fwdH
+ {0.7071067811865476,
+ -0.7071067811865476},
+ // bwdL
+ {0.7071067811865476,
+ 0.7071067811865476},
+ // bwdH
+ {-0.7071067811865476,
+ 0.7071067811865476}
+};
+
+DWTFilter DWTFilters::db2 = {
+ // fwdL
+ {0.48296291314469025,
+ 0.836516303737469,
+ 0.22414386804185735,
+ -0.12940952255092145},
+ // fwdH
+ {-0.12940952255092145,
+ -0.22414386804185735,
+ 0.836516303737469,
+ -0.48296291314469025},
+ // bwdL
+ {-0.12940952255092145,
+ 0.22414386804185735,
+ 0.836516303737469,
+ 0.48296291314469025},
+ // bwdH
+ {-0.48296291314469025,
+ 0.836516303737469,
+ -0.22414386804185735,
+ -0.12940952255092145}
+};
+
+DWTFilter DWTFilters::db3 = {
+ // fwdL
+ {0.3326705529509569,
+ 0.8068915093133388,
+ 0.4598775021193313,
+ -0.13501102001039084,
+ -0.08544127388224149,
+ 0.035226291882100656},
+ // fwdH
+ {0.035226291882100656,
+ 0.08544127388224149,
+ -0.13501102001039084,
+ -0.4598775021193313,
+ 0.8068915093133388,
+ -0.3326705529509569},
+ // bwdL
+ {0.035226291882100656,
+ -0.08544127388224149,
+ -0.13501102001039084,
+ 0.4598775021193313,
+ 0.8068915093133388,
+ 0.3326705529509569},
+ // bwdH
+ {-0.3326705529509569,
+ 0.8068915093133388,
+ -0.4598775021193313,
+ -0.13501102001039084,
+ 0.08544127388224149,
+ 0.035226291882100656}
+};
+
+DWTFilter DWTFilters::db4 = {
+ // fwdL
+ {0.23037781330885523,
+ 0.7148465705525415,
+ 0.6308807679295904,
+ -0.02798376941698385,
+ -0.18703481171888114,
+ 0.030841381835986965,
+ 0.032883011666982945,
+ -0.010597401784997278},
+ // fwdH
+ {-0.010597401784997278,
+ -0.032883011666982945,
+ 0.030841381835986965,
+ 0.18703481171888114,
+ -0.02798376941698385,
+ -0.6308807679295904,
+ 0.7148465705525415,
+ -0.23037781330885523},
+ // bwdL
+ {-0.010597401784997278,
+ 0.032883011666982945,
+ 0.030841381835986965,
+ -0.18703481171888114,
+ -0.02798376941698385,
+ 0.6308807679295904,
+ 0.7148465705525415,
+ 0.23037781330885523},
+ // bwdH
+ {-0.23037781330885523,
+ 0.7148465705525415,
+ -0.6308807679295904,
+ -0.02798376941698385,
+ 0.18703481171888114,
+ 0.030841381835986965,
+ -0.032883011666982945,
+ -0.010597401784997278}
+};
+
+DWTFilter DWTFilters::db5 = {
+ // fwdL
+ {0.160102397974125,
+ 0.6038292697974729,
+ 0.7243085284385744,
+ 0.13842814590110342,
+ -0.24229488706619015,
+ -0.03224486958502952,
+ 0.07757149384006515,
+ -0.006241490213011705,
+ -0.012580751999015526,
+ 0.003335725285001549},
+ // fwdH
+ {0.003335725285001549,
+ 0.012580751999015526,
+ -0.006241490213011705,
+ -0.07757149384006515,
+ -0.03224486958502952,
+ 0.24229488706619015,
+ 0.13842814590110342,
+ -0.7243085284385744,
+ 0.6038292697974729,
+ -0.160102397974125},
+ // bwdL
+ {0.003335725285001549,
+ -0.012580751999015526,
+ -0.006241490213011705,
+ 0.07757149384006515,
+ -0.03224486958502952,
+ -0.24229488706619015,
+ 0.13842814590110342,
+ 0.7243085284385744,
+ 0.6038292697974729,
+ 0.160102397974125},
+ // bwdH
+ {-0.160102397974125,
+ 0.6038292697974729,
+ -0.7243085284385744,
+ 0.13842814590110342,
+ 0.24229488706619015,
+ -0.03224486958502952,
+ -0.07757149384006515,
+ -0.006241490213011705,
+ 0.012580751999015526,
+ 0.003335725285001549}
+};
+
+DWTFilter DWTFilters::db6 = {
+ // fwdL
+ {0.11154074335008017,
+ 0.4946238903983854,
+ 0.7511339080215775,
+ 0.3152503517092432,
+ -0.22626469396516913,
+ -0.12976686756709563,
+ 0.09750160558707936,
+ 0.02752286553001629,
+ -0.031582039318031156,
+ 0.0005538422009938016,
+ 0.004777257511010651,
+ -0.00107730108499558},
+ // fwdH
+ {-0.00107730108499558,
+ -0.004777257511010651,
+ 0.0005538422009938016,
+ 0.031582039318031156,
+ 0.02752286553001629,
+ -0.09750160558707936,
+ -0.12976686756709563,
+ 0.22626469396516913,
+ 0.3152503517092432,
+ -0.7511339080215775,
+ 0.4946238903983854,
+ -0.11154074335008017},
+ // bwdL
+ {-0.00107730108499558,
+ 0.004777257511010651,
+ 0.0005538422009938016,
+ -0.031582039318031156,
+ 0.02752286553001629,
+ 0.09750160558707936,
+ -0.12976686756709563,
+ -0.22626469396516913,
+ 0.3152503517092432,
+ 0.7511339080215775,
+ 0.4946238903983854,
+ 0.11154074335008017},
+ // bwdH
+ {-0.11154074335008017,
+ 0.4946238903983854,
+ -0.7511339080215775,
+ 0.3152503517092432,
+ 0.22626469396516913,
+ -0.12976686756709563,
+ -0.09750160558707936,
+ 0.02752286553001629,
+ 0.031582039318031156,
+ 0.0005538422009938016,
+ -0.004777257511010651,
+ -0.00107730108499558}
+};
+
+DWTFilter DWTFilters::bior13 = {
+ // fwdL
+ {-0.08838834764831845,
+ 0.08838834764831845,
+ 0.7071067811865476,
+ 0.7071067811865476,
+ 0.08838834764831845,
+ -0.08838834764831845},
+ // fwdH
+ {0.0,
+ 0.0,
+ 0.7071067811865476,
+ -0.7071067811865476,
+ 0.0,
+ 0.0},
+ // bwdL
+ {0.0,
+ 0.0,
+ 0.7071067811865476,
+ 0.7071067811865476,
+ 0.0,
+ 0.0},
+ // bwdH
+ {0.08838834764831845,
+ 0.08838834764831845,
+ -0.7071067811865476,
+ 0.7071067811865476,
+ -0.08838834764831845,
+ -0.08838834764831845}
+};
+
+DWTFilter DWTFilters::bior15 = {
+ // fwdL
+ {0.01657281518405971,
+ -0.01657281518405971,
+ -0.12153397801643787,
+ 0.12153397801643787,
+ 0.7071067811865476,
+ 0.7071067811865476,
+ 0.12153397801643787,
+ -0.12153397801643787,
+ -0.01657281518405971,
+ 0.01657281518405971},
+ // fwdH
+ {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.7071067811865476,
+ -0.7071067811865476,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.0},
+ // bwdL
+ {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.7071067811865476,
+ 0.7071067811865476,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.0},
+ // bwdH
+ {-0.01657281518405971,
+ -0.01657281518405971,
+ 0.12153397801643787,
+ 0.12153397801643787,
+ -0.7071067811865476,
+ 0.7071067811865476,
+ -0.12153397801643787,
+ -0.12153397801643787,
+ 0.01657281518405971,
+ 0.01657281518405971}
+};
+
+DWTFilter DWTFilters::bior22 = {
+ // fwdL
+ {-0.1767766952966369,
+ 0.3535533905932738,
+ 1.0606601717798214,
+ 0.3535533905932738,
+ -0.1767766952966369,
+ 0.0},
+ // fwdH
+ {0.0,
+ 0.0,
+ 0.3535533905932738,
+ -0.7071067811865476,
+ 0.3535533905932738,
+ 0.0},
+ // bwdL
+ {0.0,
+ 0.0,
+ 0.3535533905932738,
+ 0.7071067811865476,
+ 0.3535533905932738,
+ 0.0},
+ // bwdH
+ {0.1767766952966369,
+ 0.3535533905932738,
+ -1.0606601717798214,
+ 0.3535533905932738,
+ 0.1767766952966369,
+ 0.0}
+};
+
+DWTFilter DWTFilters::bior24 = {
+ // fwdL
+ {0.03314563036811942,
+ -0.06629126073623884,
+ -0.1767766952966369,
+ 0.4198446513295126,
+ 0.9943689110435825,
+ 0.4198446513295126,
+ -0.1767766952966369,
+ -0.06629126073623884,
+ 0.03314563036811942,
+ 0.0},
+ // fwdH
+ {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.3535533905932738,
+ -0.7071067811865476,
+ 0.3535533905932738,
+ 0.0,
+ 0.0,
+ 0.0},
+ // bwdL
+ {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.3535533905932738,
+ 0.7071067811865476,
+ 0.3535533905932738,
+ 0.0,
+ 0.0,
+ 0.0},
+ // bwdH
+ {-0.03314563036811942,
+ -0.06629126073623884,
+ 0.1767766952966369,
+ 0.4198446513295126,
+ -0.9943689110435825,
+ 0.4198446513295126,
+ 0.1767766952966369,
+ -0.06629126073623884,
+ -0.03314563036811942,
+ 0.0}
+};
+
+DWTFilter DWTFilters::bior31 = {
+ // fwdL
+ {-0.3535533905932738,
+ 1.0606601717798214,
+ 1.0606601717798214,
+ -0.3535533905932738},
+ // fwdH
+ {0.1767766952966369,
+ -0.5303300858899107,
+ 0.5303300858899107,
+ -0.1767766952966369},
+ // bwdL
+ {0.1767766952966369,
+ 0.5303300858899107,
+ 0.5303300858899107,
+ 0.1767766952966369},
+ // bwdH
+ {0.3535533905932738,
+ 1.0606601717798214,
+ -1.0606601717798214,
+ -0.3535533905932738}
+};
+
+DWTFilter DWTFilters::bior33 = {
+ // fwdL
+ {0.06629126073623884,
+ -0.19887378220871652,
+ -0.15467960838455727,
+ 0.9943689110435825,
+ 0.9943689110435825,
+ -0.15467960838455727,
+ -0.19887378220871652,
+ 0.06629126073623884},
+ // fwdH
+ {0.0,
+ 0.0,
+ 0.1767766952966369,
+ -0.5303300858899107,
+ 0.5303300858899107,
+ -0.1767766952966369,
+ 0.0,
+ 0.0},
+ // bwdL
+ {0.0,
+ 0.0,
+ 0.1767766952966369,
+ 0.5303300858899107,
+ 0.5303300858899107,
+ 0.1767766952966369,
+ 0.0,
+ 0.0},
+ // bwdH
+ {-0.06629126073623884,
+ -0.19887378220871652,
+ 0.15467960838455727,
+ 0.9943689110435825,
+ -0.9943689110435825,
+ -0.15467960838455727,
+ 0.19887378220871652,
+ 0.06629126073623884}
+};
+
+DWTFilter DWTFilters::bior35 = {
+ // fwdL
+ {-0.013810679320049757,
+ 0.04143203796014927,
+ 0.052480581416189075,
+ -0.26792717880896527,
+ -0.07181553246425874,
+ 0.966747552403483,
+ 0.966747552403483,
+ -0.07181553246425874,
+ -0.26792717880896527,
+ 0.052480581416189075,
+ 0.04143203796014927,
+ -0.013810679320049757},
+ // fwdH
+ {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.1767766952966369,
+ -0.5303300858899107,
+ 0.5303300858899107,
+ -0.1767766952966369,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.0},
+ // bwdL
+ {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.1767766952966369,
+ 0.5303300858899107,
+ 0.5303300858899107,
+ 0.1767766952966369,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.0},
+ // bwdH
+ {0.013810679320049757,
+ 0.04143203796014927,
+ -0.052480581416189075,
+ -0.26792717880896527,
+ 0.07181553246425874,
+ 0.966747552403483,
+ -0.966747552403483,
+ -0.07181553246425874,
+ 0.26792717880896527,
+ 0.052480581416189075,
+ -0.04143203796014927,
+ -0.013810679320049757}
+};
diff --git a/lib/Numerical/DWTFilters.hpp b/lib/Numerical/DWTFilters.hpp
new file mode 100644
index 0000000..8ceab9a
--- /dev/null
+++ b/lib/Numerical/DWTFilters.hpp
@@ -0,0 +1,53 @@
+/*
+ * DWTFilters.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 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 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. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DWTFilters_hpp_
+#define Latan_DWTFilters_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+struct DWTFilter
+{
+ const std::vector<double> fwdL, fwdH, bwdL, bwdH;
+};
+
+namespace DWTFilters
+{
+ extern DWTFilter haar;
+ extern DWTFilter db2;
+ extern DWTFilter db3;
+ extern DWTFilter db4;
+ extern DWTFilter db5;
+ extern DWTFilter db6;
+ extern DWTFilter bior13;
+ extern DWTFilter bior15;
+ extern DWTFilter bior22;
+ extern DWTFilter bior24;
+ extern DWTFilter bior31;
+ extern DWTFilter bior33;
+ extern DWTFilter bior35;
+
+ extern std::map<std::string, const DWTFilter *> fromName;
+}
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DWTFilters_hpp_
diff --git a/utils/Makefile.am b/utils/Makefile.am
index 660649d..de68500 100644
--- a/utils/Makefile.am
+++ b/utils/Makefile.am
@@ -9,9 +9,11 @@ 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 \
@@ -25,6 +27,10 @@ latan_sample_combine_SOURCES = sample-combine.cpp
latan_sample_combine_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_combine_LDFLAGS = -L../lib/.libs -lLatAnalyze
+latan_sample_dwt_SOURCES = sample-dwt.cpp
+latan_sample_dwt_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_dwt_LDFLAGS = -L../lib/.libs -lLatAnalyze
+
latan_sample_element_SOURCES = sample-element.cpp
latan_sample_element_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_element_LDFLAGS = -L../lib/.libs -lLatAnalyze
@@ -37,6 +43,10 @@ latan_sample_ft_SOURCES = sample-ft.cpp
latan_sample_ft_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_ft_LDFLAGS = -L../lib/.libs -lLatAnalyze
+latan_sample_merge_SOURCES = sample-merge.cpp
+latan_sample_merge_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_merge_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-dwt.cpp b/utils/sample-dwt.cpp
new file mode 100644
index 0000000..958eea6
--- /dev/null
+++ b/utils/sample-dwt.cpp
@@ -0,0 +1,167 @@
+/*
+ * sample-dwt.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
+ *
+ * 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/Core/Plot.hpp>
+#include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Numerical/DWT.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+ // argument parsing ////////////////////////////////////////////////////////
+ OptParser opt;
+ bool parsed, doPlot, ss, saveAll;
+ unsigned int level;
+ string inFilename, outFilename, filterName;
+
+ opt.addOption("l", "level", OptParser::OptType::value, true,
+ "number of levels", "1");
+ opt.addOption("f", "filter", OptParser::OptType::value, true,
+ "filter name", "haar");
+ opt.addOption("s", "ss", OptParser::OptType::trigger, true,
+ "super-sampling (inverse DWT on data)");
+ opt.addOption("a", "all", OptParser::OptType::trigger, true,
+ "save all-levels");
+ opt.addOption("o", "output", OptParser::OptType::value, true,
+ "output file name, or directory name if --all is used (default: result not saved)", "");
+ opt.addOption("p", "plot", OptParser::OptType::trigger, true,
+ "show plot");
+ 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> <input file>" << endl;
+ cerr << endl << "Possible options:" << endl << opt << endl;
+
+ return EXIT_FAILURE;
+ }
+ inFilename = opt.getArgs()[0];
+ level = opt.optionValue<unsigned int>("l");
+ filterName = opt.optionValue("f");
+ ss = opt.gotOption("s");
+ saveAll = opt.gotOption("a");
+ outFilename = opt.optionValue("o");
+ doPlot = opt.gotOption("p");
+
+ // DWT /////////////////////////////////////////////////////////////////////
+ DMatSample in = Io::load<DMatSample>(inFilename), res;
+ Index nSample = in.size(), n = in[central].rows();
+ vector<DMatSample> out(ss ? 1 : level, DMatSample(nSample)),
+ outh(ss ? 0 : level, DMatSample(nSample));
+ DWT dwt(*DWTFilters::fromName.at(filterName));
+ vector<DWT::DWTLevel> dataDWT(level);
+
+ if (!ss)
+ {
+ cout << "-- compute discrete wavelet transform" << endl;
+ cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
+ FOR_STAT_ARRAY(in, s)
+ {
+ dataDWT = dwt.forward(in[s].col(0), level);
+ for (unsigned int l = 0; l < level; ++l)
+ {
+ out[l][s] = dataDWT[l].first;
+ outh[l][s] = dataDWT[l].second;
+ }
+ }
+ }
+ else
+ {
+ Index ssn = n;
+
+ cout << "-- compute inverse discrete wavelet transform" << endl;
+ cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
+ for (int l = level - 1; l >= 0; --l)
+ {
+ dataDWT[l].first.resize(ssn);
+ dataDWT[l].second.resize(ssn);
+ dataDWT[l].first.fill(0.);
+ dataDWT[l].second.fill(0.);
+ ssn *= 2;
+ }
+ FOR_STAT_ARRAY(in, s)
+ {
+ dataDWT.back().first = in[s].col(0);
+ out[0][s] = dwt.backward(dataDWT);
+ }
+ }
+ if (!outFilename.empty())
+ {
+ if (!ss and saveAll)
+ {
+ Latan::mkdir(outFilename);
+ for (unsigned int l = 0; l < level; ++l)
+ {
+ Io::save<DMatSample>(out[l], outFilename + "/L" + strFrom(l) + ".h5");
+ Io::save<DMatSample>(outh[l], outFilename + "/H" + strFrom(l) + ".h5");
+ }
+ }
+ else
+ {
+ Io::save<DMatSample>(out.back(), outFilename);
+ }
+ }
+
+ // plot ////////////////////////////////////////////////////////////////////
+ if (doPlot)
+ {
+ Plot p;
+ DVec x;
+
+ x.setLinSpaced(n, 0., n - 1.);
+ p << PlotRange(Axis::x, 0., n);
+ p << Title("original") << PlotData(x, in);
+ if (!ss)
+ {
+ Index ln = n, step = 1;
+
+ for (unsigned int l = 0; l < level; ++l)
+ {
+ ln /= 2;
+ step *= 2;
+ x.setLinSpaced(ln, 0., n - step);
+ p << Title("level " + strFrom(l + 1) + " L") << PlotData(x, out[l]);
+ p << Title("level " + strFrom(l + 1) + " H") << PlotData(x, outh[l]);
+ }
+ p.display();
+ }
+ else
+ {
+ double step = 1.;
+ DVec err;
+
+ for (unsigned int l = 0; l < level; ++l)
+ {
+ step /= 2.;
+ }
+ x.setLinSpaced(out[0][central].size(), 0., n - step);
+ err = out[0].variance().cwiseSqrt();
+ p << Title("supersampled") << Color("3") << PlotPredBand(x, out[0][central], err);
+ p << Color("3") << PlotLine(x, out[0][central]);
+ p.display();
+ }
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/utils/sample-merge.cpp b/utils/sample-merge.cpp
new file mode 100644
index 0000000..11c2e87
--- /dev/null
+++ b/utils/sample-merge.cpp
@@ -0,0 +1,94 @@
+/*
+ * sample-merge.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>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+ // argument parsing ////////////////////////////////////////////////////////
+ OptParser opt;
+ bool parsed;
+ string outFileName = "";
+ vector<string> fileName;
+ unsigned int n = 0;
+
+ opt.addOption("o", "output", OptParser::OptType::value , true,
+ "output file name (default: result not saved)");
+ 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 << " <sample 1> ... <sample n>" << endl;
+ cerr << endl << "Possible options:" << endl << opt << endl;
+
+ return EXIT_FAILURE;
+ }
+ n = opt.getArgs().size();
+ outFileName = opt.optionValue("o");
+ for (unsigned int i = 0; i < n; ++i)
+ {
+ fileName.push_back(opt.getArgs()[i]);
+ }
+
+ // figure out dimensions ///////////////////////////////////////////////////
+ Index nCol, nSample, totRows;
+ DMatSample buf;
+
+ buf = Io::load<DMatSample>(fileName[0]);
+ nSample = buf.size();
+ totRows = buf[central].rows();
+ nCol = buf[central].cols();
+ for (unsigned int i = 1; i < n; ++i)
+ {
+ buf = Io::load<DMatSample>(fileName[i]);
+ if (buf.size() != nSample)
+ {
+ LATAN_ERROR(Size, "sample size mismatch");
+ }
+ if (buf[central].cols() != nCol)
+ {
+ LATAN_ERROR(Size, "column size mismatch");
+ }
+ totRows += buf[central].rows();
+ }
+
+ // merge along rows ////////////////////////////////////////////////////////
+ DMatSample out(nSample, totRows, nCol);
+ Index rowo = 0, r;
+
+ for (unsigned int i = 0; i < n; ++i)
+ {
+ buf = Io::load<DMatSample>(fileName[i]);
+ r = buf[central].rows();
+ out.block(rowo, 0, r, nCol) = buf;
+ rowo += r;
+ }
+ if (!outFileName.empty())
+ {
+ Io::save<DMatSample>(out, outFileName);
+ }
+
+ return EXIT_SUCCESS;
+}
\ No newline at end of file
diff --git a/utils/sample-plot-corr.cpp b/utils/sample-plot-corr.cpp
index 02316c0..c900dc0 100644
--- a/utils/sample-plot-corr.cpp
+++ b/utils/sample-plot-corr.cpp
@@ -68,6 +68,7 @@ int main(int argc, char *argv[])
var = sample.varianceMatrix();
corr = sample.correlationMatrix();
+ cout << "dynamic range " << Math::svdDynamicRangeDb(corr) << " dB" << endl;
p << PlotCorrMatrix(corr);
p.display();
if (!outVarName.empty())