Lifting

This example shows how to use lifting on a 1-D signal.

Create a 1-D signal that is piecewise constant over 2 samples. Add $$N(0,0.1^2)$$ noise to the signal.

x = [1 1 2 2 -3.5 -3.5 4.3 4.3 6 6 -4.5 -4.5 2.2 2.2 -1.5 -1.5];
x = repmat(x,1,64);
rng default
x = x+ 0.1*randn(size(x));

Plot the signal and zoom in on the first 100 samples to visualize the correlation in neighboring samples.

plot(x)
xlim([0 100])

Use the lazy wavelet to obtain the even and odd polyphase components of the signal.

LS = liftwave('lazy');
[A,D] = lwt(x,LS);

If you plot the detail (wavelet) coefficients in D, you see that this transform has not decorrelated the signal. The wavelet coefficients look very much like the signal.

Add a dual lifting step that subtracts the even-indexed coefficient from the odd-coefficient one sample later, $$x(2n+1)-x(2n)$$.

els = {'d',-1,0};
LSnew = addlift(LS,els);

Because the signal is piecewise constant over consecutive samples with additive noise, the new dual lifting step should result in wavelet coefficients small in absolute value. In this case, the wavelet transform does decorrelate the data. Verify this by finding the approximation and detail coefficients with the new dual lifting step.

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

If you plot the detail (wavelet) coefficients, you see that the wavelet coefficients no longer resemble the original signal.

The approximation coefficients, A, of the previous transform constitute the even polyphase component of the signal. Therefore, the coefficients are affected by aliasing. Use a primal lifting step to update the approximation coefficients and reduce aliasing. The primal step replaces the approximation coefficients by $$x(2n)+1/2(x(2n+1)-x(2n))$$, which is equal to the average of $$x(2n)$$ and $$x(2n+1)$$. The averaging is a lowpass filtering, which helps to reduce aliasing.

els = {'p',1/2, 0};
LSnew = addlift(LSnew,els);

Use the updated lifting scheme to obtain the wavelet transform of the input signal.

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

Add the appropriate scaling to ensure perfect reconstruction. Obtain the approximation and wavelet coefficients using lifting and reconstruct the signal using ilwt. Verify perfect reconstruction.

LSnew(end,:) = {sqrt(2),sqrt(2)/2,[]};
[A,D] = lwt(x,LSnew);
x1 = ilwt(A,D,LSnew);
max(abs(x1-x))

ans = 1.7764e-15

The preceding example designed a wavelet, which effectively removed a zeroth order polynomial (constant). If the behavior of the signal is better represented by a higher-order polynomial, you can design a dual wavelet with the appropriate number of vanishing moments to decorrelate the signal.

Use the lifting scheme to design a wavelet with 2 vanishing moments. A dual wavelet with 2 vanishing moments decorrelates a signal with local behavior approximated by a first-order polynomial. Create a signal characterized by first-order polynomial behavior with additive $$N(0,0.25^2)$$ noise.

y = [1 0 0 4 0 0 -1 0 0 2 0 0 7 0 0 -4 0 0 1 0 0 -3];
x1 = 1:(21/1024):22-(21/1024);
y1 = interp1(1:22,y,x1,'linear');
rng default
y1 = y1+0.25*randn(size(y1));
plot(x1,y1)
xlim([1 22])

In this case, the wavelet coefficients should remove a first-order polynomial. If the signal value at an odd index, $$x(2n+1)$$, is well approximated by a first-order polynomial fitted to the surrounding sample values, then $$1/2(x(2n)+x(2n+2))$$ should provide a good fit for $$x(2n+1)$$. In other words, $$x(2n+1)$$ should be the midpoint between $$x(2n)$$ and $$x(2n+2)$$.

It follows that $$x(2n+1)-1/2(x(2n)+x(2n+2))$$ should decorrelate the signal.

Start with the lazy wavelet transform and add a dual lifting step which models the preceding equation.

LS = liftwave('lazy');
els = {'d',[-1/2 -1/2],1};
LSnew = addlift(LS,els);

Use the lifting scheme to obtain the approximation and detail coefficients and plot the result.

[A,D] = lwt(y1,LSnew);
subplot(2,1,1)
plot(A)
xlim([1 512])
title('Approximation Coefficients')
subplot(2,1,2)
plot(D)
xlim([1 512])
title('Detail Coefficients')

You see that the wavelet coefficients appear to only contain noise, while the approximation coefficients represent a denoised version of the original signal. Because the preceding transform uses only the even polyphase component for the approximation coefficients, you can reduce aliasing by adding a primal lifting step. Finally, add the normalization constants to produce a perfect reconstruction filter bank.

Obtain the discrete wavelet transform with the new lifting scheme and plot the results.

els = {'p',[1/4 1/4],0};
LSnew = addlift(LSnew,els);
LSnew(end,:) = {sqrt(2),sqrt(2)/2,[]};
[A,D] = lwt(y1,LSnew);
subplot(2,1,1)
plot(A)
xlim([1 512])
title('Approximation Coefficients')
subplot(2,1,2)
plot(D)
xlim([1 512])
title('Detail Coefficients')

Demonstrate that you have designed a perfect reconstruction filter bank.

y2 = ilwt(A,D,LSnew);
max(abs(y2-y1))

ans = 8.8818e-16