matlab-双树复小波变换DTCWT(转自 MathWorks)

此示例展示了双树复小波变换 (DTCWT) 如何在信号、图像和体积处理方面提供优于临界采样 DWT 的优势。DTCWT 作为两个独立的双通道滤波器组来实现。为了获得本示例中描述的优点，您不能任意选择两个树中使用的缩放和小波滤波器。一棵树的低通（缩放）和高通（小波）滤波器，{H0,H1}，必须生成一个缩放函数和小波，它们是由另一棵树的低通和高通滤波器生成的缩放函数和小波的近似希尔伯特变换，{G0,G1}。因此，由两棵树形成的复值标度函数和小波是近似解析的。因此，与只有2d冗余系数d- 维度数据。DTCWT 中的冗余明显小于未抽取（固定）DWT 中的冗余。该示例说明了 DTCWT 的近似平移不变性、2-D 和 3-D 中双树分析小波的选择性方向，以及双树复离散小波变换在图像和体积去噪中的使用。

DTCWT 的近移不变性

DWT 会受到偏移方差的影响，这意味着输入信号或图像中的微小偏移可能会导致 DWT 系数中信号/图像能量跨尺度的分布发生显着变化。DTCWT 近似是平移不变的。

为了在测试信号上证明这一点，构建两个长度为 128 个样本的移位离散时间脉冲。一个信号在样本 60 处具有单位脉冲，而另一个信号在样本 64 处具有单位脉冲。两个信号显然都具有单位能量（l^2 norm)）。

kronDelta1 = zeros(128,1);
kronDelta1(60) = 1;
kronDelta2 = zeros(128,1);
kronDelta2(64) = 1;

设置DWT扩展模式为周期。使用长度为 14 的小波和缩放滤波器获取两个信号低至 3 级的 DWT 和 DTCWT。提取 3 级细节系数进行比较。

origmode = dwtmode('status','nodisplay');
dwtmode('per','nodisp')
J = 3;
[dwt1C,dwt1L] = wavedec(kronDelta1,J,'sym7');
[dwt2C,dwt2L] = wavedec(kronDelta2,J,'sym7');
dwt1Cfs = detcoef(dwt1C,dwt1L,3);
dwt2Cfs = detcoef(dwt2C,dwt2L,3);

[dt1A,dt1D] = dualtree(kronDelta1,'Level',J,'FilterLength',14);
[dt2A,dt2D] = dualtree(kronDelta2,'Level',J,'FilterLength',14);
dt1Cfs = dt1D{3};
dt2Cfs = dt2D{3};

绘制第 3 级两个信号的 DWT 和 DTCWT 系数的绝对值，并计算能量（平方我2范数）的系数。在相同比例上绘制系数。信号中的四个样本移位导致了 3 级 DWT 系数的能量发生显着变化。3 级 DTCWT 系数中的能量仅改变了约 3%。

