Not achieving perfect reconstruction with symmlet

[Ryan_Mitchley]

I am struggling to achieve perfect (wavelet) reconstruction usingippsWTFwd_32f andippsWTInv_32f, given an order 3, length 8 symmlet. I have calculated the group delay of the forward "low" filter as 4 and the group delay of the forward "high" filter as 3. However, I can't seem to get perfect reconstruction with any offset settings.

The docs note that "Biorthogonal and orthogonal filter banks are distinguished by one specific

peculiarity, that is, inverse transform additional delays must be uniformly even and

opposite to the evenness of the decomposition delays for faultless signal reconstruction."

However, I can't seem to get reconstruction working, using even delays.

I've attached a code sample I made to demonstrate the problem, largely based on example code:

#include

#include

#include

#include

#include

#include

#include

// Filter bank for symmlet, order 3, filter

static const int fwdFltLenL = 8;

static const int fwdFltLenH = 8;

static const float fwdFltL[8] =

{ // group delay 4

0.03222310060407f,

-0.01260396726226f,

-0.09921954357694f,

0.29785779560554f,

0.80373875180522f,

0.49761866763246f,

-0.02963552764595f,

-0.07576571478934f

};

static const float fwdFltH[8] =

{ // group delay 3

0.07576571478934f,

-0.02963552764595f,

-0.49761866763246f,

0.80373875180522f,

-0.29785779560554f,

-0.09921954357694f,

0.01260396726226f,

0.03222310060407f

};

static const int invFltLenL = 8;

static const int invFltLenH = 8;

static const float invFltL[8] =

{

-0.07576571478934f,

-0.02963552764595f,

0.49761866763246f,

0.80373875180522f,

0.29785779560554f,

-0.09921954357694f,

-0.01260396726226f,

0.03222310060407f

};

static const float invFltH[8] =

{

0.03222310060407f,

0.01260396726226f,

-0.09921954357694f,

-0.29785779560554f,

0.80373875180522f,

-0.49761866763246f,

-0.02963552764595f,

0.07576571478934f

};

int filterlength = 8;

using namespace std;

/// Number of decomposition levels

int iNumLevels;

/// Number of input samples

int N = 32;

/// Count of samples used per level

int iCount;

/// Forward transform states for each level

IppsWTFwdState_32f * vpFwdState;

/// Inverse transform states for each level

IppsWTInvState_32f * vpInvState;

/// Vector of high (difference) decompositions (one per level)

float * vpfDecompLow;

/// Vector of low (average) decompositions (one per level)

float * vpfDecompHigh;

/// Vector to hold the reconstructed data (exceeds iN by extraHead samples)

vector vfReconstruction;

Ipp32f * vpfFwdDelayLow;

Ipp32f * vpfFwdDelayHigh;

Ipp32f * vpfInvDelayLow;

Ipp32f * vpfInvDelayHigh;

int viFwdOffsetLow;

int viFwdOffsetHigh;

int viInvOffsetLow;

int viInvOffsetHigh;

int main(int argc, char ** argv)

{

int i;

vfReconstruction.resize(N + filterlength - 1);

iCount = (N + filterlength - 1) / 2;

viFwdOffsetLow = -1; viFwdOffsetHigh = 0;

viInvOffsetLow = 0; viInvOffsetHigh = -1;

ippsWTFwdInitAlloc_32f (&vpFwdState, fwdFltL, fwdFltLenL, viFwdOffsetLow, fwdFltH, fwdFltLenH, viFwdOffsetHigh);

ippsWTInvInitAlloc_32f (&vpInvState, invFltL, invFltLenL, viInvOffsetLow, invFltH, invFltLenH, viInvOffsetHigh);

vpfDecompLow = new float[iCount];

vpfDecompHigh = new float[iCount];

// Set up delay lines for forward transform

vpfFwdDelayLow = new float[fwdFltLenL + viFwdOffsetLow - 1];

vpfFwdDelayHigh = new float[fwdFltLenH + viFwdOffsetHigh - 1];

memset(vpfFwdDelayLow, 0, (fwdFltLenL + viFwdOffsetLow - 1) * sizeof(float));

memset(vpfFwdDelayHigh, 0, (fwdFltLenH + viFwdOffsetHigh - 1) * sizeof(float));

ippsWTFwdSetDlyLine_32f(vpFwdState, vpfFwdDelayLow, vpfFwdDelayHigh);

// Set up delay lines for inverse transform

vpfInvDelayLow = new float[(invFltLenL + viInvOffsetLow - 1) / 2];

vpfInvDelayHigh = new float[(invFltLenH + viInvOffsetHigh - 1) / 2];

memset(vpfInvDelayLow, 0, ((invFltLenL + viInvOffsetLow - 1) / 2) * sizeof(float));

memset(vpfInvDelayHigh, 0, ((invFltLenH + viInvOffsetHigh - 1) / 2) * sizeof(float));

ippsWTInvSetDlyLine_32f(vpInvState, vpfInvDelayLow, vpfInvDelayHigh);

Ipp32f *src = new Ipp32f;

/* for (i = 0; i < 32; i++)

src = int(10000.0f * float(rand()) / float(RAND_MAX) - 5000.0f);

*/

src[0] = 1.0;

src[1] = 0.0;

printf("\\noriginal:\\t");

for (i = 0; i < N; i++)

printf("%.4f, ", src);

// Analyse

ippsWTFwd_32f (src, vpfDecompLow, vpfDecompHigh, iCount, vpFwdState);

// Synthesize

int extraHead = 4;

ippsWTInv_32f (vpfDecompLow, vpfDecompHigh+extraHead, iCount, (Ipp32f *) &vfReconstruction[0], vpInvState);

printf("\\nreconstruction:\\t");

for (i = 0; i < N; i++)

printf("%.4f, ", vfReconstruction);

printf(" \\n");

ippsWTFwdFree_32f (vpFwdState);

delete vpfDecompLow;

delete vpfDecompHigh;

delete vpfFwdDelayLow;

delete vpfFwdDelayHigh;

ippsWTInvFree_32f (vpInvState);

delete vpfInvDelayLow;

delete vpfInvDelayHigh;

return 0;

}

[Ryan_Mitchley]

*Bump* ANY ideas?

[Chao_Y_Intel]

Hello,

The border (delay line) extension does not look to set correctly. These data set as 0 in the code.

memset(vpfFwdDelayLow, 0, (fwdFltLenL + viFwdOffsetLow - 1) * sizeof(float));
memset(vpfFwdDelayHigh, 0, (fwdFltLenH + viFwdOffsetHigh - 1) * sizeof(float));

I attached one sample code, which were used before. The code used wrapped border extension. Please check if it can provide some help.

Thanks,
Chao