CN105426822A - 基于双树复小波变换的非平稳信号多重分形特征提取方法 - Google Patents

Abstract

本发明公开了一种基于双树复小波变换的非平稳信号多重分形特征提取方法，步骤为：对待分析的非平稳信号进行集成处理；对集成信号进行双树复小波变换，利用小波分解尺度系数和细节系数得到各尺度下信号的波动成分；利用得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小；结合尺度值大小，对各尺度下的波动成分进行分段；计算信号的各阶波动函数，利用波动函数和尺度值的双对数关系，经过最小二乘拟合得到广义<i>hurst</i>指数，得到各阶尺度指数；利用legendre变换得到信号的多分形奇异谱。本发明利用双树复小波变换进行信号分解，克服了传统小波变换缺乏平移不变性，保证了多重分形特征提取的准确性，运算速度快，有利于在线应用。

Description

基于双树复小波变换的非平稳信号多重分形特征提取方法

技术领域

本发明涉及非平稳信号处理方法的技术领域，具体涉及一种基于双树复小波变换的非平稳信号多重分形特征提取方法。

背景技术

在机械设备状态监控、紊流分析、心电、脑电等信号分析领域，所处理的对象多为非平稳信号，对信号进行分析的一个关键步骤是提取信号的特征，分形特征是其中的重要一类。相比单分形分析，多重分形分析更适合非平稳信号，它能够实现对信号局部尺度行为更精细的刻画，从而为进一步的分析提供更丰富的信息。

传统的多重分析工具是盒子法，但是在求信号波动函数时，其时间尺度是人为确定的，对信号没有自适应性。基于经验模式分解的方法是一种自适应的分析方法，但是其计算复杂，并且经验模式分解过程往往出现一些无关的模式。小波分析技术是信号处理领域中强有力的处理工具，与傅里叶分析技术相比，它在时域和频域具有良好的局部化性质及多分辨率特性，因而受到广泛的应用。基于连续小波变换的小波模极大方法是多重分形分析的经典方法，但是连续小波变换计算复杂。常规离散小波采用金字塔算法，运算速度快，但是其不满足平移不变性，直接用于多重分形分析会产生较大误差。

发明内容

为了解决上述技术问题，本发明提供了一种基于双树复小波变换的非平稳信号多重分形特征提取方法，采用小波变换的方法来确定自适应的信号的趋势和时间尺度，利用双树复小波变换，克服了传统小波不满足平移不变性的缺陷，提高了自适应性和运算速度，更准确、快捷的提取非平稳信号的多重分形特征。

为了达到上述目的，本发明的技术方案是：一种基于双树复小波变换的非平稳信号多重分形特征提取方法，其步骤如下：

步骤1：对待分析的非平稳信号进行集成处理以突出信号的分形特性；

步骤2：对集成信号进行双树复小波变换，利用小波分解尺度系数和细节系数得到各尺度下信号的波动成分；

步骤3：利用得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小；

步骤4：结合尺度值大小，对各尺度下的波动成分进行分段；计算信号的各阶波动函数，利用波动函数和尺度值的双对数关系，经过最小二乘拟合得到广义hurst指数，进而得到各阶尺度指数；利用legendre变换，得到信号的多分形奇异谱。

所述对信号进行集成处理的方法是：

其中，x(k)是原始信号，k＝1,...,t；＜x＞为信号的均值；N为信号的数据点数。

所述得到各尺度下信号的波动成分方法的步骤为：

a.选择双树复小波滤波器；ψ_h(t),ψ_g(t)分别为双树复小波变换采用的实值小波函数，φ_h(t),φ_g(t)分别为对应的尺度函数，小波函数与尺度函数互为希尔伯特变换对；

b.利用双树复小波滤波器对信号进行M层分解，分别得到小波系数和尺度系数其中1≤l≤M；构成信号在1≤l≤M尺度下的复小波系数

c.对第l尺度的小波系数进行单支重构，得到重构信号d_l(t)，其中1≤l≤M；对第M尺度的尺度系数的进行单支重构，得到重构信号c_M(t)；

d.每个尺度l的趋势表示为每个尺度l的波动成分表示为Fl_l(t)＝y(t)-Tr_l(t)。

所述利用得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小的方法的步骤为：

i.对第l尺度小波系数重构信号d_l(t)进行希尔伯特变换，得到该尺度下的解析信号：其中，为d_l(t)的希尔伯特变换；

ii.由解析信号z_l(t)得到信号的相位角利用相位角的微分得到瞬时频率ω_l(k),k＝1,...,N/2^l；

iii.利用瞬时频率ω_l(k)得到对应尺度l的时间尺度的大小s_l＝1/＜ω_l(k)＞，其中，＜ω_l(k)＞为ω_l(k)在该尺度的均值。

