第8章 小波包及其应用

从前面的章节已经了解到小波变换是将时频分析联系在一起的分析方法，而多分辨分析又是将信号分为高频和低频两个部分的分析方法。然而，对于低频部分时常保留，而高频部分的分析方法却往往不知所措。
小波包相对于小波的主要优点是小波包可以对信号的高频部分做更加细致地刻画，对信号的分析能力更强，因此掌握小波包分析的方法也尤为重要。
在第一节将会重点介绍什么是小波包以及小波包的性质。紧接着将介绍小波包在MATLAB中的应用，主要包括对信号和图像的压缩与去噪处理。特别是，详细地描述了小波包分解与重构函数的使用及算例。读者可以根据自身需要来选择自己的学习尺度。
学习目标：
(1)了解小波包的定义及性质
(2)熟悉一维小波包分解和重构函数
(3)熟悉二维小波包分解和重构函数
(4)熟悉小波包界面式应用
8.1 小波包
小波包分析能够为信号提供一种更精细的分析方法，它将频带进行多层次划分，对多分辨分析没有细分的高频部分进一步分解。
8.1.1 小波包的定义
短时傅立叶变换对信号的频带划分是线性等间隔的。多分辨分析可以对信号进行有效的时频分解，但由于其尺度是按二进制变化的，所以在高频频段其频率分辨率较差，而在低频频段其时间分辨率较差，即对信号的频带进行指数等间隔划分。
小波包分析能够为信号提供一种更精细的分析方法，它将频带进行多层次划分，对多分辨分析没有细分的高频部分进一步分解，并能够根据被分析信号的特征，自适应地选择相应频带，使之与信号频谱相匹配，从而提高了时—频分辨率，因此小波包具有更广泛的应用价值。
关于小波包分析的理解，这里以一个二层的分解进行说明，其小波包分解树如图8-1所示。图8-1中A表示低频，D表示高频，末尾的序号数表示小波分解的层树(也即尺度数)。
在多分辨分析中，
 ，表明多分辨分析是按照不同的尺度因子j把Hilbert空间
 分解为所有子空间Wj
 (j∈Z)的正交和。其中，Wj
 为小波函数Ψ(t)的闭包(小波子空间)。现在，我们希望对小波子空间Wj
 按照二进制进行频率的细分，以达到提高频率分辨率的目的。
图8-1 小波包分解树
一种自然的做法是将尺度空间Vj
 和小波子空间Wj
 用一个新的子空间
 统一起来表征，若令
则Hilbert空间的正交分解Vj+1
 =Vj
 ⊕Wj
 即可用 的分解统一为：
定义子空间
 是函数Un
 (t)的闭包空间，而Un
 (t)是函数U2n
 (t)的闭包空间，并令Un
 (t)满足下面的双尺度方程：
式中，
 ，即两系数也具有正交关系。当n=0时，以上两式直接给出
与在多分辨分析中，ϕ(t)、ψ(t)满足双尺度方程：
相比较，u0
 (t)和u1
 (t)分别退化为尺度函数ϕ(t)和小波基函数ψ(t)。式(8.3)是式(8.1)的等价表示。把这种等价表示推广到n∈Z+
 (非负整数)的情况，即得到式(8.2)的等价表示为：
定义(小波包)由式(8.2)构造的序列{un
 (t)}(其中n∈Z+
 )称为由基函数u0
 (t)=ϕ(t)确定的正交小波包。当n=0时，即为式(8.3)的情况。
由于ϕ(t)由hk
 唯一确定，所以又称{un
 (t)}n∈Z
 为关于序列{hk
 }的正交小波包。
8.1.2 小波包的性质
定理 8-1 设非负整数 n的二进制表示为
 ，εi
 =0或1，则小波包 的傅立叶变换由下式给出：
式中：
定理8-2 设{un
 (t)}n∈Z
 是正交尺度函数ϕ(t)的正交小波包，则
 即{un
 (t)}n∈Z
 构成L2
 (R)的规范正交基。
