博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
【小波变换】STL版 一维离散小波变换(DWT)库,完全按matlab的wavelet toolbox 的API实现的...
阅读量:5058 次
发布时间:2019-06-12

本文共 8959 字,大约阅读时间需要 29 分钟。

【转】一维离散小波变换(DWT)库,完全按matlab的wavelet toolbox 的API实现的

来源:

一维离散小波变换(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);
}

转载于:https://www.cnblogs.com/IDoIUnderstand/archive/2013/04/02/3280723.html

你可能感兴趣的文章
Linux服务器删除乱码文件/文件夹的方法
查看>>
牛腩记账本core版本源码
查看>>
Word Break II
查看>>
UVA 11082 Matrix Decompressing 矩阵解压(最大流,经典)
查看>>
jdk从1.8降到jdk1.7失败
查看>>
一些关于IO流的问题
查看>>
mongo备份操作
查看>>
8 -- 深入使用Spring -- 3...1 Resource实现类InputStreamResource、ByteArrayResource
查看>>
硬件笔记之Thinkpad T470P更换2K屏幕
查看>>
一个关于vue+mysql+express的全栈项目(六)------ 聊天模型的设计
查看>>
【知识库】-数据库_MySQL 的七种 join
查看>>
.net 写文件上传下载webservice
查看>>
noSQL数据库相关软件介绍(大数据存储时候,必须使用)
查看>>
iOS开发——缩放图片
查看>>
HTTP之URL的快捷方式
查看>>
满世界都是图论
查看>>
配置链路聚合中极小错误——失之毫厘谬以千里
查看>>
代码整洁
查看>>
蓝桥杯-分小组-java
查看>>
Java基础--面向对象编程1(类与对象)
查看>>