主要内容

小波包:分解细节

这个例子展示了小波包与离散小波变换(DWT)的区别。这个例子显示了小波包变换如何导致信号的等宽子带滤波,而不是在小波变换中发现的粗倍频带滤波。

这使得小波包在许多应用中成为DWT的一个有吸引力的替代品。这里给出的两个示例是时频分析和信号分类。您必须拥有统计学和机器学习工具箱™ 信号处理工具箱™ 进行分类。

如果将正交小波与小波包变换结合使用,则会在等宽子带之间对信号能量进行分割。

这个例子集中在一维情况下,但是许多概念直接扩展到图像的小波包变换。

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

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

图1:对于一维输入信号,DWT树降至3级。

G f 缩放(低通)分析过滤器和 H f 表示小波(高通)分析滤波器。底部的标签显示频率轴[0,1/2]划分为子带。

该图显示DWT的后续电平仅在低通(缩放)输出上运行滤波器。在每一级,DWT将信号分为倍频带。在临界采样DWT中,带通滤波器的输出在每一级下采样两次。在非抽样离散小波变换中,输出不下采样。

将小波变换树与全小波包树进行比较。

图2:完整的小波包树下降到第3级。

在小波包变换中,滤波操作也应用于小波或细节系数。结果是,小波包将输入信号的子带滤波成逐渐变细的等宽间隔。在每个级别, j ,频率轴[0,1/2]分为 2 j 子频带。电平上的子频带(以赫兹为单位) j 大约

n F 年代 2 j + 1 n + 1 F 年代 2 j + 1 n 0 1 ... 2 j - 1 在哪里 财政司司长 为采样频率。

这个带通近似有多好取决于滤波器的频率本地化程度。比如Fejér-Korovkin过滤器“fk18”(18个系数),近似值相当好。对于像Haar这样的滤波器(“哈尔”),近似是不准确的。

在临界采样小波包变换中,带通滤波器的输出由两个降采样。在非降采样小波包变换中,输出不降采样。

小波包时频分析

由于小波包将频率轴划分成比DWT更精细的间隔,小波包在时频分析中具有优越性,作为一个例子,在加性噪声中考虑两个频率为150和200 Hz的间歇性正弦波,在1 kHz处采样数据,以防止临界SAMP中固有的时间分辨率的丢失。led小波包变换,采用非抽样变换。

dt = 0.001;t = 0: dt: 1 dt;x =...因为(2 *π* 150 * t) * (t > = 0.2 & t < 0.4) +罪(2 *π* 200 * t) * 0.6 (t > & t < 0.9);y = x + 0.05 * randn(大小(t));(wpt - F) = modwpt (x,“TimeAlign”,对);等高线(t,F.*(1/dt),abs(wpt)。^2)网格xlabel(的时间(秒)) ylabel (“赫兹”)头衔(《时频分析——非抽取小波包变换》

图中包含一个轴对象。以时频分析—非抽取小波包变换为标题的轴对象包含一个轮廓型对象。

注意小波包变换能够分离150hz和200hz分量。这对于DWT是不正确的,因为150和200赫兹落在同一个八度频带。四级DWT的八度频带为(赫兹)。

  • [0,31.25)

  • [31.25,62.5)

  • [62.5,125)

  • [125,250)

  • [250,500)

证明这两个分量是被小波变换时定域的,但在频率上不是分离的。

mra=modwtmra(modwt(y,“fk18”,4),“fk18”);J = 5: 1:1;频率= 1/2 * (1000. / 2 ^ J);轮廓(t,频率,flipud (abs (mra)。^ 2))网格ylim([0 500])xlabel(的时间(秒)) ylabel (“赫兹”)头衔(“时频分析—非抽取小波变换”

图中包含一个axes对象。标题为“时频分析-未抽取小波变换”的axes对象包含一个contour类型的对象。

当然,连续小波变换(CWT)提供了比小波包变换更高分辨率的时频分析,但小波包具有作为正交变换的附加优点(当使用正交小波时)正交变换意味着信号中的能量被保留,并在系数之间进行分割,如下一节所示。

小波包变换中的能量保持

除了在每个级别将信号过滤为等宽子带外,小波包变换还在子带之间划分信号的能量。对于抽取和非抽取小波包变换都是如此。为了证明这一点,加载包含地震记录的数据集。该数据类似于ti下面的信号分类示例中使用的me系列。

负载科比绘图(神户)网格xlabel(“秒”)头衔(“神户地震数据”)轴心牢固的

图中包含一个轴对象。标题为“Kobe地震数据”的axes对象包含一个line类型的对象。

获取数据的抽取和非抽取小波包变换,直至第3级。为确保抽取小波包变换的结果一致,请将边界扩展模式设置为“周期”

wptreeDecimated =方法(科比,“水平”3.“边界”“周期”); wptUndecimated=modwpt(神户,3);

计算抽取和非抽取小波包三级系数中的总能量,并与原始信号中的能量进行比较。

抽取能量=sum(cell2mat(cellfun)(x)sum(abs(x))^2,wptre抽取,“UniformOutput”,false)))
抽取能量=2.0189e+11
undecimatedEnergy =总和(总和(abs (wptUndecimated)。^ 2,2))
undecimatedEnergy = 2.0189 e + 11
信号能量=标准(神户,2)^2
信号能量=2.0189e+11

