第12章 小波分析用于信号识别与检测

时频方法现在发展较快的是利用经验模态分解方法(EMD方法)求解模态参数的HHT方法，另外一种是连续小波模态参数辨识方法。由传感器所检测到的奇异信号往往载有设备运行状态特征的重要信息。判断状态信号的奇异点出现时刻，并对信号奇异性实现定量描述，在信号处理和故障诊断等领域有着重要的意义。
小波变换具有时－频局部化特性，能够有效地分析信号的奇异性，确定奇异点的位置与奇异度的大小，为信号奇异性分析提供了有力的工具。利用连续小波变换方法研究系统模态参数识别是一个新的研究方向。本章将结合具体实例系统介绍小波分析用于信号识别与检测(提取)的方法。
学习目标：
(1)了解信号奇异性检测原理和方法
(2)熟练掌握信号的间断点检测
(3)熟练掌握信号的自相似检测
(4)熟练掌握信号的识别与提取
(5)熟练掌握模态参数识别的应用
(6)熟悉图像的边缘检测
12.1 信号的奇异性检测理论
本节将首先介绍信号奇异性的定义，接下来介绍奇异性检测的基本理论，同时还将介绍信号的奇异性和其与小波变换的关系，以及信号与噪声的小波变换特性。
12.1.1 信号奇异性概念
宇宙射线和太阳黑子爆发，空间核磁暴，对于在太空飞行的卫星和飞船安全构成重大威胁，影响太空飞行器的使用寿命。通过处理采集的空间数据，检测到奇变点，找到太空天气异动的时刻，做出科学的判断，及时调整飞行器姿态，以保护飞行器的安全。寻找变化周期，总结规律，可为进行太空天气预报提供依据。
数学上称无限次可导函数是光滑的或没有奇异性，若函数在某处有间断或某阶导数不连续，则称函数在此处有奇异性，该点就是奇异点。信号的突变点往往包含重要的信息。
一般用Lipschitz指数来刻画信号的奇异性。Lipschitzα越大，函数越光滑。下面给出描述信号奇异性的一般定义。
定义 设有非负整数n(n≤α≤n+1)，如果存在着两个常数C>0和x0
 >0与一个n 阶多项式Pn
 (x)，对于
 使得
则称函数f在点x0
 是Lipsctitzα的。
(1)如果存在一个常数C与一个n阶的多项式Pn
 使得对于任意x0
 <x上式都成立，则称函数在区间(a，b)上是一致Lipschitzrα的。
(2)函数f(x)在点x0
 上的Lipschitz指数为所有满足上式中α的上确界。
(3)如果函数f(x)在x0
 上的Lipschitzr指数α小于1，则称函数在该点是奇异的。
Lipschitz 指数给出了信号 f(x)在 x0
 点可导性的精确信息，某一点连续可微的函数 f(x)在该点上的Lipschitz指数为1。如果某一点上函数f(x)的微分有界，但不连续则f(x)在该点的Lipschitz指数仍然为1。
我们可以证明当f(x)在x0
 上的Lipschitz指数α>1时，则f(x)在x0
 是n次可导的，而多项式Pn
 是f(x)在x0
 上泰勒展开式的前n+1项。
如果n=0，则Pn
 =f(x0
 )。由定义可以得知，当f(x)的Lipschitz指数α满足条件n<α<n+1时，f(x)是n阶可微的，它的第n阶导数是奇异的，即f(x)的n+1阶导数发散。此时如果f(x)的n阶导数有界，则f(x)的n阶导数是正则的。如果f(x)的Lipschitz指数为α，则每积分一次Lipschitz指数增加1，每微分一次Lipschitz指数减少1。
12.1.2 Fourier变换与信号奇异性的关系
Fourier变换是信号奇异性检测的一种传统工具，信号f(x)在R上的一致Lipschitz正则性(即整体正则性)与其Fourier变换的渐近衰减性有关。
定理12-1
(1)如果f满足
则函数f是有界函数且f是p次连续可微的有界变差函数。
(2)如果存在常数K及ε>0使得
则f∈CP
 。此时如果f还具有紧支集，则可以推断 f ∈C∞
 。
由定理12-1可见，信号的衰减性和信号Fourier变换的光滑程度也有一定的联系，如果信号f(x)是n阶连续可导的，则n越大其Fourier变换的衰减也越快。具体地说，如果f(x)具有连续可微的n阶导数且其n阶导数均可积，则f(x)趋于零的衰减速率为x-n
 ，反之却不然。
