From 13fddf4947c18f2981d0a1a0f9025649297cc986 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Fri, 12 Jan 2024 14:22:05 +0100
Subject: [PATCH] DWT concatenation and CDR

---
 lib/Numerical/DWT.cpp | 23 ++++++++++++++++++++++
 lib/Numerical/DWT.hpp |  2 ++
 utils/sample-dwt.cpp  | 44 +++++++++++++++++++++++++++++++++++++------
 3 files changed, 63 insertions(+), 6 deletions(-)

diff --git a/lib/Numerical/DWT.cpp b/lib/Numerical/DWT.cpp
index bc0e03a..8cf4f6b 100644
--- a/lib/Numerical/DWT.cpp
+++ b/lib/Numerical/DWT.cpp
@@ -135,3 +135,26 @@ DVec DWT::backward(const std::vector<DWTLevel>& dwt) const
 
     return res;
 }
+
+// concatenate levels //////////////////////////////////////////////////////////
+DVec DWT::concat(const std::vector<DWTLevel> &dwt, const int maxLevel, const bool dropLow)
+{
+    unsigned int level = ((maxLevel >= 0) ? (maxLevel + 1) : dwt.size());
+    Index nlast = dwt[level - 1].first.size();
+    Index n = 2*dwt.front().first.size() - ((dropLow) ? nlast : 0);
+    Index pt = n, nl;
+    DVec  res(n);
+
+    for (unsigned int l = 0; l < level; ++l)
+    {
+        nl  = dwt[l].second.size();
+        pt -= nl;
+        res.segment(pt, nl) = dwt[l].second;
+    }
+    if (!dropLow)
+    {
+        res.segment(0, nl) = dwt[level-1].first;
+    }
+
+    return res;
+}
diff --git a/lib/Numerical/DWT.hpp b/lib/Numerical/DWT.hpp
index 5745245..ebc512f 100644
--- a/lib/Numerical/DWT.hpp
+++ b/lib/Numerical/DWT.hpp
@@ -46,6 +46,8 @@ public:
     // DWT
     std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
     DVec                  backward(const std::vector<DWTLevel>& dwt) const;
+    // concatenate levels
+    static DVec concat(const std::vector<DWTLevel>& dwt, const int maxLevel = -1, const bool dropLow = false);
 private:
     DWTFilter filter_;
 };
diff --git a/utils/sample-dwt.cpp b/utils/sample-dwt.cpp
index 958eea6..522e692 100644
--- a/utils/sample-dwt.cpp
+++ b/utils/sample-dwt.cpp
@@ -17,6 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/Core/OptParser.hpp>
 #include <LatAnalyze/Core/Plot.hpp>
 #include <LatAnalyze/Io/Io.hpp>
@@ -53,7 +54,13 @@ int main(int argc, char *argv[])
         cerr << "usage: " << argv[0];
         cerr << " <options> <input file>" << endl;
         cerr << endl << "Possible options:" << endl << opt << endl;
-        
+        cerr << "Available DWT filters:" << endl;
+        for (auto &fv: DWTFilters::fromName)
+        {
+            cerr << fv.first << " ";
+        }
+        cerr << endl << endl;
+
         return EXIT_FAILURE;
     }
     inFilename  = opt.getArgs()[0];
@@ -68,22 +75,45 @@ int main(int argc, char *argv[])
     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));
+                          outh(ss ? 0 : level, DMatSample(nSample)),
+                          concath(ss ? 0 : level, DMatSample(nSample));
+    DMatSample            concat(nSample, n, 1);
     DWT                   dwt(*DWTFilters::fromName.at(filterName));
     vector<DWT::DWTLevel> dataDWT(level);
 
+    FOR_STAT_ARRAY(in, s)
+    {
+        in[s].conservativeResize(n, 1);
+    }
     if (!ss)
     {
+        DMatSample buf(nSample);
+
         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);
+            dataDWT = dwt.forward(in[s], level);
             for (unsigned int l = 0; l < level; ++l)
             {
-                out[l][s]  = dataDWT[l].first;
-                outh[l][s] = dataDWT[l].second;
+                out[l][s]     = dataDWT[l].first;
+                outh[l][s]    = dataDWT[l].second;
+                concath[l][s] = DWT::concat(dataDWT, l, true);
             }
+            concat[s] = DWT::concat(dataDWT);
+        }
+        cout << "Data CDR " << Math::cdr(in.correlationMatrix()) << " dB" << endl;
+        cout << "DWT CDR " << Math::cdr(concat.correlationMatrix()) << " dB" << endl;
+        for (unsigned int l = 0; l < level; ++l)
+        {
+            cout << "DWT level " << l << " CDR: L= ";
+            cout << Math::cdr(out[l].correlationMatrix()) << " dB / H= ";
+            cout << Math::cdr(outh[l].correlationMatrix()) << " dB" << endl;
+        }
+        for (unsigned int l = 0; l < level; ++l)
+        {
+            cout << "DWT detail level " << l << " CDR: ";
+            cout << Math::cdr(concath[l].correlationMatrix()) << " dB" << endl;
         }
     }
     else
@@ -102,7 +132,7 @@ int main(int argc, char *argv[])
         }
         FOR_STAT_ARRAY(in, s)
         {
-            dataDWT.back().first = in[s].col(0);
+            dataDWT.back().first = in[s];
             out[0][s] = dwt.backward(dataDWT);
         }
     }
@@ -115,7 +145,9 @@ int main(int argc, char *argv[])
             {
                 Io::save<DMatSample>(out[l], outFilename + "/L" + strFrom(l) + ".h5");
                 Io::save<DMatSample>(outh[l], outFilename + "/H" + strFrom(l) + ".h5");
+                Io::save<DMatSample>(concath[l], outFilename + "/concatH" + strFrom(l) + ".h5");
             }
+            Io::save<DMatSample>(concat, outFilename + "/concat.h5");
         }
         else
         {