www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/DualTreeWaveletTransformsExample.m
%% Dual-Tree Wavelet Transforms
% This example shows how the dual-tree discrete wavelet transform (DWT)
% provides advantages over the critically-sampled DWT for signal and image
% processing. The dual-tree DWT is implemented as two separate two-channel
% filter banks. To gain the advantages described in this example, you
% cannot arbitrarily choose the scaling and wavelet filters used in the two
% trees. The lowpass (scaling) and highpass (wavelet) filters of one tree,
% $\{h_0,h_1\}$ must generate a scaling function and wavelet that are
% approximate Hilbert transforms of the scaling function and wavelet
% generated by the lowpass and highpass filters of the other tree,
% $\{g_0,g_1\}$. Therefore, the complex-valued scaling function and wavelet
% formed from the two trees are approximately analytic.
%
% As a result, the dual-tree (complex) DWT exhibits less shift variance and
% more directional selectivity than the critically sampled DWT with only
% a $2^d$ redundancy factor for d-dimensional data. The redundancy in the
% dual-tree DWT is significantly less than the redundancy in the
% undecimated (stationary) DWT.
%
% This example illustrates the approximate shift invariance of the dual-
% tree DWT, the selective orientation of the dual-tree analyzing wavelets
% in 2-D, and the use of the dual-tree DWT in image denoising.
% Copyright 2015 The MathWorks, Inc.
%% Near Shift Invariance of the Dual-Tree DWT
% The DWT suffers from shift variance, meaning that small shifts in the
% input signal or image can cause significant changes in the distribution
% of signal/image energy across scales in the DWT coefficients. The complex
% dual-tree DWT is approximately shift invariant.
%
% To demonstrate this on a test signal, construct two shifted discrete-time
% impulses 128 samples in length. One signal has the unit impulse at sample
% 60, while the other signal has the unit impulse at sample 64. Both
% signals clearly have unit energy ($l^2$ norm).
kronDelta1 = zeros(128,1);
kronDelta1(60) = 1;
kronDelta2 = zeros(128,1);
kronDelta2(64) = 1;
%%
% Obtain the DWT and dual-tree DWT of the two signals down to level 3
% with wavelet and scaling filters of length 14. Extract the level-3 detail
% coefficients for comparison.
J = 3;
dwt1 = dddtree('dwt',kronDelta1,J,'sym7');
dwt2 = dddtree('dwt',kronDelta2,J,'sym7');
dwt1Cfs = dwt1.cfs{J};
dwt2Cfs = dwt2.cfs{J};
dt1 = dddtree('cplxdt',kronDelta1,J,'dtf3');
dt2 = dddtree('cplxdt',kronDelta2,J,'dtf3');
dt1Cfs = dt1.cfs{J}(:,:,1)+1i*dt1.cfs{J}(:,:,2);
dt2Cfs = dt2.cfs{J}(:,:,1)+1i*dt2.cfs{J}(:,:,2);
%%
% Plot the absolute values of the DWT and dual-tree DWT coefficients for
% the two signals at level 3 and compute the energy (squared $l^2$ norms)
% of the coefficients.
figure;
subplot(1,2,1)
stem(abs(dwt1Cfs),'markerfacecolor',[0 0 1]);
title(['DWT Squared 2-norm = ' num2str(norm(dwt1Cfs,2)^2)]);
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title(['DWT Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2)]);
figure;
subplot(1,2,1)
stem(abs(dt1Cfs),'markerfacecolor',[0 0 1]);
title(['Dual-tree DWT Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2)]);
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title(['Dual-tree DWT Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2)]);
%%
% Note the four sample shift in the signal has caused an almost 6.5% change
% in the energy of the level-3 DWT wavelet coefficients. However, the dual-
% tree wavelet coefficients show only a 0.3% change in energy.
%
% To demonstrate the utility of approximate shift invariance in real data,
% we analyze an electrocardiogram (ECG) signal. The sampling interval for
% the ECG signal is 1/180 seconds. The data are taken from Percival &
% Walden (2000), p.125 (data originally provided by William Constantine
% and Per Reinhall, University of Washington). For convenience, we take the
% data to start at t=0.
load wecg
dt = 1/180;
t = 0:dt:(length(wecg)*dt)-dt;
figure;
plot(t,wecg)
xlabel('Seconds'); ylabel('Millivolts');
%%
% The large positive peaks approximately 0.7 seconds apart are the R waves
% of the cardiac rhythm. First, decompose the signal using the critically
% sampled DWT and plot the original signal along with the level-2 and
% level-3 wavelet coefficients. The level-2 and level-3 coefficients were
% chosen because the R waves are isolated most prominently in those scales.
J = 6;
dtDWT1 = dddtree('dwt',wecg,J,'farras');
details = zeros(2048,3);
details(2:4:end,2) = dtDWT1.cfs{2};
details(4:8:end,3) = dtDWT1.cfs{3};
subplot(311)
stem(t,details(:,2),'Marker','none','ShowBaseline','off')
title('Level 2'); ylabel('mV');
subplot(312)
stem(t,details(:,3),'Marker','none','ShowBaseline','off')
title('Level 3'); ylabel('mV');
subplot(313)
plot(t,wecg); title('Original Signal');
xlabel('Seconds'); ylabel('mV');
%%
% Repeat the above analysis for the dual-tree transform. In this case, just
% plot the real part of the dual-tree coefficients at levels 2 and 3.
dtcplx1 = dddtree('cplxdt',wecg,J,'dtf3');
details = zeros(2048,3);
details(2:4:end,2) = dtcplx1.cfs{2}(:,1,1)+1i*dtcplx1.cfs{2}(:,1,2);
details(4:8:end,3) = dtcplx1.cfs{3}(:,1,1)+1i*dtcplx1.cfs{3}(:,1,2);
subplot(311)
stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off')
title('Level 2'); ylabel('mV');
subplot(312)
stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off')
title('Level 3'); ylabel('mV');
subplot(313)
plot(t,wecg); title('Original Signal');
xlabel('Seconds'); ylabel('mV');
%%
% Both the critically sampled and dual-tree wavelet transforms localize an
% important feature of the ECG waveform to similar scales. An important
% application of wavelets in 1-D signals is to obtain an analysis of
% variance by scale. It stands to reason that this analysis of variance
% should not be sensitive to circular shifts in the input signal.
% Unfortunately, this is not the case with the critically sampled DWT.
% To demonstrate this, we circularly shift the ECG signal by 4 samples,
% analyze the unshifted and shifted signals with the critically sampled
% DWT, and calculate the distribution of energy across scales.
wecgShift = circshift(wecg,4);
dtDWT2 = dddtree('dwt',wecgShift,J,'farras');
sigenrgy = norm(wecg,2)^2;
enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT1.cfs,'uni',0));
enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT2.cfs,'uni',0));
levels = {'D1';'D2';'D3';'D4';'D5';'D6';'A6'};
enr1 = enr1(:);
enr2 = enr2(:);
table(levels,enr1,enr2,'VariableNames',{'Level','enr1','enr2'})
%%
% Note that the wavelet coefficients at levels 3 and 4 show approximately
% 3% changes in energy between the original and shifted signal. Next, we
% repeat this analysis using the complex dual-tree wavelet transform.
%
% The dual-tree transform produces a consistent analysis of variance by
% scale for the original signal and its circularly shifted version.
dtcplx2 = dddtree('cplxdt',wecgShift,J,'dtf3');
cfs1 = cellfun(@squeeze,dtcplx1.cfs,'uni',0);
cfs2 = cellfun(@squeeze,dtcplx2.cfs,'uni',0);
cfs1 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs1,'uni',0);
cfs2 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs2,'uni',0);
dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0));
dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'uni',0));
dtenr1 = dtenr1(:);
dtenr2 = dtenr2(:);
table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})
%% Directional Selectivity in Image Processing
% The standard implementation of the DWT uses separable filtering of the
% columns and rows of the image. The LH, HL, and HH wavelets for
% Daubechies' least-asymmetric phase wavelet with 4 vanishing moments
% (sym4) are shown in the following figure.
figure;
J = 5;
L = 3*2^(J+1);
N = L/2^J;
Y = zeros(L,3*L);
dt = dddtree2('dwt',Y,J,'sym4');
dt.cfs{J}(N/3,N/2,1) = 1;
dt.cfs{J}(N/2,N/2+N,2) = 1;
dt.cfs{J}(N/2,N/2+2*N,3) = 1;
dwtImage = idddtree2(dt);
imagesc(dwtImage); axis xy; axis off;
title({'Critically Sampled DWT';...
'2-D separable wavelets (sym4)- LH, HL, HH'});
%%
% Note that the LH and HL wavelets have clear horizontal and vertical
% orientations respectively. However, the HH wavelet on the far right mixes
% both the +45 and -45 degree directions, producing a checkerboard
% artifact. This mixing of orientations is due to the use of real-valued
% separable filters. The HH real-valued separable filter has passbands in
% all four high frequency corners of the 2-D frequency plane.
%
% The dual-tree DWT achieves directional selectivity by using wavelets that
% are approximately analytic, meaning that they have support on only one
% half of the frequency axis. In the dual-tree DWT, there are six subbands
% for both the real and imaginary parts. The six real parts are formed by
% adding the outputs of column filtering followed by row filtering of the
% input image in the two trees. The six imaginary parts are formed by
% subtracting the outputs of column filtering followed by row filtering.
%
% The filters applied to the columns and rows may be from the same filter
% pair, $\{h_0,h_1\}$ or $\{g_0,g_1\}$, or from different filter pairs,
% $\{h_0,g_1\}$, $\{g_0,h_1\}$. The following code shows the orientation
% of the 12 wavelets corresponding to the real and imaginary parts of the
% complex oriented dual-tree DWT.
J = 4;
L = 3*2^(J+1);
N = L/2^J;
Y = zeros(2*L,6*L);
wt = dddtree2('cplxdt',Y,J,'dtf3');
wt.cfs{J}(N/2,N/2+0*N,2,2,1) = 1;
wt.cfs{J}(N/2,N/2+1*N,3,1,1) = 1;
wt.cfs{J}(N/2,N/2+2*N,1,2,1) = 1;
wt.cfs{J}(N/2,N/2+3*N,1,1,1) = 1;
wt.cfs{J}(N/2,N/2+4*N,3,2,1) = 1;
wt.cfs{J}(N/2,N/2+5*N,2,1,1) = 1;
wt.cfs{J}(N/2+N,N/2+0*N,2,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+1*N,3,1,2) = 1;
wt.cfs{J}(N/2+N,N/2+2*N,1,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+3*N,1,1,2) = 1;
wt.cfs{J}(N/2+N,N/2+4*N,3,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+5*N,2,1,2) = 1;
waveIm = idddtree2(wt);
imagesc(waveIm); axis off;
title('Complex Oriented Dual-Tree 2-D Wavelets');
%%
% The top row of the preceding figure shows the six directional wavelets
% of the real oriented dual-tree wavelet transform. The second row shows
% the imaginary parts. Together the real and imaginary parts form the
% complex oriented dual-tree wavelet transform. The real and imaginary
% parts are oriented in the same direction. You can use dddtree2 with the
% 'realdt' option to obtain the real oriented dual-tree DWT, which uses
% only the real parts. Using the real oriented dual-tree wavelet transform,
% you can achieve directional selectivity, but you do not gain the full
% benefit of using analytic wavelets such as approximate shift invariance.
%% Image Denoising
% The dual-tree DWT isolates distinct orientations is separate subbands and
% outperforms the standard separable DWT in applications like image
% denoising.
%
% This example uses a helper function, helperCompare2DDenoising, to load
% an image and adds zero-mean white Gaussian noise with sigma equal to 25.
% For a user-supplied range of thresholds, the function compares denoising
% using soft thresholding for the critically sampled DWT, the real oriented
% dual-tree DWT, and the complex oriented dual-tree DWT. For each threshold
% value, the root-mean-square (RMS) error and peak signal-to-noise ratio
% (PSNR) are displayed.
numex = 3;
helperCompare2DDenoising(numex,0:2:100,'PlotMetrics');
%%
% Both the real oriented and the complex oriented dual-tree DWTs outperform
% the standard DWT in RMS error and PSNR.
%%
% Next, obtain the denoised images for a threshold value of 25, which is
% equal to the standard deviation of the additive noise.
numex = 3;
helperCompare2DDenoising(numex,25,'PlotImage');
%%
% With a threshold value equal to the standard deviation of the additive
% noise, the complex oriented dual-tree transform provides a PSNR almost
% 4 dB higher than the standard 2-D DWT.
%% Edge Representation in Two Dimensions
% The approximate analyticity and selective orientation of the complex
% dual-tree wavelets provide superior performance over the standard 2-D DWT
% in the representation of edges in images. To illustrate this, we analyze
% test images with edges consisting of line and curve singularities in
% multiple directions using the critically sampled 2-D DWT and the 2-D
% complex oriented dual-tree transform. First, analyze an image of an
% octagon, which consists of line singularities.
load woctagon;
figure;
imagesc(woctagon); colormap gray;
title('Original Image'); axis equal; axis off;
%%
% Decompose the image down to level 4 and reconstruct an image
% approximation based on the level-4 detail coefficients.
dtcplx = dddtree2('cplxdt',woctagon,4,'dtf3');
dtDWT = dddtree2('dwt',woctagon,4,'farras');
dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1}));
dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2}));
dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3}));
dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5}));
dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1}));
dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2}));
dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3}));
dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5}));
dtImage = idddtree2(dtcplx);
dwtImage = idddtree2(dtDWT);
subplot(121)
imagesc(dtImage); axis equal; axis off; colormap gray;
title('Complex Oriented Dual-Tree');
subplot(122)
imagesc(dwtImage); axis equal; axis off; colormap gray;
title('DWT')
%%
% Next, analyze an octagon with hyperbolic edges. The edges in the
% hyperbolic octagon are curve singularities.
load woctagonHyperbolic;
figure;
imagesc(woctagonHyperbolic); colormap gray;
title('Octagon with Hyperbolic Edges'); axis equal; axis off;
%%
% Decompose the image down to level 4 and reconstruct an image
% approximation based on the level-4 detail coefficients for both the
% standard 2-D DWT and the complex oriented dual-tree DWT.
dtcplx = dddtree2('cplxdt',woctagonHyperbolic,4,'dtf3');
dtDWT = dddtree2('dwt',woctagonHyperbolic,4,'farras');
dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1}));
dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2}));
dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3}));
dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5}));
dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1}));
dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2}));
dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3}));
dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5}));
dtImage = idddtree2(dtcplx);
dwtImage = idddtree2(dtDWT);
subplot(121)
imagesc(dtImage); axis equal; axis off; colormap gray;
title('Complex Oriented Dual-Tree');
subplot(122)
imagesc(dwtImage); axis equal; axis off; colormap gray;
title('DWT')
%%
% Note that the ringing artifacts evident in the 2-D critically sampled DWT
% are absent in the 2-D complex dual-tree transforms of both images. The
% complex oriented dual-tree DWT more faithfully reproduces line and curve
% singularities.
%% Summary
% We have shown that the dual-tree DWT possesses the desirable properties
% of near shift invariance and directional selectivity not achievable with
% the critically sampled DWT. We have demonstrated how these properties can
% result in improved performance in signal analysis, the representation of
% singularities in images, and image denoising. In addition to the real-
% oriented and complex-oriented dual-tree DWT, dddtree and dddtree2 also
% support the double-density wavelet transform and dual-tree double-density
% wavelet transforms, which are additional examples of overcomplete wavelet
% filter banks (frames) with advantages over the standard DWT.
%% Further Reading
% Kingsbury, N.G. "Complex Wavelets for Shift Invariant Analysis and
% Filtering of Signals". Journal of Applied and Computational Harmonic
% Analysis. Vol 10, Number 3, May 2001, pp. 234-253.
%
% Percival, D.B. and A.T. Walden. "Wavelet Methods for Time Series
% Analysis", Cambridge University Press, 2000.
%
% Selesnick, I., Baraniuk, R.G., and N.G. Kingsbury. "The Dual-Tree
% Complex Wavelet Transform." IEEE Signal Processing Magazine. Vol. 22,
% Number 6, November, 2005, pp. 123-151.
%
% Selesnick, I. "The Double Density DWT". Wavelets in Signal and Image
% Analysis: From Theory to Practice (A.A Petrosian, F.G. Meyer, eds.),
% Norwell, MA: Kluwer Academic Publishers, 2001, pp. 39-66.
%
% Selesnick, I. "The Double-Density Dual-Tree Wavelet Transform". IEEE
% Transactions on Signal Processing. Vol. 52, Number 5, May 2004, pp.
% 1304-1314.