From 499e173bac62e2d9d561b89a41b25ac6deeb1c5d Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Nov 2021 16:20:56 +0000 Subject: [PATCH 1/7] Update Readme.md --- Readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Readme.md b/Readme.md index 9137da0..c7079c3 100644 --- a/Readme.md +++ b/Readme.md @@ -1,6 +1,6 @@ # LatAnalyze -[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/10201777.svg)](https://zenodo.org/badge/latestdoi/10201777) +[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/10201777.svg)](https://zenodo.org/badge/latestdoi/10201777) [![Build Ubuntu](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml) [![Build macOS](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml) ## Description LatAnalyze is a C++11 library for statistical data analysis based on bootstrap From 65a656f2573d61010d7e8ba8ccc94e5fa10c6831 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 17 Feb 2022 19:24:34 +0000 Subject: [PATCH 2/7] first skeleton for DWT --- lib/Makefile.am | 4 + lib/Numerical/DWT.cpp | 57 ++++ lib/Numerical/DWT.hpp | 52 ++++ lib/Numerical/DWTFilters.cpp | 528 +++++++++++++++++++++++++++++++++++ lib/Numerical/DWTFilters.hpp | 53 ++++ 5 files changed, 694 insertions(+) create mode 100644 lib/Numerical/DWT.cpp create mode 100644 lib/Numerical/DWT.hpp create mode 100644 lib/Numerical/DWTFilters.cpp create mode 100644 lib/Numerical/DWTFilters.hpp 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..b5c49a7 --- /dev/null +++ b/lib/Numerical/DWT.cpp @@ -0,0 +1,57 @@ +/* + * 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 . + */ + +#include +#include + +using namespace std; +using namespace Latan; + +/****************************************************************************** + * DWT implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +DWT::DWT(const DWTFilter &filter) +: filter_(filter) +{} + +// convolution primitive /////////////////////////////////////////////////////// +DVec DWT::filterConvolution(const DVec &data, const DWTFilter &filter, + const Index offset) +{ + DVec res(data.size()); + + return res; +} + +// DWT ///////////////////////////////////////////////////////////////////////// +std::vector +DWT::forward(const DVec &data, const unsigned int level) const +{ + std::vector dwt(level); + + return dwt; +} + +DVec DWT::backward(const std::vector& dwt) const +{ + DVec res; + + return res; +} diff --git a/lib/Numerical/DWT.hpp b/lib/Numerical/DWT.hpp new file mode 100644 index 0000000..f1ee646 --- /dev/null +++ b/lib/Numerical/DWT.hpp @@ -0,0 +1,52 @@ +/* + * 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 . + */ + +#ifndef Latan_DWT_hpp_ +#define Latan_DWT_hpp_ + +#include +#include + +BEGIN_LATAN_NAMESPACE + +/****************************************************************************** + * Discrete wavelet transform class * + ******************************************************************************/ +class DWT +{ +public: + typedef std::pair DWTLevel; +public: + // constructor + DWT(const DWTFilter &filter); + // destructor + virtual ~DWT(void) = default; + // convolution primitive + static DVec filterConvolution(const DVec &data, const DWTFilter &filter, + const Index offset); + // DWT + std::vector forward(const DVec &data, const unsigned int level) const; + DVec backward(const std::vector& 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 . + */ + +#include +#include + +// 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 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 . + */ + +#ifndef Latan_DWTFilters_hpp_ +#define Latan_DWTFilters_hpp_ + +#include + +BEGIN_LATAN_NAMESPACE + +struct DWTFilter +{ + const std::vector 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 fromName; +} + +END_LATAN_NAMESPACE + +#endif // Latan_DWTFilters_hpp_ From 9e78b96260b5b57d84d146e2832d9c1dc56f80ad Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 18 Feb 2022 14:06:52 +0000 Subject: [PATCH 3/7] DWT working and tested --- examples/Makefile.am | 5 +++ examples/exDWT.cpp | 28 +++++++++++++ lib/Numerical/DWT.cpp | 92 ++++++++++++++++++++++++++++++++++++++++--- lib/Numerical/DWT.hpp | 7 +++- 4 files changed, 124 insertions(+), 8 deletions(-) create mode 100644 examples/exDWT.cpp 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 + +using namespace std; +using namespace Latan; + +int main(void) +{ + DVec data, dataRec; + vector 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/Numerical/DWT.cpp b/lib/Numerical/DWT.cpp index b5c49a7..bc0e03a 100644 --- a/lib/Numerical/DWT.cpp +++ b/lib/Numerical/DWT.cpp @@ -32,26 +32,106 @@ DWT::DWT(const DWTFilter &filter) {} // convolution primitive /////////////////////////////////////////////////////// -DVec DWT::filterConvolution(const DVec &data, const DWTFilter &filter, - const Index offset) +void DWT::filterConvolution(DVec &out, const DVec &data, + const std::vector &filter, const Index offset) { - DVec res(data.size()); + Index n = data.size(), nf = n*filter.size(); - return res; + 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::forward(const DVec &data, const unsigned int level) const { - std::vector dwt(level); + std::vector dwt(level); + DVec *finePt = const_cast(&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& dwt) const { - DVec res; + 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 index f1ee646..5745245 100644 --- a/lib/Numerical/DWT.hpp +++ b/lib/Numerical/DWT.hpp @@ -38,8 +38,11 @@ public: // destructor virtual ~DWT(void) = default; // convolution primitive - static DVec filterConvolution(const DVec &data, const DWTFilter &filter, - const Index offset); + static void filterConvolution(DVec &out, const DVec &data, + const std::vector &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 forward(const DVec &data, const unsigned int level) const; DVec backward(const std::vector& dwt) const; From 9afd40a1ad21ce0e6eabd8cdd7da0e12596893a3 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:16:40 +0000 Subject: [PATCH 4/7] utility to compute sample DWT --- utils/Makefile.am | 5 ++ utils/sample-dwt.cpp | 167 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 172 insertions(+) create mode 100644 utils/sample-dwt.cpp diff --git a/utils/Makefile.am b/utils/Makefile.am index 660649d..6abf910 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -9,6 +9,7 @@ endif bin_PROGRAMS = \ latan-plot \ latan-sample-combine \ + latan-sample-dwt \ latan-sample-element \ latan-sample-fake \ latan-sample-ft \ @@ -25,6 +26,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 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 . + */ + +#include +#include +#include +#include + +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 << " " << endl; + cerr << endl << "Possible options:" << endl << opt << endl; + + return EXIT_FAILURE; + } + inFilename = opt.getArgs()[0]; + level = opt.optionValue("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(inFilename), res; + Index nSample = in.size(), n = in[central].rows(); + vector out(ss ? 1 : level, DMatSample(nSample)), + outh(ss ? 0 : level, DMatSample(nSample)); + DWT dwt(*DWTFilters::fromName.at(filterName)); + vector 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(out[l], outFilename + "/L" + strFrom(l) + ".h5"); + Io::save(outh[l], outFilename + "/H" + strFrom(l) + ".h5"); + } + } + else + { + Io::save(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; +} From 43dd295f946c96fced794e5c6ed14d352aad3179 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:20:38 +0000 Subject: [PATCH 5/7] data plots compatible with multi-column arrays --- lib/Core/Plot.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) 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) From 9455e2c66eaca1684aa0b6e5c6da77a672afee97 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:21:10 +0000 Subject: [PATCH 6/7] correlation matrix plot compute dynamic range --- utils/sample-plot-corr.cpp | 1 + 1 file changed, 1 insertion(+) 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()) From 4f919bc00743bf2e893b9964e0634f725ebc380e Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:21:26 +0000 Subject: [PATCH 7/7] tool to merge samples --- utils/Makefile.am | 5 +++ utils/sample-merge.cpp | 94 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+) create mode 100644 utils/sample-merge.cpp diff --git a/utils/Makefile.am b/utils/Makefile.am index 660649d..70c1b87 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -12,6 +12,7 @@ bin_PROGRAMS = \ latan-sample-element \ latan-sample-fake \ latan-sample-ft \ + latan-sample-merge \ latan-sample-plot \ latan-sample-plot-corr\ latan-sample-read \ @@ -37,6 +38,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-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 . + */ + +#include +#include + +using namespace std; +using namespace Latan; + +int main(int argc, char *argv[]) +{ + // argument parsing //////////////////////////////////////////////////////// + OptParser opt; + bool parsed; + string outFileName = ""; + vector 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 << " ... " << 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(fileName[0]); + nSample = buf.size(); + totRows = buf[central].rows(); + nCol = buf[central].cols(); + for (unsigned int i = 1; i < n; ++i) + { + buf = Io::load(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(fileName[i]); + r = buf[central].rows(); + out.block(rowo, 0, r, nCol) = buf; + rowo += r; + } + if (!outFileName.empty()) + { + Io::save(out, outFileName); + } + + return EXIT_SUCCESS; +} \ No newline at end of file