mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2024-11-10 00:45:36 +00:00
138 lines
4.1 KiB
C++
138 lines
4.1 KiB
C++
/*
|
|
* 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;
|
|
}
|