12.1.3 小波变换与信号的奇异性
定理12-2 (小波变换与一致Lipschitz 指数)设f(x)∈L2 (R)，[a，b]∈R，0<α<1，f(x)在区间[a+ε，b+ε]上是一致 Lipschitzα的。当且仅当存在任意ε>0和常数 A 对任意 s>0及x∈[a+ε，b+ε]有
成立。
上式给出了f(x)的局部而非全局的渐进衰减性。另外，定理12-1只给出了区间上的奇异性，下面的定理给出了在x0
 点上的奇异性度量。
定理12-3 (小波变换和局部Lipschitz 指数)设f∈L2
 (R)，ψ(x)具有n阶消失距且n次连续可微，令
 ，如果f(x)在点x0
 的局部Lipschitz指数为α，则存在常数A使得对任意尺度s在x0
 附近的点x都满足
反之，令α<n且α∉Z，则f(x)在点x0
 的局部Lipschitz指数为α，则如下两个条件成立：
(1)存在ε>0和常数A使得对任意尺度s，x0
 附近的点x满足不等式
(2)存在常数B使得对任意尺度s，x0
 附近的点x满足不等式
上述两个不等式说明，小波变换确实能用来估计函数的局部奇异性，但是很多实际的数值计算中很难直接运用上述两个不等式的结论来检测函数的局部奇异性和估计 Lipschitz指数。
在下一小节我们将引进小波变换的模极大值，利用小波变换的局部极大值和函数奇异性之间的关系同样可以对函数的局部奇异性进行分析，而且具有较小的运算量。
12.1.4 小波变换模极大值点同信号突变点之间的关系
设θ(t)是某一起平滑作用的低通平滑函数且满足条件
 ，为方便起见，不妨取其为高斯函数，即假定θ(t)是二次可导的并且定义ψ(1)
 (t)和ψ(2)
 (t)分别为θ(t)的一阶导数和二阶导数：
则函数ψ(1)
 (t)和ψ(2)
 (t)满足可容许性条件：
因此，ψ(1)
 (t)和ψ(2)
 (t)可作小波母函数。
对于任意函数g(t)引进记号
 表示函数g(t)在尺度因子s下的伸缩。由于小波变换是通过将原信号f(t)与小波进行卷积运算，得到ψ(1)
 (t)为小波的函数，因此将信号f(t)在尺度s下位置τ处的卷积型小波变换定义为：
对应于ψ(2)
 (t)的小波变换为：
由ψ(1)
 (t)和ψ(2)
 (t)之定义可得：
由于 fθs
 可以看成是由低通平滑函数θ(t)在尺度 s 下对函数f(t)进行平滑的结果，由上面式子可见小波变换
 分别是函数f(t)在尺度s下由θ(t)平滑后再取一阶与二阶导数。若取θ(t)为高斯函数函数 f(t)在ψ(1)
 (t)和ψ(2)
 (t)下进行二进卷积型小波变换如图12-1所示。
图12-1 f(t)在θ(t)、ψ(1)
 (t)和ψ(2)
 (t)变换下的结果
由图12-1可知，
 的模局部极值点t0
 、t1
 、t2
 既对
 的过零点，又对应于平滑后信号fθs
 的拐点。当s很小时，用θs
 对f(t)平滑的结果对f(t)的突变部分的位置与形态影响不大；当s较大时，则此平滑过程会将f(t)的一些细小的突变消去，只剩下大尺寸的突变。
由此可知，当小波函数看作是某一平滑函数的一阶导数时，信号的小波变换模的局部极值点对应于信号的突变点；当小波函数看作是某一平滑函数的二阶导数时，信号的小波变换的过零点也对应于信号的突变点。
因此采用检测小波变换系统模的过零点和局部极值点的方法可以检测信号的边缘位置。这是两种相似的方法，但是比较来说，用局部极值点进行检测更具有优越性。
由于函数fθs
 的拐点既对应它的一阶导数模的极大点，又对应于极小值点，而
 的极大值点t0
 、t2
 是对应于信号的快变位置
 的极小值t1
 对应变化最慢的位置。所以，单凭二阶导数的过零点很难区分是信号的突变点还是缓变点。此外，过零点只给出拐点的位置信息，却不能给出变化点的快慢信息。而对局部极大值来说，我们很容易确定该点一定对应信号的快变化点。另外，通过记录模极大点在各个尺度上的取值，我们还可以推测出拐点处的导数值，从而得到信号变化快慢的信息。
鉴于上述原因，用小波变换的模极大值点来测定信号的奇异性与用小波变换的零交叉点相比是有优越性的。
12.1.5 信号与噪声的小波变换特性
假设一维信号s(t)=s0
 (t)+n(t)，n(t)为噪声。s0
 (t)在t的某一邻域内二阶可导，一阶导数S'0
 (t)的极值代表了在极值某一相邻区间内信号线性变化最剧烈的时刻。
