1-d Wavelet functions.

[jimmyliu51]

I am studying the wavelet functions in IPP. Two questions:

1. How to use delay line set up functions? Is there any examples? I have tried to find one but cannot find anything yet.

2. I have some example in Wavelab, the samples are padded, I like to know how to do it in IPP. Because padding can be done either or both end of

the signal, I don't know if delay line set up function in IPP will do the same or not. At least I don't see a way to specify the padding method.

Thank you!

Jun

[jimmyliu51]

I haven't seen any reply regarding my previous posting yet. I also found a section

in IPP manual discussing boundaries extrapolation on page 7-68. I copy it here

and hopefully someone can help me to figure out:

"Typically, reversible wavelet transforms of a bounded signal require data extrapolation towards one or both sides. All internal delay lines are set to zero at the stage of initialization. To set a non-zero signal prehistory, call the function ippsWTFwdSetDlyLine. When processed an entire limited data set, data extrapolation may be performed both towards the start and the end of the data vector. For that, the source data and their initial extrapolation are used to form the delay line, the rest of the signal is subdivided into the main block and the signal end. The signal end data and their extrapolation are used to form the last block."

The last long sentence is very strange to me becauseit is not clear to me where the main block and last block should be put.

Jun

[hrvoje_intel]

I can not say that I understand the idea completely but I have managed to make perfect reconstruction.

Here is how it goes.

First of all it is the same procedure for both low pass and high pass filtering as well as for the analysis and the reconstruction phase.
One part of wavelet decomposition is FIR filtering of the signal. The problem when performing the FIR filtering is what to do with the boundaries because clasic zero padding will not work (at least not if you want to have exactly tha half the number of coefficients after one decomposition step as implemented in IPP).

The procedure with the delay lines is as follows:
You need to put some bondary extrapolated values and shift some data values into the delay line. Since you shifted some data values into the delay line you need to fill the now missing signal values at the end. You fill it again with the extrapolated boundary values.



Example:

Signal [1 2 3 4 5 6 7 8], you have 4-coeff FIR filter and using periodic boundary extension.

Your delay line is 3 elements long since your FIR filter is positioned on the first element of the signal so you need 3 boundary values. Filtering will be performed in 8 subsequent steps. Now you shift the signal by one position to the left into the delay line. The reason for this is that you want to simetrically spread the boundary influence on both sides of the signal not just the front one. For longer FIR filters larger shift is required.

So the delay line is [6 7 1] and now fill this the missing value of the signal with boundary extension value. Now you have the delay line [6 7 1] then the signal [2 3 4 5 6 7 8 1]. Use these two in your IPP wavelet function.



Regards,
Hrvoje

[jimmyliu51]

Thank you for your advice!

I will try it again. But, when you talk about delay line, do you mean both high pass and low pass delay lines? The reason I cannot figure out this is that the matlab code I have use different way for high pass and low pass filtering.

IPP seems don't care where the delay line goes, or just pad at the end of the signal.

Also, in your example, why delay line is not [7 8 1]?

Since you move 1 to the back, I thought 8 should be in the front of 1.

Jun

[hrvoje_intel]

You are right about the delay line values, I made a mistake, they should be [7 8 1].

About your question with high and low pass filters.
Yes, I mean for both, although in my case I used only the orthogonal filters so low and high pass filters were equally long. I do not know how your code differs between low and high pass filtering and I am not sure whether this will require different procedure than mine.

Regards,
Hrvoje

[Mikhail_Kulikov__Int]

Hi,

I just decided to place some more or less simple example right here.

I started from'db2' filters from MATLAB wavelet toolbox. It can be usefull if you have MATLAB and already know what to do with it.

But if you do not have MATLAB you can start from "main.cpp" code. It's simple enough.

Before starting the example,
some words about interface models in MATLAB and IPP.
And regarding use of these tools together.

MATLAB wavelet toolbox interface model (aka "dwt" function interface) assumes applying of some of fixed enumerated boundaries extension models. To define it explicity you need to pass 'zpd' (zero-padding extension), 'sym' (symmetrical extension) or some other pre-defined model specifiers to "dwt" or "idwt" function (or globally to "dwtmode" function). It works excellent until you work with non-continues signal. If you work with streaming continuous signal (stream from ADC, for example) you need to use "overlapping" bocks to work-around pre-defined boundaries. Of cousre MATLAB wavelet toolbox is not targeted to be the high-performance tool for low-latency real-time streamming processing applications, it's excellent scintific educational and research tool and it contains a LOT of usefull information.

