www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/LiftingExample.m

    %% 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);
set(gca,'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)) 

%%
% The preceding example designed a wavelet, which effectively removed a
% zero-th 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); set(gca,'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(211)
plot(A); set(gca,'xlim',[1 512]);
title('Approximation Coefficients');
subplot(212)
plot(D); set(gca,'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(211)
plot(A); set(gca,'xlim',[1 512]);
title('Approximation Coefficients');
subplot(212)
plot(D); set(gca,'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))