我们往往对信号变化最剧烈的时刻，也就是
 最大值所在时刻感兴趣。S'0
 (t)极值点对应于 s0
 (t)的拐点。s0
 (t)的拐点对应于二阶导数S''0
 (t)的零点位置。S''0
 (t)的零点不一定是S'0
 (t)的极值点。
通过S''0
 (t)的零点位置以及零点邻域内值的正负性质，可以判断出 s0
 (t)变化的凹凸特性，也可以判断该点是否是S'0
 (t)在某一邻域内的极值(极大值或极小值)。但是，S'0 (t)的极小值并不一定成为
 的最大值。
观测信号s(t)含有噪声。应用观测数据s(t)代替s0
 (t)进行分析处理时，要进行平滑和滤波处理。设θ(t)是一个适当的光滑函数，满足
 。一般情况下选θ(t)为高斯函数，即
若令
则有
显然θ' (t)和θ2
 (t)均为子波函数，现在对函数θ(t)引入尺度因子a，并采用如下表示：
然后将信号s(t)关于θ' (t)和θ2
 (t)的小波变换定义为：
这里用卷积运算来定义小波变换，与标准的小波变换定义不同，不过没有本质区别。将上式代入可得：
这里用到了卷积和微分可以交换的性质。类似地可得：
从以上分析可知，用观测信号的平滑版本s (t)*θ(t)取代原信号s0
 (t)，W1
 s (t)和
 取代 s0
 (t)的一阶和二阶导数进行分析。表12-1比较了含有噪声的信号和不含噪声突变点和拐点的区别。为了使结果更准确，综合考虑W1
 s(t)和W2
 s(t)多尺度计算的结果来判断信号的突变点。
