portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2025-11-04 08:04:32 +00:00
Code Releases Activity

Merge branch 'feature/dwt' into develop

Browse Source
This commit is contained in:
Antonin Portelli
2022-03-10 08:21:40 +00:00
parent 4f919bc007 9afd40a1ad
commit 58a355478a
9 changed files with 982 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
examples/Makefile.am
View File

@@ -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

28
examples/exDWT.cpp Normal file
View File

@@ -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;
}

4
lib/Makefile.am
View File

@@ -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\

137
lib/Numerical/DWT.cpp Normal file
View File

@@ -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;
}

55
lib/Numerical/DWT.hpp Normal file
View File

@@ -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_

528
lib/Numerical/DWTFilters.cpp Normal file
View File

@@ -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}
};

53
lib/Numerical/DWTFilters.hpp Normal file
View File

@@ -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_

5
utils/Makefile.am
View File

@@ -9,6 +9,7 @@ endif
bin_PROGRAMS = \
latan-plot \
latan-sample-combine \
latan-sample-dwt \
latan-sample-element \
latan-sample-fake \
latan-sample-ft \
@@ -26,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

167
utils/sample-dwt.cpp Normal file
View File

@@ -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;
}