所述得到广义hurst指数、各阶尺度指数、多分形奇异谱的方法是：

①对第l尺度沿信号的正反两个方向，利用时间尺度s_l对波动成分进行无覆盖的分段，共得到2N_s段，每段记为ε_v(i)，i＝1,...,s_l；

②计算每段的局部波动函数

$F^2(v,s_l)=\frac{1}{s_l}\underset{i=1}{\overset{s_l}{\&Sigma;}}{\&lsqb;{\mathrm{\&epsiv;}}_v\left(i\right)\&rsqb;}^2;$

③取q∈[-q_lim,+q_lim]，计算q＝0之外的各阶波动函数

$F_q\left(s_l\right)={\left\{\frac{1}{2N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}{\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;}^{q/2}\right\}}^{1/q};$

对q＝0，波动函数为

$F_0\left(s_l\right)=\exp \{\frac{1}{4N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}ln\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;\}$

④对logF_q(s_l)和logs_l进行最小二乘拟合，所得斜率即为广义hurst指数h(q)，计算得到尺度指数：τ(q)＝qh(q)-1；利用legendre变换，计算得到信号的奇异指数α和多分形奇异谱f(α)分别为α＝h(q)+qh'(q)、f(α)＝qα-τ(q)，其中h'(q)是广义hurst指数h(q)的legendre变换函数。

本发明利用双树复小波变换对信号进行M层分解，得到信号各尺度下的复数形式小波系数和尺度系数，利用第i+1～M层小波系数和第M层尺度系数的单支重构信号叠加，得到第i尺度对应的趋势项，进而得到该尺度的波动项；同时，利用第i尺度的小波系数，经过希尔伯特变换得到对应尺度时间尺度值的估计；利用时间尺度值对该尺度的信号波动进行分段，取不同的阶q计算信号在该尺度下的波动函数；波动函数与q的对数最小二乘拟合斜率对应信号的广义hurst指数，进一步可以得到信号的尺度指数。通过legendre变换，可以得到信号的多分形奇异谱。本发明充分利用了小波变换的信号自适应性，将信号分解为物理含义清楚的多尺度表示形式；利用双树复小波变换进行信号分解，克服了传统小波变换缺乏的平移不变性，保证了多重分形特征提取的准确性；双树复小波变换属于离散小波变换方法，具有快速算法，运算速度快，有利于在线应用。

附图说明

图1为本发明的流程图。

图2为p-模型乘性级联信号图。

图3为p-模型乘性级联集成信号图。

图4为双树复小波分解与重构的示意图。

图5为仿真信号8个尺度趋势图。

图6为仿真信号8个尺度波动图。

图7为仿真信号广义hurst指数特征。

图8为仿真信号广义尺度指数特征。

图9为仿真信号多重分形谱特征。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

一种基于双树复小波变换的非平稳信号多重分形特征提取方法，其步骤如下：

步骤1：对待分析的非平稳信号进行集成处理以突出信号的分形特性。

对信号进行集成处理的方法是：其中，x(k)是原始的非平稳信号，k＝1,...,t；＜x＞为信号x(k)的均值；N为信号的数据点数。对非平稳信号x(k)集成处理可以突出信号的分形特性。

步骤2：对集成信号进行双树复小波变换，利用小波分解尺度系数和细节系数得到各尺度下信号的波动成分。

得到各尺度下信号的波动成分方法的步骤为：

a.选择双树复小波滤波器；ψ_h(t),ψ_g(t)分别为双树复小波变换采用的实值小波函数，φ_h(t),φ_g(t)分别为对应的尺度函数，小波函数与尺度函数互为希尔伯特变换对。

b.利用双树复小波滤波器对集成信号y(t)进行M层分解，分别得到小波系数和尺度系数其中1≤l≤M。利用小波系数和尺度系数构成信号在1≤l≤M尺度下的复小波系数

c.对第l尺度的小波系数进行单支重构，得到重构信号d_l(t)，其中1≤l≤M；对第M尺度的尺度系数的进行单支重构，得到重构信号c_M(t)；

d.每个尺度l的趋势表示为每个尺度l的波动成分表示为Fl_l(t)＝y(t)-Tr_l(t)。

步骤3：利用得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小。

利用步骤2得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小的方法的步骤为：

i.对第l尺度小波系数重构信号d_l(t)进行希尔伯特变换，得到该尺度下的解析信号：其中，为d_l(t)的希尔伯特变换。

ii.由解析信号z_l(t)得到信号的相位角利用相位角的微分得到瞬时频率ω_l(k),k＝1,...,N/2^l。