8.1.3 小波包的空间分解
令{un
 (t)}n∈Z
 是关于hk 的小波包族，考虑用下列方式生成子空间族。现在令n=1，2，…；j=1，2，…，并对式(8.1)作迭代分解，则有：
因此，我们很容易得到小波子空间Wj
 的各种分解如下：
Wj
 空间分解的子空间序列可写作 子空间序列
 的标准正交基为
容易看出，当l=0和m=0时，子空间序列
 简化为
 ，相应的正交基简化为
 ，它恰好是标准正交小波族
若 n是一个倍频程细划的参数，即令 n=2l
 +m，则我们有小波包的简略记号
我们把ψj，k，n
 (t)称为既有尺度指标j、位置指标k和频率指标n的小波包。将它与前面的小波ψj，k，
 (t)作一比较知，小波只有离散尺度j和离散平移k两个参数，而小波包除了这两个离散参数外，还增加了一个频率参数n=2l
 +m。
正是这个频率新参数的作用，使得小波包克服了小波时间分辨率高时频率分辨率低的缺陷，于是，参数 n 表示
 函数的零交叉数目，也就是其波形的震荡次数。
定义(小波库)由ψn
 (t)生成的函数族ψj，k，n
 (t) (其中n∈Z+
 ；j，k∈Z)称为由尺度函数ψ(t)构造的小波库。
推论 对于每个j=0，1，2，…
这时，族
是L2
 (R)的一个正交基。
随着尺度 j的增大，相应正交小波基函数的空间分辨率越高，而其频率分辨率越低，这正是正交小波基的一大缺陷。而小波包却具有将随 j 增大而变宽的频谱窗口进一步分割变细的优良性质，从而克服了正交小波变换的不足。
小波包可以对Wj
 进一步分解，从而提高频率分辨率，是一种比多分辨分析更加精细的分解方法，具有更好的时频特性。
8.1.4 小波包算法
下面给出小波包的分解算法和重构算法。设
 可表示为：
