Main Content

Lifting a Filter Bank

This example shows how to use lifting to progressively change the properties of a perfect reconstruction filter bank. The following figure shows the three canonical steps in lifting: split, predict, and update.

The first step in lifting is simply to split the signal into its even- and odd-indexed samples. These are called polyphase components and that step in the lifting process is often referred to as the "lazy" lifting step because you really are not doing that much work. You can do this in MATLAB by creating a "lazy" lifting scheme.

LS = liftwave('lazy');

Apply the lifting scheme to some data.

x = randn(8,1);
[ALazy,DLazy] = lwt(x,LS);

MATLAB™ indexes from 1 so ALazy contains the odd-indexed samples of x and DLazy contains the even-indexed samples. Most explanations of lifting assume that the signal starts with sample 0, so ALazy would be the even-indexed samples and DLazy the odd-indexed samples. This example follows that latter convention. The "lazy" wavelet transform treats one half of the signal as wavelet coefficients, DLazy, and the other half as scaling coefficients, ALazy. This is perfectly consistent within the context of lifting, but a simple split of the data does really sparsify or capture any relevant detail.

The next step in the lifting scheme is to predict the odd samples based on the even samples. The theoretical basis for this is that most natural signals and images exhibit correlation among neighboring samples. Accordingly, you can "predict" the odd-indexed samples using the even-indexed samples. The difference between your prediction and the actual value is the "detail" in the data missed by the predictor. That missing detail comprises the wavelet coefficients.

The prediction step is also referred to as a "dual lifting step". In equation form, you can write the prediction step as dj(n)=dj-1(n)-P(aj-1) where dj-1(n) are the wavelet coefficients at the finer scale and aj-1 is some number of finer-scale scaling coefficients. P() is the prediction operator.

Add a simple (Haar) dual lifting step that subtracts the even (approximation) coefficient from the odd (detail) coefficient. In this case the prediction operator is simply (-1)aj-1(n). In other words, it predicts the odd samples based on the immediately preceding even sample.

ElemLiftStep = {'d',-1,0};

The above code says "create an elementary dual (predict) lifting step using a polynomial in z with the highest power z0. The coefficient is -1. Update the lazy lifting scheme.

LSN = addlift(LS,ElemLiftStep,'end');

Apply the new lifting scheme to the signal.

[A,D] = lwt(x,LSN);

Note that the elements of A are identical to those in ALazy. This is expected because you did not modify the approximation coefficients. If you look at the elements of D, you see that they are equal to

Dnew = DLazy-ALazy;

Compare Dnew to D. Imagine an example where the signal was piecewise constant over every two samples.

v = [1 -1 1 -1 1 -1];
u = repelem(v,2);

Apply the new lifting scheme to u.

[Au,Du] = lwt(u,LSN);

You see that all the Du are zero. This signal has been compressed because all the information is now contained in 6 samples instead of 12 samples. You can easily reconstruct the original signal

urecon = ilwt(Au,Du,LSN);

In your prediction (dual lifting) step, you predicted that the adjacent odd sample in your signal had the same value as the immediately preceding even sample. Obviously, this is true only for trivial signals. The wavelet coefficients capture the difference between the prediction and the actual values (at the odd samples). Finally, use the update step to update the even samples based on differences obtained in the prediction step. In this case, update using the following aj(n)=aj-1(n)+dj-1(n)/2. This replaces each even-indexed coefficient by the arithmetic average of the even and odd coefficients. An update step is also referred to as a primal lifting step.

elsprimal = {'p',1/2,0};
LSupdated = addlift(LSN,elsprimal,'end');

Obtain the wavelet transform of the signal with the updated lifting scheme.

[A,D] = lwt(x,LSupdated);

If you compare A to the original signal, x, you see that the signal mean is captured in the approximation coefficients.

ans = -0.0131
ans = -0.0131

In fact, the elements of A are easily obtainable from x by the following.

n = 1;
for ii = 1:2:numel(x)
    meanz(n) = mean([x(ii) x(ii+1)]);
    n = n+1;

Compare meanz and A. As always, you can invert the lifting scheme to obtain a perfect reconstruction of the data.

xrec = ilwt(A,D,LSupdated);
ans = 2.2204e-16

It is common to add a normalization step at the end so that the energy in the signal (2 norm) is preserved as the sum of the energies in the scaling and wavelet coefficients. Without this normalization step, the energy is not preserved.

ans = 11.6150
ans = 16.8091

Add the necessary normalization step.

LSscaled = LSupdated;
LSscaled(end,1:2) = {sqrt(2), sqrt(2)/2};
[A,D] = lwt(x,LSscaled);
ans = 11.6150

Now the 2 norm of the signal is equal to the sum of the energies in the scaling and wavelet coefficients. The lifting scheme you developed in this example is the Haar lifting scheme.

The Wavelet Toolbox™ supports many commonly used lifting schemes through liftwave with pre-defined dual, primal, and normalization steps. For example, you can obtain the Haar lifting scheme with the following.

lshaar = liftwave('haar');

If you compare lshaar to LSUpdated, you see that our step-by-step lifting scheme matches the Haar lifting scheme. To see that not all lifting schemes consist of single dual and primal lifting steps, examine the lifting scheme that corresponds to the 'bior3.1' wavelet.

lsbior3_1 = liftwave('bior3.1')
lsbior3_1=4×3 cell array
    {'d'     }    {[  0.3333]}    {[       1]}
    {'p'     }    {1x2 double}    {[       0]}
    {'d'     }    {[ -0.4444]}    {[       0]}
    {[0.4714]}    {[  2.1213]}    {0x0 double}

You can also use liftfilt if you want to start with a set of biorthogonal (or orthogonal) scaling and wavelet filters and "lift" them to another set. For example, start with the Haar scaling and lifting filters.

[LoD,HiD,LoR,HiR] = wfilters('haar');

Lift the Haar filters with two primal lifting steps.

twoels(1) = struct('type','p','value',...
laurpoly([0.125 -0.125],0));
twoels(2) = struct('type','p','value',...
laurpoly([0.125 -0.125],1));
[LoDN,HiDN,LoRN,HiRN] = liftfilt(LoD,HiD,LoR,HiR,twoels);

Plot the resulting scaling and wavelet functions.

[phia,psia,phis,psis,xval] = bswfun(LoDN,HiDN,LoRN,HiRN);
title('Analysis Scaling Function');
axis tight;
grid on;
axis tight;
grid on;
title('Synthesis Scaling Function');
axis tight;
grid on;
title('Analysis Wavelet');
axis tight;
grid on;
title('Synthesis Wavelet');

If you plot the analysis and synthesis scaling functions and wavelets for the 'bior1.3' wavelet, you see that lifting the Haar wavelet as in the previous example has essentially provided the 'bior1.3' wavelet to within a change of sign on the synthesis wavelet.

[LoD,HiD,LoR,HiR] = wfilters('bior1.3');
[phia,psia,phis,psis,xval] = bswfun(LoD,HiD,LoR,HiR);