iii.利用瞬时频率ω_l(k)得到对应尺度l的时间尺度的大小s_l＝1/＜ω_l(k)＞，其中，＜ω_l(k)＞为ω_l(k)在该尺度的均值。

步骤4：结合尺度值大小，对各尺度下的波动成分进行分段；计算信号的各阶波动函数，利用波动函数和尺度值的双对数关系，经过最小二乘拟合得到广义hurst指数，进而得到各阶尺度指数；利用legendre变换，得到信号的多分形奇异谱。

广义hurst指数、尺度指数、多分形奇异谱为信号的多重分形特征。得到广义hurst指数、各阶尺度指数、多分形奇异谱的方法是：

①对第l尺度沿信号的正反两个方向，利用时间尺度s_l对波动成分进行无覆盖的分段，共得到2N_s段，每段记为ε_v(i)，i＝1,...,s_l。

②计算每段的局部波动函数：

③取q∈[-q_lim,+q_lim]，计算q＝0之外的各阶波动函数

$F_q\left(s_l\right)={\left\{\frac{1}{2N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}{\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;}^{q/2}\right\}}^{1/q};$

对q＝0，波动函数为： $F_0\left(s_l\right)=\exp \{\frac{1}{4N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}ln\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;\}.$

④对各阶波动函数F_q(s_l)和时间尺度s_l取对数，对logF_q(s_l)和logs_l进行最小二乘拟合：logF_q(s_l)＝h*logs_l所得斜率h即为广义hurst指数h(q)，计算得到尺度指数：τ(q)＝qh(q)-1；利用legendre变换，计算得到信号的奇异指数α和多分形奇异谱f(α)分别为α＝h(q)+qh'(q)、f(α)＝qα-τ(q)，h'(q)是广义hurst指数h(q)的legendre变换函数。

具体实例：

一种基于双树复小波变换的非平稳信号多重分形特征提取方法，对典型的非平稳信号p-模型乘性级联信号，模型参数p₁＝0.3和p₂＝0.7，信号长度为2¹⁶，如图1所示，用如下步骤进行处理：

步骤1：对信号进行集成处理，如下式所示：

其中，x(k)是原始信号,k＝1,...,t；＜x＞为信号的均值。原始信号与集成信号分别如图2和图3所示。

步骤2：选择双树复小波滤波器，第一层两个树分解均采用(13，19)阶近似对称的双正交滤波器，滤波器系数为：

其余各层分析选用14阶线性相位Q平移滤波器，滤波器系数分别为:

对图3中的集成信号进行8层双数复小波分解，得到小波系数和尺度系数其中1≤l≤8，构成信号在1≤l≤8尺度下的复小波系数

对1≤l≤8尺度进行小波系数单支重构，得到信号细节成分d_l(t)；对第8尺度进行尺度系数的单支重构，得到信号近似成分c₈(t)。

对1≤l≤8尺度，计算每个尺度的趋势成分：如图5所示；对应的波动成分为Fl_l(t)＝y(t)-Tr_l(t)，如图6所示。

步骤3：对1≤l≤8尺度的信号细节成分d_l(t)进行希尔伯特变换，得到对应尺度下的解析信号其中，

那么瞬时频率

对1≤l≤8尺度，计算每个尺度对应的尺度大小其中，

$<{\mathrm{\&omega;}}_l\left(k\right)>=\frac{1}{K}\underset{k=1}{\overset{K}{\&Sigma;}}{\mathrm{\&omega;}}_l\left(k\right).$

步骤4：对1≤l≤8的每个尺度，沿波动成分的正反两个方向利用s_l对其进行无覆盖的分段，共得到2N_s段，每段记为ε_v(i)，i＝1,...,s_l。

利用计算每段的波动函数。在区间[-10,10]内取101个值作为阶数q的取值，利用下式计算各阶波动函数