小波包分解算法：由
小波包重构算法：由
8.2 一维小波包在MATLAB中的应用
本节将上一节的理论知识转化成实现形式，通过MATLAB提供的一些具体函数和应用实例来了解小波包的实际用途。
8.2.1 一维小波包函数
上一节已经了解到小波包分析分解中，细节部分也进行相同的分解。小波包分解具有任意多尺度特点，避免了小波变换固定时频分解的缺陷(如高频段频率分辨率低)，为时频分析提供了极大的选择余地，更能反映信号的本质和特征。
(1)小波包一维分解函数wpdec，其调用格式为：
T=wpdec(X，N，'wname'，E，P)
T=wpdec(X，N，'wname')
其中：X是信号向量；N是分解层数；wname是小波基；E是熵的类型；P是可选参数，取决于返回值T；T是返回的小波包树。
(2)小波包能量谱wenergy，其调用格式为：
E=wenergy(T)
其中：T是小波包树形结构；E代表能量谱。
(3)小波包的一维重构函数wprec，其调用格式为：
X=wprec(T)
其中：T是小波包树形结构；X是信号向量。
【例 8.1】利用小波包分析函数，将噪声信号分解为三层的小波包树，并且选取 Haar作为小波基，熵的类型为shannon。
MATLAB代码设置如下：
load noisdopp;
x=noisdopp;
wpt=wpdec(x，3，'db1'，'shannon');
plot(wpt)
画出的结果如图8-2所示。
图8-2 小波包树【例8.2】利用小波包分解函数及能量谱函数，将噪声信号分别分解为第一到第四层，并且分析每层的能量谱。MATLAB代码设置如下：
load noisdopp;
s=noisdopp;
% 画出初始信号
plot(s)
for i=1:4
wpt=wpdec(s，i，'db1');
e=wenergy(wpt);
E=zeros(1，length(e));
for j=1:2^i
E(j)=sum(abs(wprcoef(wpt，[i，j-1])).^2);
end
% 画出第一到第四层分解图的能量谱
figure(i+1)
bar(e);
axis([0 length(e) 0 130]);
title(['第'，num2str(i)，'层']);
for j=1:length(e)
text(j-0.2，e(j)+20，num2str(e(j)，'%4.2f'));
end
end
结果分别如图8-3和图8-4所示。
图8-3 噪声信号
图8-4 第一层到第四层分解的能量谱
图8-4 第一层到第四层分解的能量谱(续)
【例8.3】利用小波包分析函数，将噪声信号分解为三层的小波包树，并且选取db3作为小波基，熵的类型为shannon。最后，利用小波包重构函数使其重构信号，并且计算绝对误差值。
MATLAB代码设置如下：
加载信号：
load noisdopp;
s=noisdopp;
plot(s)
title('原始信号')
结果如图8-5所示。
t=wpdec(s，3，'db3'，'shannon');
x=wprec(t);
figure(2)
plot(x)
title('重构信号')
结果如图8-6所示。
图8-5 原始信号
图8-6 重构信号
err=abs(s-x);
figure(3)
plot(err)
title('绝对误差')
结果如图8-7所示。
图8-7 原始与重构信号的绝对误差
8.2.2 一维小波包界面式应用——信号压缩
前一小节主要介绍了小波包分析函数进行小波包分析，MATLAB又提供了图形接口界面，让使用者更加简单方便地通过按键来进行小波包分析。本小节将通过一个具体的实例介绍这些操作流程。
(1)启动小波工具箱。
启动小波工具箱有两种方法：第一方法是在MATLAB窗口中键入wavemenu命令；第二种方法是单击APPS，然后找到wavelet单击即可，具体如图8-8所示。
按照上述两种方法操作会弹出小波分析界面。
图8-8 小波工具箱启动
说明：在小波工具箱主界面中分别有一维工具、二维工具、三维工具、一维多重工具、显示工具、一维专用工具、二维专用工具、小波设计工具及扩展工具等。每种工具中又分为各种具体类型的应用，例如一维工具箱中包含有一维小波、一维小波包、一维连续小波及一维复连续小波。
(2)启动一维小波包的主界面。
单击Wavelet Packet 1-D按钮，弹出一维小波包的主界面。
(3)下载信号源。
选择 File→Load Signal命令，在弹出的对话框中，选择 noissin.mat 文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到了此工具中，如图8-9所示。
图8-9 加载信号
从图8-9中可以看出，一维小波包图形工具窗口在加载完信号后的工具界面，左边是信号的显示区域，右边是对信号进行小波包分析的各种按键和参数选择设置框。在左边的显示区域中有4个图形区如下。
①Analyzed Signal 在右上角，显示原始信号波形。
②Decomposition Tree 在左上角，显示小波包分解树结构。右移上端的滚动条可以放大图形，左移滚动条则缩小图形。
③Node Action Result 在左下角，显示某个节点分析的结果。
④Colored Coefficients for Terminal Nodes在右下角，显示信号小波包分解后系数的灰度图像。在一维小波包图形工具中，对一个信号是按照完整二叉树的方式来进行n层的小波包分解的，这样终节点有2n
 个。图形工具对第n层的2n
 组分解系数以灰度图像的方式来显示其结果。在其下方还有一个标为Scale of Colors from MIN to MAX的颜色条，用来表示系数大小与颜色的对应关系。
(4)小波包分析。
在图8-9的右边可以对信号进行设置。选择haar小波，并且分角层数设置为3，选择的熵类型为shannon。单击Analyze按钮，执行对信号分解，显示在窗口中，如图8-10所示。
图8-10 分析信号
这里对Entorpy(熵标准)与上一节介绍的wpdec小波分解函数中的输入参数是一致的。在MATLAB中，可选用的熵标准类型如下。
①Shannon熵：

 (约定0log0=0)。
②Threshold阈值熵：设阈值为ε，如果
 。所以此时E(s)等于信号所有采样值大于阈值ε的采样点个数。
