Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Ryan_Mitchley
Beginner
65 Views

Not achieving perfect reconstruction with symmlet

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;
}
0 Kudos
2 Replies
Ryan_Mitchley
Beginner
65 Views

*Bump* ANY ideas?
Chao_Y_Intel
Employee
65 Views


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

Reply