$\left\{\begin{array}{cc}{F}_q\left(s_l\right)={\left\{\frac{1}{2N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}{\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;}^{q/2}\right\}}^{1/q}& q\&NotEqual;0\\ F_0\left(s_l\right)=\exp \{\frac{1}{4N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}ln\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;\}& q=0\end{array}\right..$

对logF_q(s_l)和logs_l进行最小二乘法拟合，直线的斜率即为广义hurst指数h(q)，如图7所示为所得结果与理论值的对比。由τ(q)＝qh(q)-1得尺度指数τ(q)，如图8所示为所得结果与理论值的对比。利用α＝h(q)+qh'(q)、f(α)＝qα-τ(q)得信号的奇异指数α和多分形奇异谱f(α)，如图9所示为所得结果与理论值的对比。从图中可以看出利用基于双树复小波变换的多重分形方法得到的多重分形特征，包括广义hurst指数、尺度指数、多分形奇异谱，与理论值非常接近，在特征曲线的大部分区域，二者几乎是重合的。

本发明利用双树复小波分解能够对信号进行自适应分解，具有信号分解的平移不变性，保证了多重分形特征提取的准确性；双树复小波分解利用金字塔快速算法，比传统的平稳小波变换方法和连续小波变换方法效率更高。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (5)

1.一种基于双树复小波变换的非平稳信号多重分形特征提取方法，其特征在于，其步骤如下：

步骤1：对待分析的非平稳信号进行集成处理以突出信号的分形特性；

步骤2：对集成信号进行双树复小波变换，利用小波分解尺度系数和细节系数得到各尺度下信号的波动成分；

步骤3：利用得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小；

步骤4：结合尺度值大小，对各尺度下的波动成分进行分段；计算信号的各阶波动函数，利用波动函数和尺度值的双对数关系，经过最小二乘拟合得到广义hurst指数，进而得到各阶尺度指数；利用legendre变换，得到信号的多分形奇异谱。

2.根据权利要求1所述的基于双树复小波变换的非平稳信号多重分形特征提取方法，其特征在于，所述对信号进行集成处理的方法是： $y\left(t\right)=\underset{k=1}{\overset{t}{\&Sigma;}}\&lsqb;x\left(k\right)-<x>\&rsqb;,t=1,...,N;$

其中，x(k)是原始信号，k＝1,...,t；<x>为信号的均值；N为信号的数据点数。

3.根据权利要求1所述的基于双树复小波变换的非平稳信号多重分形特征提取方法，其特征在于，所述得到各尺度下信号的波动成分方法的步骤为：

a.选择双树复小波滤波器；ψ_h(t),ψ_g(t)分别为双树复小波变换采用的实值小波函数，φ_h(t),φ_g(t)分别为对应的尺度函数，小波函数与尺度函数互为希尔伯特变换对；

b.利用双树复小波滤波器对信号进行M层分解，分别得到小波系数和尺度系数其中1≤l≤M；构成信号在1≤l≤M尺度下的复小波系数

c.对第l尺度的小波系数进行单支重构，得到重构信号d_l(t)，其中1≤l≤M；对第M尺度的尺度系数的进行单支重构，得到重构信号c_M(t)；

d.每个尺度l的趋势表示为每个尺度l的波动成分表示为Fl_l(t)＝y(t)-Tr_l(t)。

4.根据权利要求1所述的基于双树复小波变换的非平稳信号多重分形特征提取方法，其特征在于，所述利用得到的各尺度小波系数，估算各尺度的瞬时频率，得到各尺度的时间尺度值的大小的方法的步骤为：

i.对第l尺度小波系数重构信号d_l(t)进行希尔伯特变换，得到该尺度下的解析信号：其中，的希尔伯特变换；

ii.由解析信号z_l(t)得到信号的相位角利用相位角的微分得到瞬时频率ω_l(k),k＝1,...,N/2^l；

iii.利用瞬时频率w_l(k)得到对应尺度l的时间尺度的大小s_l＝1/<w_l(k)>，其中，<ω_l(k)>为ω_l(k)在该尺度的均值。

5.根据权利要求1所述的基于双树复小波变换的非平稳信号多重分形特征提取方法，其特征在于，所述得到广义hurst指数、各阶尺度指数、多分形奇异谱的方法是：

①对第l尺度沿信号的正反两个方向，利用时间尺度s_l对波动成分进行无覆盖的分段，共得到2N_s段，每段记为ε_v(i)，i＝1,...,s_l；

②计算每段的局部波动函数

$F^2(v,s_l)=\frac{1}{s_l}\underset{i=1}{\overset{s_l}{\&Sigma;}}{\&lsqb;{\mathrm{\&epsiv;}}_v\left(i\right)\&rsqb;}^2;$

③取q∈[-q_lim,+q_lim]，计算q＝0之外的各阶波动函数

$F_q\left(s_l\right)={\left\{\frac{1}{2N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N_s}{\&Sigma;}}{\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;}^{q/2}\right\}}^{1/q};$

对q＝0，波动函数为

$F_0\left(s_l\right)=\exp \{\frac{1}{4N_s}\underset{\mathrm{\&nu;}=1}{\overset{2N}{\&Sigma;}}ln\&lsqb;F^2(\mathrm{\&nu;},s_l)\&rsqb;\}$

④对logF_q(s_l)和logs_l进行最小二乘拟合，所得斜率即为广义hurst指数h(q)，计算得到尺度指数：τ(q)＝qh(q)-1；利用legendre变换，计算得到信号的奇异指数α和多分形奇异谱f(α)分别为α＝h(q)+qh'(q)、f(α)＝qα-τ(q)，其中h'(q)是广义hurst指数h(q)的legendre变换函数。