小波变换与小波包变换具有相同的重要性质。

wt = modwt(科比,“fk18”3);undecimatedWTEnergy =总和(总和(abs (wt)。^ 2,2))
undecimatedWTEnergy = 2.0189 e + 11

由于每个级别的小波包变换将频率轴划分为等宽的间隔,并在这些间隔之间划分信号能量,因此通常只需保留每个小波包中的相对能量即可构造有用的特征向量进行分类。下一个示例演示了这一点。

小波包分类——地震还是爆炸?

地震记录从许多来源提取活动。能够根据其来源对活动进行分类非常重要。例如,地震和爆炸可能呈现类似的时域特征,但区分这两个事件非常重要。数据集EarthQuakeExplosionData包含16个录音和2048个样本。前8个记录(栏)是地震的记录,后8个记录(栏)是爆炸的地震记录。载入数据并绘制地震和爆炸记录以作比较。数据取自R包astsa Stoffer[4]这是沙姆韦(Shumway)和斯托夫(Stoffer)的补充[3]. 作者欣然同意在这个例子中使用它。

绘制每组的时间序列以作比较。

负载EarthQuakeExplosionData子地块(2,1,1)图(EarthQuakeExplosionData(:,3))xlabel(“时间”)头衔(“地震记录”)网格牢固的次要情节(2,1,2)情节(EarthQuakeExplosionData(:, 9)包含(“时间”)头衔(“爆炸记录”)网格牢固的

图中包含2个轴对象。标题为地震记录的轴对象1包含line类型的对象。标题为爆炸记录的轴对象2包含line类型的对象。

通过使用“fk6”具有非抽样小波包变换的小波。这将产生8个子带,宽度约为1/16个循环/样本。使用每个子带中的相对能量创建特征向量。作为示例,获得第一次地震记录的小波包相对能量特征向量。

(wpt ~ F E RE) = modwpt (EarthQuakeExplosionData (: 1), 3,“fk6”);

再保险包含了8个子带中的每个子带的相对能量。如果你把所有的元素加起来再保险,等于1。请注意,这是数据的显著减少。长度为2048的时间序列减少为一个包含8个元素的向量,每个元素表示3级小波包节点中的相对能量。辅助函数helperEarthQuakeExplosionClassifier获得每16个记录的相对能量在水平3使用“fk6”小波。这将产生一个16×8的矩阵,并使用这些特征向量作为kmeans分类器的输入。分类器使用间隙统计准则来确定1到6之间的特征向量的最佳聚类数,并对每个向量进行分类。此外,分类器对使用该方法获得的第3级未抽取小波变换系数执行完全相同的分类“fk6”每个时间序列的小波和功率谱。非抽取小波变换的结果是一个16乘4的矩阵(3个小波子带和1个缩放子带)。功率谱的结果是一个16乘1025的矩阵。的源代码helperEarthQuakeExplosionClassifier附录中列出了。您必须拥有统计和机器学习工具箱™ 信号处理工具箱™ 运行分类器。

级别=3;[WavPacketCluster、WtCluster、PSpeCluster]=...helperEarthQuakeExplosionClassifier (EarthQuakeExplosionData级别)
WavPacketCluster = GapEvaluation with properties: NumObservations: 16 InspectedK: [1 2 3 4 5 6] criteria values: [-0.0039 0.0779 0.1314 0.0862 0.0296 -0.0096] OptimalK: 2
WtCluster=GapEvaluation with properties:NumObservations:16检查K:[1 2 3 4 5 6]标准值:[-0.0095 0.0474 0.0706 0.0260-0.0280-0.0346]最佳值:2
pspecluster=GapEvaluation with properties:numobervations:16检查k:[1 2 3 4 5 6]标准值:[0.3804 0.3574 0.3466 0.2879 0.3256 0.3326]最佳k:1

波包簇为未抽取小波包特征向量的聚类输出。WtCluster聚类输出是否使用未抽取DWT特征向量,和PspecCluster为基于功率谱的聚类输出。小波包分类识别出两个簇(K:2)为最优数。检验小波包聚类的结果。

res = [WavPacketCluster.OptimalY(1:8)';WavPacketCluster.OptimalY(9:结束)')
res=2×81 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 2

你可以看到,只有两个错误被发现。在八次地震记录中,有七次被归为一类(第二组)一个被错误地归类为第1组。8个爆炸记录也是如此——7个被归类在一起,1个被错误归类。基于三级未抽样DWT和功率谱的分类返回一个簇作为最佳解决方案。

如果你重复上面的水平等于4,非抽取小波变换和小波包变换执行相同。两者都将两个聚类定义为最优,都将第一个和最后一个时间序列错误地归为另一个组。

结论

在本例中,您了解了小波包变换与离散小波变换的不同之处。具体而言,小波包树还对小波(细节)系数进行过滤,而小波变换仅对缩放(近似)系数进行迭代。

您了解到小波包变换具有小波变换的重要能量保持特性,同时提供更高的频率分辨率。在某些应用中,这种优越的频率分辨率使小波包变换成为一种有吸引力的替代小波变换。

附录

helperEarthQuakeExplosionClassifier

作用[WavPacketCluster,WtCluster,pspecluster]=帮助器EarthQuakeExplosionClassifier(数据,级别)%此函数仅用于支持特征示例万博1manbetx% waveletpacketsdemo.mlx。这个函数可能会在未来的版本中被更改或删除。% Data -以单个时间序列为列向量的数据矩阵。% Level -小波包和小波变换的级别NP = 2 ^水平;西北= + 1级;NP WpMatrix = 0(16日);西北WtMatrix = 0(16日);对于jj = 1:8 [~,~,~, RE] = modwpt(数据(:,jj),水平,“fk6”);wt=modwt(数据(:,jj),电平,“fk6”);WtMatrix (jj:) =总和(abs (wt)。^ 2,2)。/规范(数据(:,jj), 2) ^ 2;WpMatrix (jj:) =再保险;终止对于kk=9:16[,,,,,,,RE]=modwpt(数据(:,kk),电平,“fk6”);wt=modwt(数据(:,kk),电平,“fk6”)矩阵(kk,:)=sum(abs(wt)。^2,2)。/norm(Data(:,kk),2)^2;矩阵(kk,:)=RE;终止%的可重复性rng (“默认”);WavPacketCluster = evalclusters (WpMatrix,“kmeans”“差距”“中”1:6,...“距离”“cityblock”);WtCluster=评估聚类(wt矩阵,“kmeans”“差距”“中”1:6,...“距离”“cityblock”);Pxx =周期图(数据,汉明(元素个数(数据(:1))),[],1,“权力”); pspecluster=evalclusters(Pxx’,“kmeans”“差距”“中”1:6,...“距离”“cityblock”);终止

参考文献

[1] 维克豪斯,姆拉登·维克托。从理论到软件的自适应小波分析.马萨诸塞州韦尔斯利:A.K.彼得斯,1994年。

Percival, Donald B.和Andrew T. Walden。时间序列分析的小波方法.统计和概率数学的剑桥系列。剑桥 ;纽约:剑桥大学出版社,2000。

[3] 罗伯特·H·舒姆韦和大卫·S·斯托弗。时间序列分析及其应用:附R个例子第三版,《统计学中的斯普林格文本》。纽约:斯普林格,2011年。

Stoffer, d.h。asta:应用统计时间序列分析。R包版本1.3。https://CRAN.R-project.org/package=astsa, 2014.

另见

|||