IPP signal processing functions interface model is intended for streaming continuous signal processing (at least), and IPP wavelet transform are not exclusion. This kind of interfaces is based on delay-line paradigm (tendentiously saying "to work-around" boundary extensions). It is intended mostly for processing of continuous signal, at the same time it is flexible enough (but low-level) to adopt different wavelet filters and arbitrary border extensions as well. But you need to care about some details, like filter values, 'anchor' positions for filters, synchronizing of low-pass and high-pass component delays, that you do not care about in MATLAB just specifying "db2" or other text value to chose "default" helper structures (and it's hard to change some of them for example if you really need to refuse from non-centered filter anchors and expanding length of decomposition).
IPP model gives the way to implement effectively custom WT design that hard to implement in MATLAB, at least without performance loss or additional programming. At the same time you need to specify all details about your filters and implement border model on both edges, if you work with non-streamming data (for streamming data zero-pre-history is implemented implicitly and usually that's enough).
And you can get a lot of data about filter taps, best extension model from MATLAB to use with IPP on the stage of design,
and improve and adjust some details using IPP implementation.

The following example illustrates correspondence of MATLAB and IPP,
together with their co-operative usage.

I. Let's start with MATLAB dwt short example:

>> dwtmode 'ppd';
>> s = [1 -10 324 48 -483 4 7 -5532 34 8889 -57 54];
>> [approx,detail] = dwt(s,'db2')

The output is:

approx =

Columns 1 through 6

1.9161e+001 5.8529e+001 8.7854e+001 4.8754e+002 -5.7669e+003 7.4324e+003

Column 7

1.9161e+001

detail =

Columns 1 through 6

9.3872e-001 2.4996e+0 02 -4.5866e+002 2.7392e+003 -3.0256e+003 -2.0706e+003

Column 7

9.3872e-001

To obtain perfect reconstruction we just tape:
>> idwt(approx, detail, 'db2')

we ommit output here, because it repeatsoriginal numbers (with float-point calculations precision)
and reconstruction is perfect for this filter bank.

Just to be accurate (taking into account filter normalization model an do so on),
let's keep in memory MATLAB 'db2' filter values to use it with IPP

>> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db2')

Lo_D =

Columns 1 through 3

-1.294095225509215e-001 2.241438680418574e-001 8.365163037374690e-001

Column 4

4.829629131446903e-001

Hi_D =

Columns 1 through 3

-4.829629131446903e-001 8.365163037374690e-001 -2.241438680418574e-001

Column 4

-1.294095225509215e-001

Lo_R =

Columns 1 through 3

4.829629131446903e-001 8.365163037374690e-001 2.241438680418574e-001

Column 4

-1.294095225509215e-001

Hi_R =

Columns 1 through 3

-1.294095225509215e-001 -2.241438680418574e-001 8.365163037374690e-001

Column 4

-4.829629131446903e-001

II. Let's try to repeat it with IPP:

start of .main.cpp
-----------
#include
#include "ipp.h"

// Filter bank for Db, order 2, filter

static const int fwdFltLenL = 4;
static const int fwdFltLenH = 4;

static const float fwdFltL[4] =
{
-1.294095225509215e-001f,
2.241438680418574e-001f,
8.365163037374690e-001f,
4.829629131446903e-001f
};

static const float fwdFltH[4] =
{
-4.829629131446903e-001f,
8.365163037374690e-001f,
-2.241438680418574e-001f,
-1.294095225509215e-001f
};

static const int invFltLenL = 4;
static const int invFltLenH = 4;

static const float invFltL[4] =
{
4.829629131446903e-001f,
8.365163037374690e-001f,
2.241438680418574e-001f,
-1.294095225509215e-001
};

static const float invFltH[4] =
{
-1.294095225509215e-001,
-2.241438680418574e-001,
8.365163037374690e-001,
-4.829629131446903e-001
};

// minimal values, that corresponds to MATLAB calculations
static const int fwdFltOffsL = -1;
static const int fwdFltOffsH = -1;

// minimal values, that corresponds to MATLAB calculations
static const int invFltOffsL = 0;
static const int invFltOffsH = 0;

void main(void)
{
Ipp32f src[] = {1, -10, 324, 48, -483, 4, 7, -5532, 34, 8889, -57, 54};

Ipp32f low[7];
Ipp32f high[7];

int i;

printf("original: ");
for(i = 0; i < 12; i++) printf("%.0f; ", src);
printf(" ");

{
// actually it can be single pointer for this offsets,
// but in general case it need to be different
Ipp32f srcExtBegL[2] = {src[10], src[11]}; // wrap [fwdFltLenL + fwdFltOffsL - 1] elements
Ipp32f srcExtBegH[2] = {src[10], src[11]}; // wrap [fwdFltLenH + fwdFltOffsH - 1] elements< /DIV>

Ipp32f srcExtEnd[] = { src[0], src[1] }; // to emulate MATLAB excessive 7's element, that repeat first

IppsWTFwdState_32f* state;

ippsWTFwdInitAlloc_32f (&state, fwdFltL, fwdFltLenL, fwdFltOffsL, fwdFltH, fwdFltLenH, fwdFltOffsH);

ippsWTFwdSetDlyLine_32f( state, srcExtBegL, srcExtBegH);

ippsWTFwd_32f (src, low, high, 6, state);
ippsWTFwd_32f (srcExtEnd, &low[6], &high[6], 1, state);

ippsWTFwdFree_32f (state);

// print decomposition result
printf("approx: ");
for(i = 0; i < 7; i++) printf("%.4f; ", low);
printf(" ");
printf("details: ");
for(i = 0; i < 7; i++) printf("%.4f; ", high);
printf(" ");

}

// note actually we simply do not need to keep MATLAB excessive low[6] & high[6] elements
// it can be applied as wrapped extension from first element
// let's change it to 'trash' to check it
low[6] = 1e20f;
high[6] = -1e10f;

{

// here is first data itself, that we need to place to delay line
// this is for non-shifted reconstruction for given conditions
// (filt. len. and offs.)
// see ippsWTInv_32f first call with &low[1] and &high[1]
// but in other conditions you may need to place wrapped data here also

Ipp32f extBegL[1] = {low[0]}; // [(invFltLenL + invFltOffsL - 1) / 2] elements
Ipp32f extBegH[1] = {high[0]}; // [(invFltLenH + invFltOffsH - 1) / 2] elements

// we decide do not use low[6] and high[6], so
// we need to extend borders honestly by wrapping
// and, please, note
// we do not use BIG buffer and copying, just input data itself
// and little extension in other memory place

Ipp32f extEndL[1] = {low[0]};
Ipp32f extEndH[1] = {high[0]};

Ipp32f dst[12];

IppsWTInvState_32f* state;

ippsWTInvInitAlloc_32f (&state, invFltL, invFltLenL, invFltOffsL, invFltH, invFltLenH, invFltOffsH);

ippsWTInvSetDlyLine_32f( state, extBegL, extBegH);

ippsWTInv_32f (&low[1], &high[1], 5, dst, state);
ippsWTInv_32f (extEndL, extEndH, 1, &dst[10], state);

ippsWTInvFree_32f (state);

// print reconstruction result
int i;
printf("reconstruction: ");
for(i = 0; i < 12; i++) printf("%.0f; ", dst);
printf(" ");
}
}

------
end of .main.cpp

Below is console output of this program:
-------
original:
1; -10; 324; 48; -483; 4; 7; -5532; 34; 8889; -57; 54;
approx:
19.1612; 58.5288; 87.8536; 487.5375; -5766.9277; 7432.4502; 19.1612;
details:
0.9387; 249.9611; -458.6568; 2739.2146; -3025.5576; -2070.5762; 0.9387;
reconstruction:
1; -10; 324; 48; -483; 4; 7; -5532; 34; 8889; -57; 54;
Press any key to continue
---------

The numbers of decomposition looks similar to MATLAB result,
and reconstruction is perfect.

Note,
we actually do not keep excessive expanding decomposition result
(but succesfully simulate their calculations to show compatibility),
at the same time we also do not use additional buffer to copy data (that take some time)
and wrap borders in small amount of non-a djacent memory.

And we successfully use MATLAB filter values,
but we need to care about parameters, that are 'hidden' in MATLAB
like extension model and filter anchors (offsets).
But you can see that we escape some of MATLAB artefacts,
like expanding decomposition.
Really, it's a bit strange to specify border-extension 'mode' and produce addtional data,
(and require it for idwt! )

***

Thank You & Regards,

Mikhail

[jimmyliu51]

Hi I have read the long program and tested it. I think it is a bit strange to call

the wavelet function twice to generate the result (for both forward and inverse). But it is a good idea indeed!

I modified the code a bit and attach below, by extending the data itself at the end

with the samples from the beginning, it "simulates" the "ppd" mode(periodical extention) of matlab's toolbox. My understand of the delay lineis more like padding the samples to the front of the data sequence. Once the padding of both ends of the data is done, I only need to call wavelet function in IPP once.

Just have to remember to select the correct number of results data samples in both forward WT and IWT. The variable resultSample is useful to control the number of forward WT, and to re-construct the result.

The original code that is not in use are still kept in the list for reference.

void main()
{
// TODO: Add your control notification handler code here
// Filter bank for Db, order 2, filter
static const int fwdFltLenL = 4;
static const int fwdFltLenH = 4;
static const float fwdFltL[4] =
{
-1.294095225509215e-001f,
2.241438680418574e-001f,
8.365163037374690e-001f,
4.829629131446903e-001f
};
static const float fwdFltH[4] =
{
-4.829629131446903e-001f,
8.365163037374690e-001f,
-2.241438680418574e-001f,
-1.294095225509215e-001f
};
static const int invFltLenL = 4;
static const int invFltLenH = 4;
static const float invFltL[4] =
{
4.829629131446903e-001f,
8.365163037374690e-001f,
2.241438680418574e-001f,
-1.294095225509215e-001f
};

static const float invFltH[4] =
{
-1.294095225509215e-001f,
-2.241438680418574e-001f,
8.365163037374690e-001f,
-4.829629131446903e-001f
};
// minimal values, that corresponds to MATLAB calculations
static const int fwdFltOffsL = -1;
static const int fwdFltOffsH = -1;
// minimal values, that corresponds to MATLAB calculations
static const int invFltOffsL = 0;
static const int invFltOffsH = 0;

int filterlength=4;
int sampleNumber=12;
int resultSamples=(filterlength+sampleNumber-1)/2;

Ipp32f src[] = {1, -10, 324, 48, -483, 4, 7, -5532, 34, 8889, -57, 54};
Ipp32f src1[16];
Ipp32f low[7];
Ipp32f high[7];
int i;
FILE *out;
out=fopen("ippoutput0.txt", "w");
fprintf(out, "original: ");
for(i = 0; i < 12; i++) fprintf(out, "%.0f; ", src);
fprintf(out, " ");
{
// actually it can be single pointer for this offsets,
// but in general case it need to be different
Ipp32f srcExtBegL[3] = {src[9], src[10], src[11]}; // wrap [fwdFltLenL + fwdFltOffsL - 1] elements
Ipp32f srcExtBegH[3] = {src[9], src[10], src[11]}; // wrap [fwdFltLenH + fwdFltOffsH - 1] elements
Ipp32f srcExtEnd[] = { src[0], src[1] }; // to emulate MATLAB excessive 7's element, that repeat first
IppsWTFwdState_32f* state;

ippsWTFwdInitAlloc_32f (&state, fwdFltL, fwdFltLenL, fwdFltOffsL, fwdFltH, fwdFltLenH, fwdFltOffsH);
ippsWTFwdSetDlyLine_32f( state, srcExtBegL, srcExtBegH);
#if 0
ippsWTFwd_32f (src, low, high, 6, state);
ippsWTFwd_32f (srcExtEnd, &low[6], &high[6], 1, state);
#endif
for (i=0; i<12; i++)
src1=src;
for (i=0; i
src1[12+i]=src[i ];
// (12+4-1)/2 =7
ippsWTFwd_32f (src1, low, high, resultSamples, state);
//ippsWTFwd_32f (srcExtEnd, &low[6], &high[6], 1, state);
ippsWTFwdFree_32f (state);
// print decomposition result
fprintf(out, "approx: ");
for(i = 0; i < resultSamples; i++) fprintf(out, "%.4f; ", low);
fprintf(out, " ");
printf("details: ");
for(i = 0; i < resultSamples; i++) fprintf(out, "%.4f; ", high);
fprintf(out, " ");
}
// note actually we simply do not need to keep MATLAB excessive low[6] & high[6] elements
// it can be applied as wrapped extension from first element
// let's change it to 'trash' to check it
#if 0
low[6] = 1e20f;
high[6] = -1e10f;
#endif
{
// here is first data itself, that we need to place to delay line
// this is for non-shifted reconstruction for given conditions
// (filt. len. and offs.)
// see ippsWTInv_32f first call with &low[1] and &high[1]
// but in other conditions you may need to place wrapped data here also
Ipp32f extBegL[1] = {low[0]}; // [(invFltLenL + invFltOffsL - 1) / 2] elements
Ipp32f extBegH[1] = {high[0]}; // [(invFltLenH + invFltOffsH - 1) / 2] elements

// we decide do not use low[6] and high[6], so
// we need to extend borders honestly by wrapping
// and, please, note
// we do not use BIG buffer and copying, just input data itself
// and little extension in other memory place
#if 0
Ipp32f extEndL[1] = {low[0]};
Ipp32f extEndH[1] = {high[0]};
#endif
Ipp32f dst2[14];
IppsWTInvState_32f* state;
ippsWTInvInitAlloc_32f (&state, invFltL, invFltLenL, invFltOffsL, invFltH, invFltLenH, invFltOffsH);
ippsWTInvSetDlyLine_32f( state, extBegL, extBegH);
#if 0
ippsWTInv_32f (&low[1], &high[1], 5, dst, state);
ippsWTInv_32f (extEndL, extEndH, 1, &dst[10], state);
#endif
ippsWTInv_32f (low, high, resultSamples, dst2, state);
ippsWTInvFree_32f (state);
// print reconstruction result
int i;
int extraHead=2*resultSamples-sampleNumber;
fprintf(out, "reconstruction: ");
for(i = 0; i < 12; i++) fprintf(out, "%.0f; ", dst2[extraHead+i]);
fprintf(out, " ");
}

}

[Mikhail_Kulikov__Int]

Hi!
Thank you for development of original code,
and, that's true, the code may be simpler, and if you do not need extremely low performance time this is the best and the right way to implement different extensions models (like it was done in MATLAB).
But if you calculate this way it requires the use of excessive copying to the additional buffer as for example:

Code:

  for (i=0; i<12; i++)
   src1=src;
  for (i=0; i
   src1[12+i]=src;

And it takes some performance time for large vectors. You may escape it, but in such case it requires to care about some details (like two calls of ippsWTFwd or ippsWTInv), that were illustrated in the first version.
In general this pair of examples shows how WT functionality may be used in both cases, when you need minimize performance time or you need minimize implementation resources. This alternative is really actual in the development of software.

Thank You and Regards!

Mikhail

[jimmyliu51]

Thank you for your reply!
I don't really use the way in the following code to implement my application. I just show the need to pad the extra data samples. But since src[12] is not long enough, I put another src1[] in use. All I need to do is to have a slightly large buffer (size = signal length +wavelet filter length-1), and always remember to pad the additional data samples at the end. My application is using data from a round disk, so it is naturally I need to use periodical padding.

My point is that it is possible to pad the data in such way that we can have the same result as "ppd" mode in wavelet toobox of Matlab.

In fact, I still do not really understand why the offset for forward WT is -1 and inverse WT is 0. I did change these values and I end up with very different value. if you can, could you give me some idea? Thank you!

Code:

  for (i=0; i<12; i++)
   src1=src;
  for (i=0; i
   src1[12+i]=src;

[Mikhail_Kulikov__Int]

Regarding offset values and minimal offset value (-1).

Because so called WT low-pass and high-pass filter responses may have maximum (or optimal in any other criteria) value in different position, there is a task to align them each with other. The delay line kind of interface is mostly intuitive for streaming data processing (as it was considered above) and in the streaming applications it can be done by additional delay to the least inertial filter to align with the most inertial and its a half of the reason to introduce offsets (additional delays) for filters.
If we look into perfect reconstruction nature (in the other words in the particular criteria), we may find there some conditions for filters alignment to successfully reconstruct original data (even if we talk about delayed perfect reconstruction), they may be different for different types of filter banks and its another half of reason to introduce filter alignment flexibility in interface.

The minimal value of aligning delay or offset (-1) allows to obtain minimal possible latency for WT analysis and subsequent reconstruction. Its just adjustment to specify minimal portion of pre-history that is required for multi-rate filtering.
Lets consider some example. The minimal length of filter-bank, which allows successful reconstruction, is 2 (Haar). If it would be usual FIR filter (without resampling) you are restricted by latency of len -1 that is 1 sample. But because of resampling with factor 2 you can throw out the first RESULTING sample and take only second (fourth and so on), so you simply do not need delay line at all, if you work with even size of buffering. It just shows that its possible for Haar to decompose data by pairs (l=src[0]+ src[1], h=src[1]-src[0]), without any actual data in mind (i.e. in state structure delay line) between blocks of length 2 x N. And the interface provide such a possibility by offset = -1 and you do not need pre-history for Haar-analyzed data block.
Its not a way-back machine:) and if buffer size restricted by single sample (or any other odd number) the external one-sample delay line will be required to provide even number of input data samples. But for buffer of even size its possible to implement 1-level Haar analysis without additional latency. Any analysis eats some samples.

WBR,
Mikhail

[Mikhail_Kulikov__Int]

P.S. some addition info regarding offset values for Db2 and other orthogonal filter banks with "wraped" borders extension.

If you need minimal latency start from -1 value for Fwd transforms for both filters, like in code examples above. It should be enough. For Inv transform you need to provide offset to have result in the right place. The zero value is the most suitable for example above.

But actually offsets are just additional delays to the input dataof Fwd and for output data of Inv filtering (and L & H values may be specified independently).

Regards,

Mikhail