本文介绍了一种基于二次样条小波的二进小波变换及其逆变换的实现方法，包括前向变换和逆向变换的具体算法流程，并提供了一个测试实例来验证算法的有效性。
头文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program 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 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program 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. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* bwt.h
*
* Dyadic Wavelet Transform.
*
* These routines are designed for computing the dyadic wavelet transform
* and it's inverse transform using quadratic spline wavelet.
*
* To distinguish with the "dwt" (discrete wavelet transform), we call this
* file as "bwt", but in fact, it should be dyadic wavelet transform. *
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef BWT_H
#define BWT_H
#include <vector.h>
#include <utilities.h>
namespace splab
{
template<typename Type> Vector< Vector<Type> > bwt( const Vector<Type>&,
int );
template<typename Type> Vector<Type> ibwt( const Vector< Vector<Type> >&,
int );
#include <bwt-impl.h>
}
// namespace splab
#endif
// BWT_H实现文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program 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 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program 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. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* bwt-impl.h
*
* Implementation for Dyadic Wavelet Transform.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
/**
* Forward Transform.
* The decomposition levels is specified by integer "J". The decomposed
* coefficients are stroed in a "Vector< vector<Type> >" structure.
* Detial coefficients are stored from 1st to Jth row, and approximation
* coefficients are stored at the last row, i.e. the (J+1)th row.
*/
template <typename Type>
Vector< Vector<Type> > bwt( const Vector<Type> &xn, int J )
{
// lowpass decomposition filter
Vector<Type> ld(4);
ld[0] = 0.125; ld[1] = 0.375;
ld[2] = 0.375; ld[3] = 0.125;
ld = Type(RT2) * ld;
int ldZeroStart = 1;
// highpass decomposition filter
Vector<Type> hd(2);
hd[0] = -0.5; hd[1] = 0.5;
hd = Type(RT2) * hd;
int hdZeroStart = 0;
// initializing the coefficients
int N = xn.size();
Vector< Vector<Type> > coefs(J+1);
for( int i=0; i<J; ++i )
coefs[i].resize(N);
Vector<Type> approx(N);
Vector<Type> a(xn);
// get the inversion of filters
Vector<Type> ll = flip(ld);
Vector<Type> hh = flip(hd);
int llZeroStart = 0,
hhZeroStart = 0,
p = 1;
for( int j=0; j<J; ++j )
{
// compute the 0 position of the new filters
llZeroStart = ll.size()-1 - p*ldZeroStart;
hhZeroStart = hh.size()-1 - p*hdZeroStart;
for( int i=0; i<N; ++i )
{
Type sum = 0;
// compute the approximation coefficients
for( int k=0; k<ll.size(); k+=p )
{
int index = mod( i+llZeroStart-k, N );
sum += ll[k]*a[index];
}
approx[i] = sum;
// compute the detial coefficients
sum = 0;
for( int k=0; k<hh.size(); k+=p )
{
int index = mod( i+hhZeroStart-k, N );
sum += hh[k]*a[index];
}
coefs[j][i] = sum;
}
a = approx;
// dyadic upsampling
ll = dyadUp( ll, 1 );
hh = dyadUp( hh, 1 );
p *= 2;
}
coefs[J] = approx;
return coefs;
}
/**
* Backword Transform.
* The reconstruction livel is specified by integer "level", and
* "levle" should between "0"(the approximation component) and "J"
* (the original signal).
*/
template <typename Type>
Vector<Type> ibwt( const Vector< Vector<Type> > &coefs, int level )
{
Vector<Type> lr(4);
lr[0] = 0.125; lr[1] = 0.375;
lr[2] = 0.375; lr[3] = 0.125;
lr = Type(RT2) * lr;
int lrZeroStart = 1;
Vector<Type> hr(6);
hr[0] = -0.03125; hr[1] = -0.21875; hr[2] = -0.6875;
hr[3] = 0.6875; hr[4] = 0.21875; hr[5] = 0.03125;
hr = Type(RT2) * hr;
int hrZeroStart = 2;
int J = coefs.dim() - 1;
if( (level < 0) || (level > J) )
{
cout << "invalid reconstruction level!" << endl;
return Vector<Type>(0);
}
int N = coefs[0].dim();
Vector<Type> a = coefs[J];
Vector<Type> xn(N);
Vector<Type> ll = lr;
Vector<Type> hh = hr;
int llZeroStart = 0;
int hhZeroStart = 0;
int p = 1;
// get the Jth level filters
for( int j=0; j<level-1; ++j )
{
p *= 2;
ll = dyadUp( ll, 1 );
hh = dyadUp( hh, 1 );
}
for( int j=level-1; j>=0; --j )
{
// compute the 0 position of the new filters
llZeroStart = p*lrZeroStart;
hhZeroStart = p*hrZeroStart;
// compute the jth approximation coefficients
for( int i=0; i<N; ++i )
{
Type sum = 0;
for( int k=0; k<ll.size(); k+=p )
{
int index = mod( i+llZeroStart-k, N );
sum += ll[k]*a[index];
}
for( int k=0; k<hh.size(); k+=p )
{
int index = mod( i+hhZeroStart-k, N );
sum += hh[k]*coefs[j][index];
}
xn[i] = sum/2;
}
a = xn;
// dyadic downsampling
ll = dyadDown( ll, 0 );
hh = dyadDown( hh, 0 );
p /= 2;
}
return xn;
}测试代码：
/*****************************************************************************
* bwt_test.cpp
*
* Dyadic wavelet transform testing.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <bwt.h>
#include <timing.h>
using namespace std;
using namespace splab;
const int Ls = 100;
int main()
{
/************************** [ signal ] *************************/
Vector<double> s(Ls);
for(int i=0; i<Ls; i++)
{
if(i<Ls/4)
s[i] = 0.0;
else if(i<2*Ls/4)
s[i] = 1.0;
else if(i<3*Ls/4)
s[i] = 3.0;
else
s[i] = 0.0;
}
/*************************** [ BWT ] ***************************/
int level = 3;
Timing cnt;
double runtime = 0.0;
cout << "Taking dyadic wavelet transform." << endl;
cnt.start();
Vector< Vector<double> > coefs = bwt( s, level );
cnt.stop();
runtime = cnt.read();
cout << "The running time = " << runtime << " (s)" << endl << endl;
/*************************** [ IBWT ] **************************/
cout << "Taking inverse dyadic wavelet transform." << endl;
cnt.start();
Vector<double> x = ibwt( coefs, level );
cnt.stop();
runtime = cnt.read();
cout << "The running time = " << runtime << " (s)" << endl << endl;
cout << "The relative error is : norm(s-x) / norm(s) = "
<< norm(s-x)/norm(s) << endl << endl;
return 0;
}运行结果：
Taking dyadic wavelet transform.
The running time = 0 (s)
Taking inverse dyadic wavelet transform.
The running time = 0 (s)
The relative error is : norm(s-x) / norm(s) = 3.76025e-016
Process returned 0 (0x0) execution time : 0.062 s
Press any key to continue.
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
