【小波变换】STL版 一维离散小波变换（DWT)库,完全按matlab的wavelet toolbox 的API实现的
【转】一维离散小波变换（DWT)库,完全按matlab的wavelet toolbox 的API实现的
来源：http://hi.baidu.com/anatacollin/item/69fdab74ca7d045c0d0a07b4
一维离散小波变换（DWT)库,完全按matlab的wavelet toolbox 的API实现的2008-12-01 20:37最近项目中需要用，就自己写了个，发在这里算是备忘。需要的朋友也可以拿去试试，经测试没有发现bug，基于STL实现。如果发现bug或有什么建议请通知我，谢谢。
/************************************************************************/
/* wavelet.h
* Author: Collin
* Date: 2008/12/01
*/
/************************************************************************/
#ifndef WAVELET_H
#define WAVELET_H
#include <vector>
namespace Wavelet{
using std::vector;
struct C_L
{
vector<double> C;
vector<int> L;
};
struct WaveFilter{
vector<double> Low;
vector<double> High;
};
struct WaveCoeff{
vector<double> app;
vector<double> det;
};
const double sym4_Lo_D[] = {-0.0758, -0.0296, 0.4976, 0.8037, 0.2979, -0.0992, -0.0126, 0.0322};
const double sym4_Hi_D[] = {-0.0322, -0.0126, 0.0992, 0.2979, -0.8037, 0.4976, 0.0296, -0.0758};
const double sym4_Lo_R[] = {0.0322, -0.0126, -0.0992, 0.2979, 0.8037, 0.4976, -0.0296, -0.0758};
const double sym4_Hi_R[] = {-0.0758, 0.0296, 0.4976, -0.8037, 0.2979, 0.0992, -0.0126, -0.0322};
const static WaveFilter sym4_d = {vector<double>(sym4_Lo_D, sym4_Lo_D + 8), vector<double>(sym4_Hi_D, sym4_Hi_D + 8)};
const static WaveFilter sym4_r = {vector<double>(sym4_Lo_R, sym4_Lo_R + 8), vector<double>(sym4_Hi_R, sym4_Hi_R + 8)};
const WaveFilter& WFilters(
const char* strWaveName,
const char d_or_r
);
C_L WaveDec(
const vector<double>& signal,
const int nMaxLevel,
const char* strWaveName
);
WaveCoeff DWT(
const vector<double>& signal,
const vector<double>& Lo_D,
const vector<double>& Hi_D
);
vector<double> WRCoef(
const char a_or_d,
const vector<double>& C,
const vector<int>& L,
const char* strWaveName,
const int nLevel
);
vector<double> AppCoef(
const vector<double>& C,
const vector<int>& L,
const char* strWaveName,
const int nLevel
);
vector<double> DetCoef(
const vector<double>& C,
const vector<int>& L,
const int nLevel
);
//upsample and convolution
vector<double> UpsConv1(
const vector<double>& signal,
const vector<double>& filter,
const int nLen,
const char* strMode = "sym"
);
vector<double> Conv(
const vector<double>& vecSignal,
const vector<double>& vecFilter
);
vector<double> IDWT(
const vector<double>& app,
const vector<double>& det,
const vector<double>& Lo_R,
const vector<double>& Hi_R,
const int nLenCentral
);
vector<double> WExtend(
const vector<double>& signal,
const int nLenExt,
const char* mode = "sym"
);
vector<double> WConv1(
const vector<double>& signal,
const vector<double>& filter,
const char* shape = "valid"
);
}
#endif
/************************************************************************/
/* wavelet.cpp
* Author: Collin
* Date: 2008/12/01
*/
/************************************************************************/
#include <vector>
#include <string>
#include <iostream>
#include "wavelet.h"
using namespace std;
using namespace Wavelet;
C_L Wavelet::WaveDec(const vector<double>& signal,
const int nMaxLevel,
const char* strWaveName
)
{
const WaveFilter& filters = WFilters(strWaveName, 'd');
int len = signal.size();
C_L cl;
cl.L.push_back(len);
WaveCoeff waveCoeff;
waveCoeff.app = signal;
vector<double>::iterator itC;
vector<int>::iterator itL;
for (int i = 0; i < nMaxLevel; ++i){
waveCoeff = DWT(waveCoeff.app, filters.Low, filters.High);
itC = cl.C.begin();
cl.C.insert(itC, waveCoeff.det.begin(), waveCoeff.det.end());
itL = cl.L.begin();
cl.L.insert(itL, waveCoeff.det.size());
}
itC = cl.C.begin();
cl.C.insert(itC, waveCoeff.app.begin(), waveCoeff.app.end());
itL = cl.L.begin();
cl.L.insert(itL, waveCoeff.app.size());
return cl;
}
vector<double> Wavelet::WRCoef(const char a_or_d,
const vector<double>& C,
const vector<int>& L,
const char* strWaveName,
const int nLevel
)
{
vector<double> Coef;
const WaveFilter& filter = WFilters(strWaveName, 'r');
int nMax = L.size() - 2;
int nMin;
char type = tolower(a_or_d);
if ('a' == type)
nMin = 0;
else if ('d' == type)
nMin = 1;
else {
cerr << "bad parameter: a_or_d: "<< a_or_d << "\n";
exit(1);
}
if (nLevel < nMin || nLevel > nMax){
cerr << "bad parameter for level\n";
exit(1);
}
vector<double> F1;
switch (type){
case 'a':
Coef = AppCoef(C, L, strWaveName, nLevel);
if (0 == nLevel)
return Coef;
F1 = filter.Low;
break;
case 'd':
Coef = DetCoef(C, L, nLevel);
F1 = filter.High;
break;
default:
;
}
int iMin = L.size() - nLevel;
Coef = UpsConv1(Coef, F1, L[iMin], "sym");
for (int k = 1; k < nLevel; ++k){
Coef = UpsConv1(Coef, filter.Low, L[iMin + k], "sym");
}
return Coef;
}
vector<double> Wavelet::UpsConv1(const vector<double>& signal,
const vector<double>& filter,
const int nLen,
const char* strMode
)
{
//implement dyadup(y,0)
vector<double> y(2 * signal.size() - 1);
y[0] = signal[0];
for (int i = 1; i < signal.size(); ++i){
y[2*i - 1] = 0;
y[2*i] = signal[i];
}
y = Conv(y, filter);
//extract the central portion
vector<double>::iterator it = y.begin();
return vector<double>(it + (y.size() - nLen) / 2, it + (y.size() + nLen) / 2);
}
vector<double> Wavelet::Conv(const vector<double>& vecSignal, const vector<double>& vecFilter){
vector<double> signal(vecSignal);
vector<double> filter(vecFilter);
if (signal.size() < filter.size())
signal.swap(filter);
int lenSignal = signal.size();
int lenFilter = filter.size();
vector<double> result(lenSignal + lenFilter - 1);
for (int i = 0; i < lenFilter; i++){
for (int j = 0; j <= i; j++)
result[i] += signal[j] * filter[i - j];
}
for (int i = lenFilter; i < lenSignal; i++){
for (int j = 0; j <lenFilter; j++)
result[i] += signal[i - j] * filter[j];
}
for (int i = lenSignal; i < lenSignal + lenFilter - 1; i++){
for (int j = i - lenSignal + 1; j < lenFilter; j++)
result[i] += signal[i - j] * filter[j];
}
return result;
}
vector<double> Wavelet::DetCoef(const vector<double>& C,
const vector<int>& L,
const int nLevel
)
{
if (nLevel < 1 || nLevel > L.size() - 2){
cerr << "bad level parameter\n";
exit(1);
}
int nlast = 0, nfirst = 0;
vector<int>::const_reverse_iterator it = L.rbegin();
++it;
for (int i = 1; i < nLevel; ++i){
nlast += *it;
++it;
}
nfirst = nlast + *it;
return vector<double>(C.end() - nfirst, C.end() - nlast);
}
WaveCoeff Wavelet::DWT(const vector<double>& signal,
const vector<double>& Lo_D,
const vector<double>& Hi_D
)
{
int nLenExt = Lo_D.size() - 1;
vector<double> y;
y = WExtend(signal, nLenExt, "sym");
vector<double> z;
z = WConv1(y, Lo_D, "valid");
WaveCoeff coeff;
for (int i = 1; i < z.size(); i += 2){
coeff.app.push_back(z[i]);
}
z = WConv1(y, Hi_D, "valid");
for (int i = 1; i < z.size(); i += 2){
coeff.det.push_back(z[i]);
}
return coeff;
}
const WaveFilter& Wavelet::WFilters(const char* strWaveName,
const char d_or_r
)
{
char type = tolower(d_or_r);
if (!strcmp(strWaveName, "sym4")){
switch(type){
case 'd':
return Wavelet::sym4_d;
break;
case 'r':
return Wavelet::sym4_r;
break;
default:
cerr << "bad parameter for d_or_r\n";
exit(1);
}
}
else {
cerr << "not implement \n";
exit(1);
}
}
vector<double> Wavelet::AppCoef(
const vector<double>& C,
const vector<int>& L,
const char* strWaveName,
const int nLevel
)
{
int nMaxLevel = L.size() - 2;
if (nLevel < 0 || nLevel > nMaxLevel){
cerr << "bad parameter for level\n";
exit(1);
}
const WaveFilter& filters = WFilters(strWaveName, 'r');
vector<double> app(C.begin(), C.begin() + L[0]); //app for the last level
vector<double> det;
for (int i = 0; i < nMaxLevel - nLevel; ++i){
det = DetCoef(C, L, nMaxLevel - i);
app = IDWT(app, det, filters.Low, filters.High, L[i + 2]);
}
return app;
}
vector<double> Wavelet::IDWT(
const vector<double>& app,
const vector<double>& det,
const vector<double>& Lo_R,
const vector<double>& Hi_R,
const int nLenCentral
)
{
vector<double> app1, app2;
app1 = UpsConv1(app, Lo_R, nLenCentral, "sym");
app2 = UpsConv1(det, Hi_R, nLenCentral, "sym");
for (int i = 0; i < nLenCentral; ++i){
app1[i] += app2[i];
}
return app1;
}
vector<double> Wavelet::WExtend(
const vector<double>& signal,
const int nLenExt,
const char* mode
)
{
int signalLen = signal.size();
vector<double> result(signalLen + 2 * nLenExt);
for (int i = 0, idx = nLenExt; idx < signalLen + nLenExt; ++i, ++idx){
result[idx] = signal[i];
}
for (int idx = nLenExt - 1, bFlag = 1, signalIdx = 0; idx >= 0; --idx){
result[idx] = signal[signalIdx];
if (bFlag && ++signalIdx == signalLen){
bFlag = 0;
signalIdx = signalLen - 1;
}
else if (!bFlag && --signalIdx == -1) {
bFlag = 1;
signalIdx = 0;
}
}
for (int idx = nLenExt + signalLen, bFlag = 0, signalIdx = signalLen - 1; idx < 2 * nLenExt + signalLen; ++idx){
result[idx] = signal[signalIdx];
if (bFlag && ++signalIdx == signalLen){
bFlag = 0;
signalIdx = signalLen - 1;
}
else if (!bFlag && --signalIdx == -1) {
bFlag = 1;
signalIdx = 0;
}
}
return result;
}
vector<double> Wavelet::WConv1(
const vector<double>& signal,
const vector<double>& filter,
const char* shape
)
{
vector<double> y;
y = Conv(signal, filter);
int nLenExt = filter.size() - 1;
return vector<double>(y.begin() + nLenExt, y.end() - nLenExt);
}