③nor范数熵。
④logenerge对数能量熵：
⑤sure 熵：类似阈值熵，只不过这里的阈值
 ，n为信号长度。
⑥user用户熵：用户以M文件格式自己定义的熵。
关于熵更详细的信息可以利用系数帮助Help来查看计算小波包的熵的函数用法。
(5)计算最佳树。
在图8-10中的右上方单击Best Tree按钮，其功能就是在执行信号的小波包分解之后用熵选择框确定的熵标准计算出最佳的小波包分解树。在它上方还有一个Initial Tree按钮，这是个初始树显示。用户可以利用这两个按钮对比，如图8-11所示。
图8-11 初始树与最佳树对比
(6)选择阈值进行压缩。
单击Compress按钮，弹出如图8-12所示的一维小波包压缩窗口。
在阈值框里可以键入用户自己设定的阈值，同时观察压缩后保留的能量与置零的系数百分比的变换。
图8-12 一维小波包压缩窗口
再单击View Compress Signal按钮可以看到压缩信号与初始信号的对比图，如图8-13所示。
图8-13 压缩信号与初始信号的对比图
可以看出，压缩信号(黑色)和原始信号(红色)在同一个显示框中。
说明：此功能是在信号小波分解后对各层的高频系数进行阈值量化处理，然后再根据处理后的系数进行信号重构，从而达到对信号用小波包方法进行压缩的目的。
注：在一维小波包图形工具中，只有全局阈值。
8.2.3 一维小波包界面式应用——信号去噪
本小节将采用一维小波包图形工具来对一个噪声信号进行去噪。
(1)启动小波工具箱。方法同上一小节。
(2)启动一维小波包的主界面。方法同上一小节。
(3)下载信号源。
选择File→Load Signal命令，在弹出的对话框中选择noisbloc.mat文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到了此工具中，如图8-14所示。
图8-14 加载信号
计算时采用SURE熵标准。
(4)小波包分析。
在图8-14的右边可以对信号进行设置。选择db3小波，并且分角层数设置为3，选择的熵类型为sure，阈值框内设置为5.82。单击Analyze按钮，执行对信号分解，显示在窗口中，如图8-15所示。
注意：图8-15的右边有4组列表框，分别介绍如下。
①Cut Tree at Level 列表框。
此项用于将小波包分解树剪切到某一层。如果对某个信号进行若干层分解时，该选择框里的数字从0到n。如果该信号只需要分解到第k层就可以满足分析要求，则在该列表框里输入数字k即可。
图8-15 信号分析
②Node Label 列表框。
此项用于将节点标志方式，共有7种表示方式，如下。
•Depth-Pos(深度-位置)方式：采用二维数组进行标注，如(1，2)等。
•Index(索引)方式：以节点的索引值来标志节点，根节点索引值为0，然后从上到下、从左到右依次为1，2，3，….
•Entropy(熵值)方式：在各个节点处标志其对应的熵值。
•OptEnt(最优熵)方式：根据该节点是否是最优熵来标志。
•Length 方式：标志各节点的数据长度。
•None是无标志方式。
•Type(类型)方式：根据各节点是高频还是低频相应标注。
③Node Action 列表框。
此项用于对节点操作进行选择，共有7种方式，如下。
•Visualize(显示)方式时，单击左边树结构图中的任意一节点，则在左边的 Node Action Result显示框中就会显示该节点的系数曲线。
•Split/Merge(分解/合成)方式时，单击树结构中的任一节点，如果该节点是终节点，则对该节点进行下一层的小波包单层二叉分解；如果该节点不是终节点，就将该节点以下的所有子节点剪切。
•Recons(重构)方式时，对所选节点的系数进行重构，使重构后的长度等于原始信号长度。
•Select On(选中)方式时，选中树结构中的几个节点，单击出现在该列表框下的Reconstruct按钮，可以对这几个节点的分解系数进行重构。
 

