Wavelet Packets: Decomposing the Details

学习笔记：

本例展示了小波包与离散小波变换（DWT）的不同之处。该示例展示了小波包变换如何对信号进行等宽子带滤波，而不是 DWT 中更粗糙的倍频程带滤波。
这使得小波包在许多应用中成为 DWT 的替代品。这里介绍的两个例子是时频分析和信号分类。您必须拥有 Statistics and Machine Learning Toolbox™ 和 Signal Processing Toolbox™ 才能进行分类。如果使用小波包变换的正交小波，还可以在等宽子带之间对信号能量进行分区。
本示例侧重于一维情况，但许多概念可直接扩展到图像的小波包变换。

离散小波变换和离散小波包变换

下图显示了一维输入信号的 DWT 树。

 图 1：1-D 输入信号的 DWT 树，下至第 3 级。
 G表示缩放（低通）分析滤波器，H表示小波（高通）分析滤波器。底部的标签表示将频率轴 [0,1/2] 划分为子带。
从图中可以看出，后续各级 DWT 只对低通（缩放）滤波器的输出进行处理。在每一级，DWT 都会将信号划分为倍频程带。在临界采样的 DWT 中，带通滤波器的输出在每一级都被降低两个采样率。在未衰减离散小波变换中，输出不进行降采样。
将 DWT 树与全小波包树进行比较。

图 2：完整的小波包树直到第 3 层。
在小波包变换中，滤波操作也应用于小波系数或细节系数。其结果是，小波包将输入信号过滤成逐渐细化的等宽区间。在每一级，频率轴 [0,1/2] 被划分为子带。以赫兹为单位的子带级别约为 这种带通近似的效果如何，取决于滤波器的频率定位程度。对于像 "fk18"（18 个系数）这样的 Fejér-Korovkin 滤波器，近似效果相当不错。而对于像哈尔（'haar'）这样的滤波器，近似值并不准确。在严格采样的小波包变换中，带通滤波器的输出被降低了两个采样率。而在未衰减小波包变换中，输出不进行降采样。

使用小波包进行时频分析

由于小波包将频率轴划分为比 DWT 更细的区间，因此小波包在时频分析方面更胜一筹。举例来说，考虑两个频率分别为 150 和 200 Hz 的间歇性正弦波。数据采样频率为 1 kHz。为了避免临界采样小波包变换固有的时间分辨率损失，可以使用未衰减变换。

dt = 0.001;
t = 0:dt:1-dt;
x = ...
cos(2*pi*150*t).*(t>=0.2 & t<0.4)+sin(2*pi*200*t).*(t>0.6 & t<0.9);
y = x+0.05*randn(size(t));
[wpt,~,F] = modwpt(x,'TimeAlign',true);
contour(t,F.*(1/dt),abs(wpt).^2)
grid on
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Packet Transform')

请注意，小波包变换能够分离 150 赫兹和 200 赫兹的分量。而 DWT 则不然，因为 150 和 200 Hz 属于同一个倍频程带。四级 DWT 的倍频程带为（单位：赫兹）

mra  = modwtmra(modwt(y,'fk18',4),'fk18');
J = 5:-1:1;
freq = 1/2*(1000./2.^J);
contour(t,freq,flipud(abs(mra).^2))
grid on
ylim([0 500])
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Transform')

当然，与小波包变换相比，连续小波变换（CWT）能提供更高分辨率的时频分析，但小波包还有一个额外的好处，那就是它是一种正交变换（当使用正交小波时）。正交变换意味着信号中的能量得以保留，并在各系数之间进行分配，正如下一节所演示的那样。

小波包变换中的能量保护

小波包变换除了将信号过滤成各级等宽的子带外，还将信号的能量划分到各个子带中。小波包变换的抽取和不抽取都是如此。为了演示这一点，请加载一个包含地震记录的数据集。该数据与下面信号分类示例中使用的时间序列类似。

load kobe
plot(kobe)
grid on
xlabel('Seconds')
title('Kobe Earthquake Data')
axis tight

同时获取数据的小波包变换（最小到第 3 层）的去微分和未去微分。为确保小波包十进制变换结果的一致性，请将边界扩展模式设置为 "周期"。

wptreeDecimated = dwpt(kobe,'Level',3,'Boundary','periodic');
wptUndecimated = modwpt(kobe,3);

 计算小波包三级系数中经细分和未细分的总能量，并与原始信号中的能量进行比较。

decimatedEnergy = sum(cell2mat(cellfun(@(x) sum(abs(x).^2),wptreeDecimated,'UniformOutput',false)))

 

undecimatedEnergy = sum(sum(abs(wptUndecimated).^2,2))

 

signalEnergy = norm(kobe,2)^2

DWT 与小波包变换共享这一重要特性。

wt = modwt(kobe,'fk18',3);
undecimatedWTEnergy = sum(sum(abs(wt).^2,2))

由于每一级的小波包变换都会将频率轴划分为等宽的区间，并将信号能量划分到这些区间中，因此通常只需保留每个小波包中的相对能量，就能构建出有用的特征向量进行分类。下一个示例演示了这一点。

小波包分类 -- 地震还是爆炸？

地震记录会捕捉到许多来源的活动。根据其来源对这些活动进行分类非常重要。例如，地震和爆炸可能呈现类似的时域特征，但区分这两种事件非常重要。数据集 EarthQuakeExplosionData 包含 16 条记录，共 2048 个样本。前 8 个记录（列）是地震的地震记录，后 8 个记录（列）是爆炸的爆炸记录。加载数据并绘制地震和爆炸记录对比图。数据来自 R 软件包 astsa Stoffer 4，它是 Shumway 和 Stoffer 3 的补充。作者允许在本示例中使用该数据包。
绘制每组数据的时间序列图进行比较。