表12-1 突变点与拐点的检测
12.2 信号的间断点检测
上一节主要从理论的角度介绍信号的奇异性分析理论基础，本节将重点介绍奇异性检测的实例。
12.2.1 第一类间断点检测
当小波函数可看作某一平滑函数的一阶导数时，信号小波变换模的局部极值点对应于信号的突变点；当小波函数可看作某一平滑函数的二阶导数时，信号小波变换的过零点对应于信号的突变点。因此，采用小波变换模的过零点和局部极值点的方法可以检测信号的突变点。比较来说，用局部极值点的方法进行检测更具优越性。
一般信号奇异性分为两种情况：①信号在某一时刻其幅值发生突变，引起信号的不连续，这种类型的突变称为第一种类型的间断点；②信号外观上很光滑，幅值没有发生突变，但是信号的一阶微分有突变发生且一阶不连续，这种类型的突变称为第二种类型的间断点。
应用小波分析可以检测出信号中突变点的位置、类型以及变化的幅度。下面的例子将介绍用小波分析检测第一类间断点。
【例12.1】检测信号幅值变化的准确时间或者是间断点的准确位置。
MATLAB代码设置如下：
load freqbrk;%装载文件名为freqbrk的信号
s=freqbrk;
ls=length(s);
[c，l]=wavedec(s，7，'db4');
plot(s);
title('原始信号');
结果如图12-2所示。
ylabel('s');
a7=wrcoef('a'，c，l，'db4'，7);
figure(2)
plot(a7);
ylabel('a7');
title('用db4小波分解7层的近似系数');
结果如图12-3所示。
图12-2 原始信号
图12-3 近似系数
figure(3)
for i=1:7
decmp=wrcoef('d'，c，l，'db4'，8-i);
subplot(7，1，i);
plot(decmp);
ylabel(['d'，num2str(8-i)]);
end
title('用db4小波分解7层的细节系数');
结果如图12-4所示。
图12-4 细节系数
注：在这个例子中信号的不连续是由于低频特征的正弦信号在后半部分突然加入高频特征的正弦信号的缘故，分析的目的是将加入高频特征的正弦信号的时间检测出来。
从结果波形图12-4上可以看出，信号的不连续点在使用db4小波将信号进行7层分解来检测第一种类型的间断点，可以非常清楚地观察到信号的不连续点，即高频特征的正弦信号的加入点，这是因为间断点包含了高频信息。
原始信号是由两个独立的满足指数方程的信号在 t=500 处连接起来的，因此它看上去是光滑的，但它的一阶微分有突变。
采用 db4 小波对信号分解后，在信号的第一层高频系数 d1中可以明显地看到 t≈500的间断点。要注意的是，在信号奇异点的检测中，选择小波的正则性非常重要，因为这时小波可实现一个长的冲激响应滤波器。
接下来的例子是利用傅立叶分析与小波分析对同一信号进行奇异性检测。上一节已经从理论上论述了两种方法。
【例12.2】用FFT和db4小波实现对信号在频域内的表现形式，并且比较这两种方法的优劣。MATLAB代码设置如下：
load freqbrk;%装载文件名为freqbrk的信号
s=freqbrk;
ls=length(s);
[c，l]=wavedec(s，7，'db4');
plot(s);
title('原始信号');
ylabel('s');
结果如图12-5所示。
%对信号进行FFT变换
fs=fft(s，1000);
fs=abs(fs);
figure;
plot(fs);
ylabel('FFT');
grid;
结果如图12-6所示。
图12-5 原始信号
图12-6 FFT变换
%信号用db4小波分解
[c，l]=wavedec(s，3，'db4');
%对分解结构的第三层低频进行重构
a3=wrcoef('a'，c，l，'db4'，3);
%对分解结构的各高频进行重构
figure
for i=1:3
decmp=wrcoef('d'，c，l，'db4'，4-i);
subplot(3，1，i);
plot(decmp);
ylabel(['d'，num2str(4-i)]);
end
结果如图12-7所示。
图12-7 小波变换
从图12-6中可以看出，由于FFT将信号变换成纯频域中的信号，使得不具有时间分辨能力，所以对信号在时域中的突变点根本无法检测出来。而db4小波分解后的信号则可以很明显地分别出间断点。
对信号进行多尺度分析，在信号出现突变时，其小波变换后的系数具有模量极大值，借此可检测故障发生的时间。
【例12.3】利用小波分解系数检测正弦信号的突变时间起始点。MATLAB代码设置如下：
t=0:pi/100:2*pi;
s1=sin(t);
s2=sin(10*t);
s3=sin(t);
%整个信号
s=[s1，s2，s3];
plot(s)
title('原始信号');
ylabel('s');
结果如图12-8所示。
[c，l]=wavedec(s，7，'sym4');
plot(s);
a7=wrcoef('a'，c，l，'sym4'，7);
figure(2)
plot(a7);
ylabel('a7');
title('用sym4小波分解7层的近似系数');
结果如图12-9所示。
图12-8 原始信号
图12-9 近似系数
figure(3)
for i=1:7
decmp=wrcoef('d'，c，l，'sym4'，8-i);
subplot(7，1，i);
plot(decmp);
ylabel(['d'，num2str(8-i)]);
end
title('用sym4小波分解7层的细节系数');
结果如图12-10所示。
图12-10 细节系数
从图12-10中的小波分解层系数可以明显看出，t=200到t=400时，系统工作出现了异常情况。
12.2.2 第二类间断点检测
前面讨论了信号本身的间断点问题即第一类间断点。在实际应用中，还有一类重要的间断问题，就是导数的间断，即第二类间断点。这种信号的本身可能是连续的，但信号的导数代表了信号的另一重要信息。比如对位移信号一阶导数可能表明了速度信息，二阶导数可能表明了受力信息等。
【例 12.4】对于一给定的线性信号，利用小波分析来检测信号的第二类间断点。MATLAB代码设置如下：
clear
clc
load nearbrk;
s=nearbrk;
plot(s)
title('原始信号');
ylabel('s');
结果如图12-11所示。
[c，l]=wavedec(s，7，'sym4');
plot(s);
a7=wrcoef('a'，c，l，'sym4'，7);
figure(2)
plot(a7);
ylabel('a7');
title('用sym4小波分解7层的近似系数');
结果如图12-12所示。
图12-11 原始信号
图12-12 近似系数
figure(3)
for i=1:7
decmp=wrcoef('d'，c，l，'sym4'，8-i);
subplot(7，1，i);
plot(decmp);
ylabel(['d'，num2str(8-i)]);
end
title('用sym4小波分解7层的细节系数');
结果如图12-13所示。
图12-13 细节系数
此例中，原始信号是一条光滑直线，但是它一阶微分是突变的。利用 sym4 小波对信号进行7层分解后，便将该信号的第二类间断点检测出来(t=500)。
另外，所选择的小波基非常重要，本处使用的是正则性的小波。读者可以选择一个不具有正则性的小波进行分析，看看是否可以检测出来。同时，也可以选择不同的分解层数，若分解为3层，看看是否可以检测出来。
【例 12.5】对一个给定信号，在外观上是光滑的，由两个指数方程接合而成，利用小波分析检测其突变的位置。MATLAB代码设置如下：
t=1:.002:4;
s1=exp(t);
s2=exp(3*t);
s=[s1，s2];
plot(s)
title('原始信号');
ylabel('s');
结果如图12-14所示。
%计算信号一阶导数
ds=diff(s);
figure
plot(ds);
ylabel('信号导数');
结果如图12-15所示。
图12-14 原始信号
图12-15 原始信号的导数
[c，l]=wavedec(s，4，'db4');
figure
a4=wrcoef('a'，c，l，'db4'，4);
plot(a4);
ylabel('a4');
figure(3)
for i=1:4
decmp=wrcoef('d'，c，l，'db4'，5-i);
subplot(4，1，i);
plot(decmp);
ylabel(['d'，num2str(5-i)]);
end
结果如图12-16所示。
图12-16 结果图
从图12-15中可以看出，信号的一阶导数在t=1 500 时发生了突变，而原始信号却是光滑的。但经过小波分解后，从图12-16中的小波分解层系数可以明显看出，t=1 500时，信号出现了异常情况。
12.3 信号的自相似检测
在连续小波变换中，理论上我们可以求得在信号上任意小尺度的小波分解系数，因此连续小波变换可以用来研究信号的自相似或者分形的性质。
在信号处理的意义上，自相似性指的是信号的某一部分在较高的分辨率下采样的结果同信号整体在较低分辨率下的采样结果有一定的相似性。这很符合连续小波变换对任意对尺度进行分解的特点，通过连续小波变换给出的相空间系数，很容易通过系数的分布找到信号的自相似特性。
从小波分析的角度来解释这个问题需要引入相似程度的概念，从直观意义上说，小波分解就是计算一系列信号和小波函数之间的相似系数，相似程度越高则小波系数的绝对值越大，反之则越小。
【例12.6】给定信号，利用小波分析检测信号的自相似性。MATLAB代码设置如下：
load scddvbrk;
s=scddvbrk;
plot(s)
title('原始信号');
ylabel('s');
结果如图12-17所示。
figure
f=cwt(s，[1:2:128]，'db4'，'plot');
title('小波分解自相似指数')
xlabel('时间')
ylabel('尺度变换')
结果如图12-18所示。
图12-17 原始信号
图12-18 小波分解系数图
目前，许多人正在做这项工作，其结果表明，采用小波分解可以很好地研究信号或图像的分形特征。当分形的特征开始时随着时间的发展而变换，然后又不变时，则这种信号被称为多分形。小波分析工具非常适合于分形的实际研究和分形的生成。
12.4 信号识别与信号提取
一些受噪声影响的信号，其发展趋势很难分辨，还有一类信号，经常需要提取其特征来识别其突变程度，本节将列举这方面的实例。
12.4.1 信号发展趋势的识别
由于噪声的干扰，对于实用信号的发展趋势在时域中难以看出，但是，通过小波分解后可以去除那些干扰信号。
【例12.7】利用小波分析以下信号的发展趋势。MATLAB代码设置如下：
t=0:pi/200:1;
s1=sin(t);
l=length(s1);
s2=randn(1，l);
s=s1+s2;
plot(s);
title('原始信号');
ylabel('s');
结果如图12-19所示。
图12-19 原始信号
[c，l]=wavedec(s，7，'db4');
figure
for i=1:7
decmp=wrcoef('a'，c，l，'db4'，8-i);
subplot(7，1，i);
plot(decmp);
ylabel(['a'，num2str(8-i)]);
end
结果如图12-20所示。
图12-20 近似分解
从原始信号中可以看出，由于噪声的污染，信号的发展趋势是不可见的。而在此例子中利用db4小波分解到第7层，则信号的发展趋势随着近似变换，而且变得越来越清晰。信号中的低频部分(a7)代表着信号的发展趋势。在小波分析中，则对应着最大尺度小波变换的低频系数。因而，随着尺度的增加，时间分辨率的降低，信号的发展趋势会表现得更加明显。此外，还可以在频率中理解它的含义，即尺度分解中的低频部分随着层次的增加，其含的高频信息会随之减小。当分解到下一层时，就有更高一些频率信息被滤掉，而剩下的就是信号的发展趋势。
在此可以看出，在展示信号的发展趋势时，小波分析是有用的。这种分析所具有的另一个目的是将隐藏在噪声或其他高频信号中的信号显示出来。这里要强调的是，所有识别的信号本身不具有大的突变。这是因为信号的发展趋势是由信号的低频部分所表征的，如果在信号本身中包含有很大突变，那么在多尺度小波变换的低频部分中，显示出来的信号会和原始信号有很大差别，因为这种变换将信号本身的突变当作高频信息给滤掉了。
12.4.2 某一频率区间上信号的识别
在傅立叶分析中，如果一个信号是由几个不同频率的正弦信号组成，则傅立叶变换可以很有效地分辨这些不同频率的正弦信号。
【例12.8】利用小波分析将给定信号的各频率分开。MATLAB代码设置如下：
t=0:pi/100:4*pi;
s1=sin(t);
s2=sin(5*t);
s3=sin(50*t);
s=s1+s2+s3;
plot(s);
title('原始信号');
ylabel('s');
结果如图12-21所示。
figure;
[c，l]=wavedec(s，6，'db4');
%对分解结构中各低频部分进行重构
for i=1:6
decmp=wrcoef('a'，c，l，'db4'，7-i);
subplot(6，2，2*i-1);
plot(decmp);
ylabel(['a'，num2str(7-i)]);
end
%对分解结构中各高频部分进行重构
for i=1:6
decmp=wrcoef('d'，c，l，'db4'，7-i);
subplot(6，2，2*i);
plot(decmp);
ylabel(['d'，num2str(7-i)]);
end
结果如图12-22所示。
图12-21 原始信号
图12-22 分解系数
%画出d1放大波形图
figure
d1=wrcoef('d'，c，l，'db4'，1);
plot(d1(1:100));
结果如图12-23所示。
在显示的分解结果中，可以看到低频第五层将正弦信号中的最低频率组成清晰地分离出来。该信号是由3个不同频率的正弦信号所组成，其频率在小波分解下的相对频率分别为：高、中、低。由信号频率的组成部分和小波分解下各层频率的分布可知，高频正弦定位于d1层。
图12-23 d1放大
总之，我们能够用小波分析将合成信号中的单纯正弦信号的频率提取出来。因为在小波分解下，不同的尺度具有不同的时间和频率分辨率，因而小波分解能够将信号的不同频率成分分开。
12.4.3 信号的特征提取
如何提取信号特征即检测信号，在信号处理时是一个关键的技术难点。信号的突变点通常是其重要特征。此外，信号的频率谱和它的幅值也表征了信号的许多信息。因此，在进行信号特征提取的研究中，信号的连续性分析、频率谱分析和幅值分析是不可或缺的。
应用小波理论进行信号特征提取时主要有两种方法，即处理边界和滤波。对于边界处理使用了延拓法。从理论上来说，对于有限区间上的数据应构造平方连续可积空间的尺度函数与小波，从而得到此空间上的小波基，这时的小波基往往带有边界条件。
在信号分析中，当对信号进行采样后，就得到在一个大的有限频带中的信号，对这个信号进行小波分解，其实质就是把采集的信号分成两个部分，即高频与低频，而低频部分通常包含了信号的主要信息。对于小波变换而言，高频部分不再分解。若想分解，则需要用小波包分析方法。这样下去，就可以提取信号特征，完成检波过程。
下面的实例介绍有关小波分析在检测方面的优越性。
【例12.9】利用离散小波变换检测信号的突变点。MATLAB代码设置如下：
load nearbrk;
s=nearbrk;
plot(s);
title('原始信号');
ylabel('s');
结果如图12-24所示。
figure;
[c，l]=wavedec(s，4，'db4');
%对分解结构中各低频部分进行重构
for i=1:4
decmp=wrcoef('a'，c，l，'db4'，5-i);
subplot(4，2，2*i-1);
plot(decmp);
ylabel(['a'，num2str(5-i)]);
end
%对分解结构中各高频部分进行重构
for i=1:4
decmp=wrcoef('d'，c，l，'db4'，5-i);
subplot(4，2，2*i);
plot(decmp);
ylabel(['d'，num2str(5-i)]);
end
结果如图12-25所示。
图12-24 原始信号
图12-25 系数分解
从图12-25中可以清晰看到，对信号小波分解将原信号的近似特征和细节特征提取出来。在原始信号上，无法得知原信号导数的不连续性。
由于小波包分析可以分解高频信息，因此下面的实例将利用小波包分解。
【例12.10】利用小波包检测信号的突变点。MATLAB代码设置如下：
load nearbrk;
s=nearbrk;
plot(s);
title('原始信号');
ylabel('s');
结果如图12-26所示。
%小波包分解
t=wpdec(s，2，'db4');
plot(t);
figure;
结果如图12-27所示。
%对第二层分解结构中各系数进行重构
for i=1:4
decmp=wprcoef(t，[2，i-1]);
subplot(4，1，i);
plot(decmp);
ylabel(['a'，num2str(7-i)]);
end
结果如图12-28所示。
图12-26 原始信号
图12-27 小波树
图12-28 结果图
从结果表明，这种基于小波包分析得到的图形可以作为故障特征进行故障诊断。这种诊断方法无需被诊断系统的数学模型，就可以迅速地进行故障检测，并能准确地进行故障定位。
12.5 模态参数识别介绍
模态参数识方法是近年来新出现的。该方法是将信号在时频平面上展开求解特征参数，而不是单独地在时域或频域进行处理。
12.5.1 模态分析的时频辨识方法概述
模态分析的时频方法最早出现的是Hilbert变换方法，由于这种方法存在较大缺陷，只存在了很短的时间就被其他时频方法所取代。
时频方法现在发展较快的是利用经验模态分解方法(EMD方法)求解模态参数的HHT方法，另外一种是连续小波模态参数识别方法。
经验模态分析方法最早是用于非平稳信号处理，信号经 EMD 分解后各分量 IMF (Intrinsic Mode Function)都是平稳的，可以进一步进行Hilbert 变换得到Hilbert 谱，由此得到的Hilbert谱能准确地反应出物理过程中能量在各种频率尺度及时间上的分布。
EMD 方法为非平稳信号进行变换奠定了基础，美国宇航中心 NASA 将其称为HHT (Hilbert-Huang Transform)方法。
HHT方法本质上是对一个信号进行平稳化处理分解，产生具有信号的不同特征尺度的本征模函数IMF分量。对于非平稳信号，对其直接进行Hilbert变换本身没有意义。
谱是一个三维(时间—频率—局部振幅)谱形。在线性框架下，HHT谱与小波谱有相同的表现特性，与FFT谱相比，从Hilbert谱中不仅可以看出幅值，而且可以看出频率随时间的变化情况，这是FFT方法所不具有的优点。
HHT方法能客观地处理一类非线性问题，所得到的三维谱形能准确地用于波内调制机制反映出的系统非线性变化特征，这也是其优点之一。但是，利用HHT方法进行模态参数识别时，在系统内阻尼比较大的情况下，会由于Hilbert方法本身的处理不力而出现较大的误差，而且该方法存在边缘效应问题。
系统模态参数识别的连续小波变换方法研究在国外是从1997年开始发展起来的一个模态分析领域的新方向。该方法从本质上讲，是利用了小波函数的衰减特性，该特性与机械系统自由响应信号的衰减在形式上有相似的地方。根据平稳相位理论，对自由响应信号进行连续小波变换可以在时间—尺度(时间—频率)平面上形成一条清晰的脊线。提取该脊线上的小波变换系数，就可以分离出系统的模态参数。该方法最早见于1997年Staszewski发表的文章。其最先采用的方法是离散小波，但提出了利用连续小波进行参数识别的理论，并提出了3种辨识方法。
随后采用该方法进行了多自由度系统的自由响应参数辨识，并指出了该方法在瞬时频率辨识上优于Hilbert变换。对于线性系统的阻尼、固有频率、振型的辨识都有学者进行了理论和实验研究。
对于非线性系统的参数辨识，Staszewski 采用了瞬时相位的概念识别了系统的瞬时频率和幅值变化。对于非线性系统的研究仍处在起步阶段，未见到更多学者的研究情况发表，这也是这一方法未来应用的主要方面。
对于小波变换母小波的选择，最开始采用的是Morlet小波，有学者采用Cauchy小波和Harmonic小波进行数学运算。
Simonovski对多种母小波进行了比较研究，指出利用Morlet小波进行模态参数辨识的优势。
12.5.2 信号的小波脊提取及计算方法
在小波分析的研究中已经发现了许多合适的小波函数，如Haar小波、样条函数等，它们都有自己适合的应用领域。由测不准定理知：不可能求得一个窗函数具有小于或等于Gaussian函数作为窗函数时窗的面积。
信号在连续小波平面上的分布会呈现出类似地形图中山脊的形状，称为小波脊线。理论分析表明，沿着小波变换时频平面中脊线上分布的参数与原始信号之间有着很强的相似性，能够用来描述原始信号的重要参数。
脊线上分布数据的起伏变化、脊线的位置都有实际的物理意义，直接对应着信号的幅值与频率的变化。
对时频平面上的小波脊线，可以从小波变换的幅值或相位中进行提取。理论上讲，以小波变换的相位变化提取脊线更为精确，但是由于实际应用中该方法涉及相位求导，数值计算的结果可能导致结果并不十分理想。
以小波变换的幅值或相位信息为参考进行脊线提取的方法，已有多位学者进行了研究，下面首先介绍典型的提取方法。
根据相位信息提取小波脊线的一种较为精确的方法称为Marseille方法。该方法的缺点主要是对噪声非常敏感，只有在比较高的信噪比的情况下才能准确地提取信号小波变换的脊线。该方法的另一个缺点是只适用于信号中含有单一成分的情况，当信号含有多个成分的情况时，即使各个成分在小波变换平面上能够清楚地分辨出来，也不能用该方法同时进行多条脊线的提取。只能通过设定时频平面一定的提取范围，使该范围内只含有单一成分的脊线的办法进行改进。显然，这种改进也是不符合实际需要的。
根据模值信息提取小波脊线该方法又称为“Carmona”法，是由Carmona等人于1997年提出的，并于1998年对该方法进行了详细的阐述。应用该方法进行小波脊线提取时，精确性相对较差(相对于Marseille 方法)，但是对于含有噪声的信号有更好的鲁棒性。该方法一个显著的特点是它不需要提供任何噪声信息。然而它允许结合脊线的先验信息，如脊线的平滑性，因此可以认为它是一种贝叶斯方法。
12.5.3 基于小波包和改进HHT的瞬时特征分析
希尔伯特—黄变换(Hilbert-Huang Transfrom，HHT)时频分析方法是Norden E.Huang 等人于1998年提出的一种新的信号分析方法，特别适于非线性、非平稳信号的时频分析。
HHT的主要目的是研究信号的瞬时频率，它的最大特色是通过对信号的经验模态分解(Empirical Mode Decomposition，EMD)使非平稳信号平稳化，从而使瞬时频率变得有意义，进而导出有意义的希尔伯特变换时频谱，这种时频分析方法无需采用信号的先验知识，基函数本身就是自适应地从原信号中分解而得，是后验、自适应。
有的文献提出用小波去噪作为预处理，再运用 EMD和希尔伯特变换方法提取信号瞬时特征。相应的文献提出对本征模态函数(Intrinsic Mode Function，IMF)分量进行小波分析的一种信号瞬时特征分析方法。还有人提出了一种EMD与WPT 相结合的新方法，该方法首先对信号进行 EMD 分解，然后对出现混成模态的固有模态函数进行小波包分解，从而消除模态混叠的影响。
以上方法各有各的优缺点。HHT像其他信号分析方法一样，也存在一些不足之处：通过EMD 得到的第一个高频IMF分量可能覆盖着很宽的频率范围，并不是单分量信号；在低频区域可能产生虚假 IMF 分量。尽管如此，HHT 也是在一定理论基础上提出来的，只是其理论依据尚不完备，有些方法还仅仅是依据实验提出来的。
因此，既使 HHT 已在许多工程应用实践中证明能够获得奇妙的效果，但其理论依据的进一步研究和理论框架的进一步完善是不可回避的问题。
12.5.4 模态参数识别的应用
上面几节是针对时频分析方法的综述，本节将通过几个具体算例来了解这些方法的应用。
【例12.11】为了验证EMD方法处理高频数据，利用MATLAB对原始信号进行EMD分解，可实现数据的各个IMF分量和瞬时频率，并能对Hilbert时谱进行刻画。
MATLAB代码设置如下。
首先，由5个正弦信号合成为一个信号：
clear
clc
f1=400;f21=310;f22=100;f3=200;f4=250;f5=150;
Fs=1000;
for i=1:Fs+1
t(i)=(i-1)/Fs;
if 0.25<t(i) & t(i)<=0.5
Ut1=1;
else Ut1=0;end
if 0.5<t(i) & t(i)<=0.75
Ut2=1;
else Ut2=0;end
if 0.75<t(i) & t(i)<=1
Ut3=1;
else Ut3=0;end
if 0<t(i) & t(i)<=0.25
Ut4=1;
else Ut4=0;end
x1(i)=sin(2*pi*f1*t(i));
x2(i)=(sin(2*pi*f21*t(i)))*Ut1;
x3(i)=4*sin(2*pi*f3*t(i))*Ut2;
x4(i)=sin(2*pi*f4*t(i))*Ut3;
x5(i)=sin(2*pi*f5*t(i))*Ut4;
end
subplot(6，1，1);
plot(t，x1);xlabel('Time')，ylabel('x1');
subplot(6，1，2)
plot(t，x2);xlabel('Time')，ylabel('x2');
subplot(6，1，3)
plot(t，x3);xlabel('Time')，ylabel('x3');
subplot(6，1，4)
plot(t，x4);xlabel('Time')，ylabel('x4');
subplot(6，1，5)
plot(t，x5);xlabel('Time')，ylabel('x5');