•Select off(取消)方式时，将已经选中的节点进行取消操作。
•Statistics(统计)方式时，单击树结构中的任一节点会弹出一个统计窗口，将所选节点的数学特征进行统计计算并以直方图形式显示出来。
•View Col.Cfs(查看终节点系数值的颜色)方式时，单击树结构中的任一终节点，则在左边图中显式该节点系数值对应的灰度图像。
④Cfs Col.列表框。
染色方式分为FRQ和NAT两类，FRQ表示系数值以频率分布来度量，NAT表示系数值以奈特来度量。
(5)计算最佳树。
在图8-15中的右上方单击Best Tree按钮，结果如图8-16所示。
(6)选择阈值进行去噪。
单击De-noise按钮，弹出如图8-17所示的一维小波包去噪窗口。
图8-16 最佳树
图8-17 去噪窗
再单击View Denoise Signal按钮可以看到，去噪信号与初始信号的对比图，如图8-18所示。可以看出，压缩信号(黑色)和原始信号(红色)在同一个显示框中。
图8-18 去噪信号与初始信号的对比图
8.3 二维小波包在MATLAB中的应用
本节是上一节的续篇，主要针对二维信号(图像)通过利用小波包分析解决实际问题，并且介绍相关MATLAB函数。
8.3.1 二维小波包函数
上一节已经了解到一维小波包分析分解中，细节部分也进行相同的分解。对于二维信号同样也可以利用小波包分析来分析图像的特征。
(1)小波包二维分解函数wpdec2，其调用格式为：
T=wpdec2(X，N，'wname'，E，P)
T=wpdec2(X，N，'wname')
其中：X是信号二维数组；N是分解层数；wname是小波基；E是熵的类型；P是可选参数取决于返回值T；T是返回的小波包树。
(2)小波包的二维重构函数wprec，其调用格式为：
X=wprec(T)
其中：T是小波包树形结构；X是信号矩阵。
【例8.4】利用二维小波包分析函数，将二维车胎图像分解为2层的小波包树，并且选取Haar作为小波基，熵的类型为shannon。
MATLAB代码设置如下：
load tire
image(X);
colormap(map);
如图8-19所示。
t=wpdec2(X，2，'db1');
%画出图像
plot(t)
结果如图8-20所示。
图8-19 车胎图
图8-20 小波包树
【例8.5】利用二维小波包分解函数及能量谱函数，将噪声信号分别分解为第一到第四层，并且分析每层的能量谱。
MATLAB代码设置如下：
close all
%加载信号
load tire;
for i=1:4
wpt=wpdec(s，i，'db1');
e=wenergy(wpt);
E=zeros(1，length(e));
for j=1:2^i
E(j)=sum(abs(wprcoef(wpt，[i，j-1])).^2);
end
% 画出第一到第四层分解图的能量谱
figure(i+1)
bar(e);
axis([0 length(e) 0 130]);
title(['第'，num2str(i)，'层']);
for j=1:length(e)
text(j-0.2，e(j)+20，num2str(e(j)，'%4.2f'));
end
end
结果如图8-21所示。
【例8.6】利用二维小波包分析函数，将噪声信号分解为三层的小波包树，并且选取db3作为小波基，熵的类型为shannon。最后，利用二维小波包重构函数使其重构信号，并且计算误差值。
MATLAB代码设置如下。首先，加载信号：
close all
load tire;
image(X);
colormap(map);
title('原始信号')
axis off
图8-21 第一层到第四层分解的能量谱
结果如图8-22所示。
t=wpdec2(X，3，'db3'，'shannon');
s=wprec2(t);
figure(2)
image(s);
colormap(map);
title('重构信号')
axis off
结果如图8-23所示。
%计算误差
err=X-s;
mesh(err)
结果如图8-24所示。
如图8-22 原始信号
图8-23 重构信号
图8-24 误差分布
8.3.2 二维小波包界面式应用——图像压缩
前一小节主要介绍了一维小波包分析进行信号的压缩与去噪，本小节将通过一个具体的实例来介绍利用二维小波包界面进行图像压缩。
(1)启动小波工具箱。
(2)启动二维小波包的主界面。
单击Wavelet Packet 2-D按钮，弹出二维小波包的主界面。可以看出，二维小波包图形工具窗口左边是信号的显示区域，右边是对信号进行小波包分析的各种按钮和参数选择设置框。在左边的显示区域中有4个图形区，与上一节介绍的功能种类都是一样的，这里就不再介绍了。
(3)下载图像。
选择File→Load Signal命令，在弹出的对话框中选择woman.mat 文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到了此工具中，如图8-25所示。
(4)二维小波包分析。
对信号进行设置。选择 haar 小波，并且分角层数设置为3，选择的熵类型为shannon如图8-25所示。单击Analyze按钮，执行对信号分解，显示在窗口中，如图8-26所示。
图8-25 加载图像
图8-26 分析图像
这里对Entorpy(熵标准)的介绍与上一节是一致的。可选用的熵标准类型有6种。在4 幅图的右下方有个Colored Coefficients for Terminal Nodes(终节点系数染色)图，它显示了哪些高频部分被分解或没被分解。大的正方形表示该部分没有像小的正方形一样进行多层的分解。
(5)计算最佳树。
在图8-26中的右上方单击Best Tree按钮，其功能就是在执行图像的小波包分解之后用熵选择框确定的熵标准计算出最佳的小波包分解树。
(6)进行压缩。
单击Compress按钮，弹出如图8-27所示的二维小波包压缩窗口。
图8-27 二维小波包压缩窗口
再单击Compress按钮可以看到压缩图像与初始图像的对比图，如图8-28所示。
图8-28 前后压缩对比
可以看出，图像保留了97.25%的能量。读者还可以通过试验找出更佳的阈值，使得在图像不失真的前提下压缩比达到最大。
在图8-28中的阈值编辑框中可以输入读者自己定义的阈值，从而获得更好的压缩效果，也可以采用小波变换进行图像压缩与之进行对比，对于同一幅图像，究竟是小波包分析好还是小波压缩好。
最后单击Close按钮，会出现更新合成图像(Update the syethesized image)对话框，单击yes按钮结束压缩过程，回到二维小波包图形工具箱。
8.3.3 二维小波包界面式应用——图像去噪
前一小节主要介绍了利用二维小波包分析进行图像的压缩，本节将通过一个具体的实例来介绍利用二维小波包界面进行图像的去噪。
(1)启动小波工具箱。
(2)启动二维小波包的主界面。单击Wavelet Packet 2-D按钮，弹出二维小波包的主界面。
(3)下载图像。选择File→Load Signal命令，在弹出的对话框中选择detfingr.mat 文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到了此工具中，如图8-29所示。
图8-29 加载图像
(4)对信号进行设置。选择db3小波，并且分角层数设置为3，选择的熵类型为shannon，如图8-29所示。单击Analyze按钮，执行对信号分解，显示在窗口中，如图8-30所示。
图8-30 分析图像
(5)计算最佳树。在图8-30中的右上方单击Best Tree按钮，其功能就是在执行图像的小波包分解之后用熵选择框确定的熵标准计算出最佳的小波包分解树。
(6)进行去噪。单击De-noise按钮，弹出如图8-31所示的二维小波包去噪窗口。
图8-31 二维小波包去噪
再单击De-noise按钮可以看到去噪的图像与初始图像的对比图，如图8-32所示。
图8-32 前后对比
最后单击Close按钮，会出现更新合成图像(Update the syethesized image)对话框，单击yes按钮结束去噪过程，回到二维小波包图形工具箱。
8.4 小波包分析的综合应用实例
本节将通过小波包分析来综合运行其理论，介绍一个信号分析的实例。
【例8.7】通过小波包分析带有随机噪声的信号，分析其障碍信号的功率谱。其噪声信号为：
MATLAB代码设置如下：
% 加载信号
t=1:1000;
s1=10*sin(pi*50*t*0.001)+sin(pi*120*t*0.001)+rand(1，length(t));
for t=1:500;
s2(t)=sin(pi*50*t*0.001)+sin(pi*120*t*0.001)+rand(1，length(t));
end
for t=501:1000;
s2(t)=sin(pi*200*t*0.001)+sin(pi*120*t*0.001)+rand(1，length(t));
end
subplot(9，2，1)
plot(s1)
title('原始信号')
ylabel('S1')
subplot(9，2，2)
plot(s2)
title('故障信号')
ylabel('S2')
wpt=wpdec(s1，3，'db4'，'shannon');
%plot(wpt);
s130=wprcoef(wpt，[3，0]);
s131=wprcoef(wpt，[3，1]);
s132=wprcoef(wpt，[3，2]);
s133=wprcoef(wpt，[3，3]);
s134=wprcoef(wpt，[3，4]);
s135=wprcoef(wpt，[3，5]);
s136=wprcoef(wpt，[3，6]);
s137=wprcoef(wpt，[3，7]);
s10=norm(s130);
s11=norm(s131);
s12=norm(s132);
s13=norm(s133);
s14=norm(s134);
s15=norm(s135);
s16=norm(s136);
s17=norm(s137);
st10=std(s130);
st11=std(s131);
st12=std(s132);
st13=std(s133);
st14=std(s134);
st15=std(s135);
st16=std(s136);
st17=std(s137);
disp('正常信号的特征向量');
snorm1=[s10，s11，s12，s13，s14，s15，s16，s17]
std1=[st10，st11，st12，st13，st14，st15，st16，st17]
subplot(9，2，3);plot(s130);
ylabel('S130');
subplot(9，2，5);plot(s131);
ylabel('S131');
subplot(9，2，7);plot(s132);
ylabel('S132');
subplot(9，2，9);plot(s133);
ylabel('S133');
subplot(9，2，11);plot(s134);
ylabel('S134');
subplot(9，2，13);plot(s135);
ylabel('S135');
subplot(9，2，15);plot(s136);
ylabel('S136');
subplot(9，2，17);plot(s137);
ylabel('S137');
wpt=wpdec(s2，3，'db4'，'shannon');
%plot(wpt);
s230=wprcoef(wpt，[3，0]);
s231=wprcoef(wpt，[3，1]);
s232=wprcoef(wpt，[3，2]);
s233=wprcoef(wpt，[3，3]);
s234=wprcoef(wpt，[3，4]);
s235=wprcoef(wpt，[3，5]);
s236=wprcoef(wpt，[3，6]);
s237=wprcoef(wpt，[3，7]);
s20=norm(s230);
s21=norm(s231);
s22=norm(s232);
s23=norm(s233);
s24=norm(s234);
s25=norm(s235);
s26=norm(s236);
s27=norm(s237);
st20=std(s230);
st21=std(s231);
st22=std(s232);
st23=std(s233);
st24=std(s234);
st25=std(s235);
st26=std(s236);
st27=std(s237);
disp('故障信号的特征向量');
snorm2=[s20，s21，s22，s23，s24，s25，s26，s27]
std2=[st20，st21，st22，st23，st24，st25，st26，st27]
subplot(9，2，4);plot(s230);
ylabel('S230');
subplot(9，2，6);plot(s231);
ylabel('S231');
subplot(9，2，8);plot(s232);
ylabel('S232');
subplot(9，2，10);plot(s233);
ylabel('S233');
subplot(9，2，12);plot(s234);
ylabel('S234');
subplot(9，2，14);plot(s235);
ylabel('S235');
subplot(9，2，16);plot(s236);
ylabel('S236');
subplot(9，2，18);plot(s237);
ylabel('S237');
结果如图8-33所示。
图8-33 原始与故障信号
%fft
figure
y1=fft(s1，1024);
py1=y1.*conj(y1)/1024;
y2=fft(s2，1024);
py2=y2.*conj(y2)/1024;
y130=fft(s130，1024);
py130=y130.*conj(y130)/1024;
y131=fft(s131，1024);
py131=y131.*conj(y131)/1024;
y132=fft(s132，1024);
py132=y132.*conj(y132)/1024;
y131=fft(s133，1024);
py131=y131.*conj(y131)/1024;
y132=fft(s134，1024);
py132=y132.*conj(y132)/1024;
y131=fft(s135，1024);
py131=y131.*conj(y131)/1024;
y132=fft(s136，1024);
py132=y132.*conj(y132)/1024;
y131=fft(s137，1024);
py131=y131.*conj(y131)/1024;
y230=fft(s230，1024);
py230=y230.*conj(y230)/1024;
y231=fft(s231，1024);
py231=y231.*conj(y231)/1024;
y232=fft(s232，1024);
py232=y232.*conj(y232)/1024;
y231=fft(s233，1024);
py231=y231.*conj(y231)/1024;
y232=fft(s234，1024);
py232=y232.*conj(y232)/1024;
y231=fft(s235，1024);
py231=y231.*conj(y231)/1024;
y232=fft(s236，1024);
py232=y232.*conj(y232)/1024;
y231=fft(s237，1024);
py231=y231.*conj(y231)/1024;
f=1000*(0:511)/1024;
subplot(1，2，1);
plot(f，py1(1:512));
ylabel('P1');
title('原始信号的功率谱')
subplot(1，2，2);
plot(f，py2(1:512));
ylabel('P2');
title('故障信号的功率谱')
结果如图8-34所示。
图8-34 原始与故障信号的功率谱
figure
subplot(4，2，1);
plot(f，py130(1:512));
ylabel('P130');
title('S130的功率谱')
subplot(4，2，2);
plot(f，py131(1:512));
ylabel('P131');
title('S131的功率谱')
subplot(4，2，3);
plot(f，py132(1:512));
ylabel('P132');
subplot(4，2，4);
plot(f，py131(1:512));
ylabel('P133');
subplot(4，2，5);
plot(f，py132(1:512));
ylabel('P134');
subplot(4，2，6);
plot(f，py131(1:512));
ylabel('P135');
subplot(4，2，7);
plot(f，py132(1:512));
ylabel('P136');
subplot(4，2，8);
plot(f，py131(1:512));
ylabel('P137');
结果如图8-35所示。
图8-35 S130/S131的功率谱
figure
subplot(4，2，1);
plot(f，py230(1:512));
ylabel('P230');
title('S230的功率谱')
subplot(4，2，2);
plot(f，py231(1:512));
ylabel('P231');
title('S231的功率谱')
subplot(4，2，3);
plot(f，py232(1:512));
ylabel('P232');
subplot(4，2，4);
plot(f，py233(1:512));
ylabel('P233');
subplot(4，2，5);
plot(f，py234(1:512));
ylabel('P234');
subplot(4，2，6);
plot(f，py235(1:512));
ylabel('P235');
subplot(4，2，7);
plot(f，py236(1:512));
ylabel('P236');
subplot(4，2，8);
plot(f，py237(1:512));
ylabel('P237');
结果如图8-36所示。
图8-36 S230/S231的功率谱
8.5 本章小结
本章详细介绍了小波包分析的定义及性质，并且针对相关小波包的分解与重构算法，介绍了MATLAB中有关一维、二维小波包的分解和重构函数(如wpdec和wprec)，还介绍了相应的小波包能量谱函数，同时给出了相应函数的算例。
在最后两节，利用MATLAB中小波工具箱提供的GUI(图形用户界面)，详细地运用这个一维、二维小波包开发工具，进行了信号和图像处理，描述了常用工具按钮及参数设置，让一些对理论并不是非常熟悉的小波初学者可以立即运用此工具，即小波包变换做信号与图像相关的处理。
根据人们认知规律和笔者多年的经验，每一位学习者在学习一门新的学科时都要经过理论与实践上的反复磨练。因此，本章仍然遵循理论与实践相结合的模式来编写，让熟悉小波的用户可以用好小波，让不熟悉的用户可以了解小波的基本知识然后去简单实践。希望学习小波的朋友都能因此而受惠！