figure
tiledlayout(1,2)
nexttile
stem(abs(dwt1Cfs),'markerfacecolor',[0 0 1])
title({'DWT';['Squared 2-norm = ' num2str(norm(dwt1Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])
nexttile
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title({'DWT';['Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])

figure
tiledlayout(1,2)
nexttile
stem(abs(dt1Cfs),'markerfacecolor',[0 0 1])
title({'Dual-tree CWT';['Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])
nexttile
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title({'Dual-tree CWT';['Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])

 为了证明近似平移不变性在实际数据中的实用性，我们分析了心电图 (ECG) 信号。ECG 信号的采样间隔为 1/180 秒。数据取自 Percival & Walden [3]，第 125 页（数据最初由华盛顿大学的 William Constantine 和 Per Reinhall 提供）。为了方便起见，我们将数据从 t=0 开始。

load wecg
dt = 1/180;
t = 0:dt:(length(wecg)*dt)-dt;
figure
plot(t,wecg)
xlabel('Seconds')
ylabel('Millivolts')

 相距约 0.7 秒的大正峰值是心律的 R 波。首先，使用临界采样 DWT 和 Farras 近对称滤波器分解信号。绘制原始信号以及 2 级和 3 级小波系数。选择 2 级和 3 级系数是因为在给定采样率的这些尺度中 R 波被最显着地隔离。

figure
J = 6; 
[df,rf] = dtfilters('farras');
[dtDWT1,L1] = wavedec(wecg,J,df(:,1),df(:,2));
details = zeros(2048,3);
details(2:4:end,2) = detcoef(dtDWT1,L1,2);
details(4:8:end,3) = detcoef(dtDWT1,L1,3);
tiledlayout(3,1)
nexttile
stem(t,details(:,2),'Marker','none','ShowBaseline','off')
title('Level 2')
ylabel('mV')
nexttile
stem(t,details(:,3),'Marker','none','ShowBaseline','off')
title('Level 3')
ylabel('mV')
nexttile
plot(t,wecg)
title('Original Signal')
xlabel('Seconds')
ylabel('mV')

对二叉树变换重复上述分析。在这种情况下，只需绘制第 2 层和第 3 层双树系数的实部。

[dtcplxA,dtcplxD] = dualtree(wecg,'Level',J,'FilterLength',14);
details = zeros(2048,3);
details(2:4:end,2) = dtcplxD{2};
details(4:8:end,3) = dtcplxD{3};
tiledlayout(3,1)
nexttile
stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off')
title('Level 2')
ylabel('mV')
nexttile
stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off')
title('Level 3')
ylabel('mV')
nexttile
plot(t,wecg)
title('Original Signal')
xlabel('Seconds')
ylabel('mV')

临界采样和双树小波变换都将 ECG 波形的重要特征定位到相似的尺度

小波在一维信号中的一个重要应用是获得按尺度的方差分析。按理说，这种方差分析应该对输入信号的循环移位不敏感。不幸的是，临界采样 DWT 的情况并非如此。为了证明这一点，我们将 ECG 信号循环平移 4 个样本，使用临界采样 DWT 分析未平移和平移的信号，并计算跨尺度的能量分布。

wecgShift = circshift(wecg,4);
[dtDWT2,L2] = wavedec(wecgShift,J,df(:,1),df(:,2));

detCfs1 = detcoef(dtDWT1,L1,1:J,'cells');
apxCfs1 = appcoef(dtDWT1,L1,rf(:,1),rf(:,2),J);
cfs1 = horzcat(detCfs1,{apxCfs1});
detCfs2 = detcoef(dtDWT2,L2,1:J,'cells');
apxCfs2 = appcoef(dtDWT2,L2,rf(:,1),rf(:,2),J);
cfs2 = horzcat(detCfs2,{apxCfs2});

sigenrgy = norm(wecg,2)^2;
enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0));
enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'uni',0));
levels = {'D1';'D2';'D3';'D4';'D5';'D6';'A6'};
enr1 = enr1(:);
enr2 = enr2(:);
table(levels,enr1,enr2,'VariableNames',{'Level','enr1','enr2'})

 请注意，第 3 级和第 4 级的小波系数显示原始信号和移位信号之间的能量变化约为 3%。接下来，我们使用复杂的双树离散小波变换重复此分析。

[dtcplx2A,dtcplx2D] = dualtree(wecgShift,'Level',J,'FilterLength',14);
dtCfs1 = vertcat(dtcplxD,{dtcplxA});
dtCfs2 = vertcat(dtcplx2D,{dtcplx2A});

dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtCfs1,'uni',0));
dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtCfs2,'uni',0));
dtenr1 = dtenr1(:);
dtenr2 = dtenr2(:);
table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})

双树变换对原始信号及其循环移位版本按比例产生一致的方差分析。

图像处理中的方向选择性

二维 DWT 的标准实现使用图像的列和行的可分离过滤。使用辅助函数helperPlotCritSampDWT绘制具有 4 个消失矩的 Daubechies 最小不对称相位小波的 LH、HL 和 HH 小波sym4。\

 helperPlotCritSampDWT

请注意，LH 和 HL 小波分别具有清晰的水平和垂直方向。然而，最右侧的 HH 小波混合了 +45 度和 -45 度方向，产生了棋盘伪影。这种方向的混合是由于使用了实值可分离滤波器。HH 实值可分离滤波器在 2-D 频率平面的所有四个高频角都有通带。

双树 DWT 通过使用近似解析的小波来实现方向选择性，这意味着它们仅在频率轴的一半上有支持。在双树 DWT 中，实部和虚部都有六个子带。六个实部是通过将两棵树中的输入图像的列过滤的输出和随后的行过滤相加而形成的。六个虚部是通过减去列滤波和行滤波的输出而形成的。

应用于列和行的过滤器可以来自同一过滤器对，{H0,H1}或者{G0,G1}，或来自不同的滤波器对，{H0,G1} , {G0,H1}。使用辅助函数helperPlotWaveletDTCWT绘制与 DTCWT 实部和虚部相对应的 12 个小波的方向。上图的第一行显示了六个小波的实部，第二行显示了虚部。

helperPlotWaveletDTCWT

二维边缘表示

复杂双树小波的近似解析性和选择性方向在图像边缘表示方面提供了优于标准二维 DWT 的性能。为了说明这一点，我们使用临界采样的二维 DWT 和二维复数定向双树变换来分析边缘由多个方向上的直线和曲线奇点组成的测试图像。首先，分析由线奇点组成的八边形图像。

load woctagon
imagesc(woctagon)
colormap gray
title('Original Image')
axis equal
axis off

使用辅助函数helperPlotOctagon将图像分解为第 4 级，并根据第 4 级细节系数重建图像近似值。

helperPlotOctagon(woctagon)

 接下来，分析具有双曲边的八边形。双曲八边形中的边是曲线奇点。

load woctagonHyperbolic
figure
imagesc(woctagonHyperbolic)
colormap gray
title('Octagon with Hyperbolic Edges')
axis equal
axis off

 再次，使用辅助函数helperPlotOctagon将图像分解为第 4 级，并基于标准 2-D DWT 和面向复杂的双树 DWT 的第 4 级细节系数重建图像近似值。

helperPlotOctagon(woctagonHyperbolic)

请注意，二维临界采样 DWT 中明显的振铃伪影在两幅图像的 2-D DTCWT 中都不存在。DTCWT 更忠实地再现直线和曲线奇点。

图像去噪

由于能够在单独的子带中隔离不同的方向，因此在图像去噪等应用中，DTCWT 通常能够优于标准的可分离 DWT。为了演示这一点，请使用辅助函数helperCompare2DDenoisingDTCWT。辅助函数加载图像并添加零均值高斯白噪声σ = 25。对于用户提供的阈值范围，该函数会比较使用软阈值处理的临界采样 DWT 和 DTCWT 的去噪效果。对于每个阈值，都会显示均方根 (RMS) 误差和峰值信噪比 (PSNR)。

numex = 3;
helperCompare2DDenoisingDTCWT(numex,0:2:100,'PlotMetrics');

DTCWT 在 RMS 误差和 PSNR 方面优于标准 DWT。

接下来，获取阈值为 25 的去噪图像，该阈值等于加性噪声的标准差。

numex = 3;
helperCompare2DDenoisingDTCWT(numex,25,'PlotImage');

由于阈值等于加性噪声的标准偏差，DTCWT 提供的 PSNR 比标准 2-D DWT 高出近 4 dB

3D 方向选择性

当将小波分析扩展到更高维度时，在二维中使用可分离 DWT 观察到的振铃伪影会加剧。DTCWT 使您能够以最小的冗余保持 3-D 方向选择性。在 3-D 中，双树变换中有 28 个小波子带。

为了演示 3-D 双树小波变换的方向选择性，可视化 3-D 双树和可分离 DWT 小波的 3-D 等值面示例。首先，分别可视化两个双树子带的实部和虚部。

helperVisualize3D('Dual-Tree',28,'separate')

helperVisualize3D('Dual-Tree',25,'separate')

等值面图的红色部分表示小波从零开始的正偏移，而蓝色部分表示负偏移。您可以清楚地看到双树小波的实部和虚部在空间中的方向选择性。现在可视化双树子带之一，将实图和虚图绘制在一起作为一个等值面。

helperVisualize3D('Dual-Tree',25,'real-imaginary')

 上图表明实部和虚部在空间中是彼此移动的版本。这反映了复小波的虚部是实部的近似希尔伯特变换。接下来，以 3 维形式可视化真实正交小波的等值面以进行比较。

helperVisualize3D( 'DWT' ,7)

在 2-D DWT 中观察到的方向混合在 3-D 中甚至更加明显。正如 2-D 情况一样，3-D 方向的混合会导致明显的振铃或块效应。为了证明这一点，请检查球形体积的 3-D DWT 和 DTCWT 小波细节。球体尺寸为 64 x 64 x 64。

load sphr
[A,D] = dualtree3(sphr,2,'excludeL1');
A = zeros(size(A));
sphrDTCWT = idualtree3(A,D);
I = imtile(rescale(reshape(sphrDTCWT,[64 64 1 64])),GridSize=[8 8]);
imshow(I)
title('DTCWT Level 2 Details')

将前面的图与基于可分离 DWT 的第二级详细信息进行比较。 

sphrDEC = wavedec3(sphr,2,'sym4','mode','per');
sphrDEC.dec{1} = zeros(size(sphrDEC.dec{1}));
for kk = 2:8
    sphrDEC.dec{kk} = zeros(size(sphrDEC.dec{kk}));
end
sphrrecDWT = waverec3(sphrDEC);
figure
I = imtile(rescale(reshape(sphrrecDWT,[64 64 1 64])),GridSize=[8 8]);
imshow(I)
title('DWT Level 2 Details')

放大 DTCWT 和 DWT 蒙太奇中的图像，您将看到与 DTCWT 相比，DWT 细节中的块伪影有多么突出。

体积去噪

与 2-D 情况类似，3-D DTCWT 的方向选择性通常会导致体积去噪的改进。

为了证明这一点，请考虑由 16 个切片组成的 MRI 数据集。标准差为10的高斯噪声已添加到原始数据集中。显示噪声数据集。

load MRI3D
I = imtile(rescale(reshape(noisyMRI,[128 128 1 16])),GridSize=[4 4]);
imshow(I)

请注意，去噪之前的原始 SNR 约为 11 dB。

20*log10(norm(origMRI(:),2)/norm(origMRI(:)-noisyMRI(:),2))

答案=11.2997

使用 DTCWT 和 DWT 将 MRI 数据集去噪至 4 级。两种情况都使用类似的小波滤波器长度。将得到的 SNR 绘制为阈值的函数。显示在最佳 SNR 下获得的 DTCWT 和 DWT 的去噪结果。

[imrecDTCWT,imrecDWT] = helperCompare3DDenoising(origMRI,noisyMRI);

figure
I = imtile(rescale(reshape(imrecDTCWT,[128 128 1 16])),GridSize=[4 4]);
imshow(I)
title('DTCWT Denoised Volume')

figure
I = imtile(rescale(reshape(imrecDWT,[128 128 1 16])),GridSize=[4 4]);
imshow(I)
title('DWT Denoised Volume')

 恢复原来的扩展模式。

dwtmode(origmode,'nodisplay')

概括

我们已经证明，双树 CWT 具有临界采样 DWT 无法实现的近平移不变性和方向选择性的理想特性。我们已经演示了这些属性如何提高信号分析、图像和体积中的边缘表示以及图像和体积去噪的性能。

参考