load EarthQuakeExplosionData
subplot(2,1,1)
plot(EarthQuakeExplosionData(:,3))
xlabel('Time')
title('Earthquake Recording')
grid on
axis tight
subplot(2,1,2)
plot(EarthQuakeExplosionData(:,9))
xlabel('Time')
title('Explosion Recording')
grid on
axis tight

使用 "fk6 "小波对每个时间序列进行分解到第三级，并进行未递减的小波包变换，从而形成小波包特征向量。这将产生 8 个子带，宽度约为 1/16 个周期/样本。利用每个子带的相对能量创建特征向量。例如，获取第一次地震记录的小波包相对能量特征向量。

[wpt,~,F,E,RE] = modwpt(EarthQuakeExplosionData(:,1),3,'fk6');

 

 RE 包含 8 个子带中每个子带的相对能量。如果将 RE 中的所有元素相加，等于 1。请注意，这大大减少了数据量。长度为 2048 的时间序列被缩减为一个包含 8 个元素的向量，每个元素代表第 3 层小波包节点中的相对能量。辅助函数 helperEarthQuakeExplosionClassifier 使用 "fk6 "小波获取第 3 层 16 个记录中每个记录的相对能量。这样就得到了一个 16 乘 8 的矩阵，并将这些特征向量作为 kmeans 分类器的输入。分类器使用 Gap 统计标准来确定特征向量的最佳簇数（1-6），并对每个向量进行分类。此外，分类器还对使用 "fk6 "小波和功率谱获得的每个时间序列的第 3 级未衰减小波变换系数执行完全相同的分类。未估计小波变换产生一个 16 x 4 矩阵（3 个小波子带和 1 个缩放子带）。功率谱产生一个 16 x 1025 的矩阵。helperEarthQuakeExplosionClassifier 的源代码见附录。您必须安装统计和机器学习工具箱（Statistics and Machine Learning Toolbox™）和信号处理工具箱（Signal Processing Toolbox™）才能运行分类器。

Level = 3;
[WavPacketCluster,WtCluster,PspecCluster] = ...
helperEarthQuakeExplosionClassifier(EarthQuakeExplosionData,Level)

WavPacketCluster 是未细分小波包特征向量的聚类输出。WtCluster 是使用未细分 DWT 特征向量的聚类输出，PspecCluster 是基于功率谱的聚类输出。小波包分类将两个聚类（OptimalK：2）确定为最佳数量。检查小波包聚类的结果。

res = [ WavPacketCluster.OptimalY(1:8)' ; WavPacketCluster.OptimalY(9:end)']

 

您可以看到只有两个错误。在 8 个地震记录中，7 个被归为一组（第 2 组），1 个被错误地归为第 1 组。 8 个爆炸记录也是如此 -- 7 个被归为一组，1 个被错误地归为第 1 组。基于第三级未细分 DWT 和功率谱的分类返回一个群组作为最优解。
如果重复上述操作，将级别设为 4，则未估计小波变换和小波包变换的表现完全相同。两者都将两个集群确定为最优集群，并且都将第一个和最后一个时间序列错误地归类为属于另一组。

结论

通过本例，您了解了小波包变换与离散小波变换的不同之处。具体来说，小波包树还会过滤小波（细节）系数，而小波变换只会迭代缩放（近似）系数。
您将了解到，小波包变换与 DWT 一样具有重要的能量保存特性，同时还能提供出色的频率分辨率。在某些应用中，这种出色的频率分辨率使得小波包变换成为 DWT 的一个有吸引力的替代方案。

function [WavPacketCluster,WtCluster,PspecCluster] = helperEarthQuakeExplosionClassifier(Data,Level)  
%   This function is provided solely in support of the featured example
%   waveletpacketsdemo.mlx. 
%   This function may be changed or removed in a future release.
%   Data - The data matrix with individual time series as column vectors.
%   Level - The level of the wavelet packet and wavelet transforms

NP = 2^Level;
NW = Level+1;
WpMatrix = zeros(16,NP);
WtMatrix = zeros(16,NW);

for jj = 1:8
    [~,~,~,~,RE] = modwpt(Data(:,jj),Level,'fk6');
    wt = modwt(Data(:,jj),Level,'fk6');
    WtMatrix(jj,:) = sum(abs(wt).^2,2)./norm(Data(:,jj),2)^2;
    WpMatrix(jj,:) = RE;
end

for kk = 9:16
    [~,~,~,~,RE] = modwpt(Data(:,kk),Level,'fk6');
    wt = modwt(Data(:,kk),Level,'fk6');
    WtMatrix(kk,:) = sum(abs(wt).^2,2)./norm(Data(:,kk),2)^2;
    WpMatrix(kk,:) = RE;
end

% For repeatability
rng('default');

WavPacketCluster = evalclusters(WpMatrix,'kmeans','gap','KList',1:6,...
    'Distance','cityblock');

WtCluster = evalclusters(WtMatrix,'kmeans','gap','KList',1:6,...
'Distance','cityblock');

Pxx = periodogram(Data,hamming(numel(Data(:,1))),[],1,'power');

PspecCluster = evalclusters(Pxx','kmeans','gap','KList',1:6,...
    'Distance','cityblock');
end