CN105426822B - 基于双树复小波变换的非平稳信号多重分形特征提取方法 - Google Patents
基于双树复小波变换的非平稳信号多重分形特征提取方法 Download PDFInfo
- Publication number
- CN105426822B CN105426822B CN201510744284.XA CN201510744284A CN105426822B CN 105426822 B CN105426822 B CN 105426822B CN 201510744284 A CN201510744284 A CN 201510744284A CN 105426822 B CN105426822 B CN 105426822B
- Authority
- CN
- China
- Prior art keywords
- scale
- signal
- dual
- wavelet
- fractal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于双树复小波变换的非平稳信号多重分形特征提取方法,步骤为:对待分析的非平稳信号进行集成处理;对集成信号进行双树复小波变换,利用小波分解尺度系数和细节系数得到各尺度下信号的波动成分;利用得到的各尺度小波系数,估算各尺度的瞬时频率,得到各尺度的时间尺度值的大小;结合尺度值大小,对各尺度下的波动成分进行分段;计算信号的各阶波动函数,利用波动函数和尺度值的双对数关系,经过最小二乘拟合得到广义hurst指数,得到各阶尺度指数;利用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尺度的小波系数进行单支重构,得到重构信号dl(t),其中1≤l≤M;对第M尺度的尺度系数的进行单支重构,得到重构信号cM(t);
d.每个尺度l的趋势表示为每个尺度l的波动成分表示为Fll(t)=y(t)-Trl(t)。
所述利用得到的各尺度小波系数,估算各尺度的瞬时频率,得到各尺度的时间尺度值的大小的方法的步骤为:
i.对第l尺度小波系数重构信号dl(t)进行希尔伯特变换,得到该尺度下的解析信号:其中,为dl(t)的希尔伯特变换;
ii.由解析信号zl(t)得到信号的相位角利用相位角的微分得到瞬时频率ωl(k),k=1,...,N/2l;
iii.利用瞬时频率ωl(k)得到对应尺度l的时间尺度的大小sl=1/<ωl(k)>,其中,<ωl(k)>为ωl(k)在该尺度的均值。
所述得到广义hurst指数、各阶尺度指数、多分形奇异谱的方法是:
①对第l尺度沿信号的正反两个方向,利用时间尺度sl对波动成分进行无覆盖的分段,共得到2Ns段,每段记为εv(i),i=1,...,sl;
②计算每段的局部波动函数
③取q∈[-qlim,+qlim],计算q=0之外的各阶波动函数
对q=0,波动函数为
④对logFq(sl)和logsl进行最小二乘拟合,所得斜率即为广义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尺度的小波系数进行单支重构,得到重构信号dl(t),其中1≤l≤M;对第M尺度的尺度系数的进行单支重构,得到重构信号cM(t);
d.每个尺度l的趋势表示为每个尺度l的波动成分表示为Fll(t)=y(t)-Trl(t)。
步骤3:利用得到的各尺度小波系数,估算各尺度的瞬时频率,得到各尺度的时间尺度值的大小。
利用步骤2得到的各尺度小波系数,估算各尺度的瞬时频率,得到各尺度的时间尺度值的大小的方法的步骤为:
i.对第l尺度小波系数重构信号dl(t)进行希尔伯特变换,得到该尺度下的解析信号:其中,为dl(t)的希尔伯特变换。
ii.由解析信号zl(t)得到信号的相位角利用相位角的微分得到瞬时频率ωl(k),k=1,...,N/2l。
iii.利用瞬时频率ωl(k)得到对应尺度l的时间尺度的大小sl=1/<ωl(k)>,其中,<ωl(k)>为ωl(k)在该尺度的均值。
步骤4:结合尺度值大小,对各尺度下的波动成分进行分段;计算信号的各阶波动函数,利用波动函数和尺度值的双对数关系,经过最小二乘拟合得到广义hurst指数,进而得到各阶尺度指数;利用legendre变换,得到信号的多分形奇异谱。
广义hurst指数、尺度指数、多分形奇异谱为信号的多重分形特征。得到广义hurst指数、各阶尺度指数、多分形奇异谱的方法是:
①对第l尺度沿信号的正反两个方向,利用时间尺度sl对波动成分进行无覆盖的分段,共得到2Ns段,每段记为εv(i),i=1,...,sl。
②计算每段的局部波动函数:
③取q∈[-qlim,+qlim],计算q=0之外的各阶波动函数
④对各阶波动函数Fq(sl)和时间尺度sl取对数,对logFq(sl)和logsl进行最小二乘拟合:logFq(sl)=h*logsl所得斜率h即为广义hurst指数h(q),计算得到尺度指数:τ(q)=qh(q)-1;利用legendre变换,计算得到信号的奇异指数α和多分形奇异谱f(α)分别为α=h(q)+qh'(q)、f(α)=qα-τ(q),h'(q)是广义hurst指数h(q)的legendre变换函数。
具体实例:
一种基于双树复小波变换的非平稳信号多重分形特征提取方法,对典型的非平稳信号p-模型乘性级联信号,模型参数p1=0.3和p2=0.7,信号长度为216,如图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尺度进行小波系数单支重构,得到信号细节成分dl(t);对第8尺度进行尺度系数的单支重构,得到信号近似成分c8(t)。
对1≤l≤8尺度,计算每个尺度的趋势成分:如图5所示;对应的波动成分为Fll(t)=y(t)-Trl(t),如图6所示。
步骤3:对1≤l≤8尺度的信号细节成分dl(t)进行希尔伯特变换,得到对应尺度下的解析信号其中,
那么瞬时频率
对1≤l≤8尺度,计算每个尺度对应的尺度大小其中,
步骤4:对1≤l≤8的每个尺度,沿波动成分的正反两个方向利用sl对其进行无覆盖的分段,共得到2Ns段,每段记为εv(i),i=1,...,sl。
利用计算每段的波动函数。在区间[-10,10]内取101个值作为阶数q的取值,利用下式计算各阶波动函数
对logFq(sl)和logsl进行最小二乘法拟合,直线的斜率即为广义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所述的基于双树复小波变换的非平稳信号多重分形特征提取方法,
其特征在于,所述对信号进行集成处理的方法是:
其中,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尺度的小波系数进行单支重构,得到重构信号dl(t),其中1≤l≤M;对第M尺度的尺度系数的进行单支重构,得到重构信号cM(t);
d.每个尺度l的趋势表示为每个尺度l的波动成分表示为Fll(t)=y(t)-Trl(t)。
4.根据权利要求1所述的基于双树复小波变换的非平稳信号多重分形特征提取方法,其特征在于,所述利用得到的各尺度小波系数,估算各尺度的瞬时频率,得到各尺度的时间尺度值的大小的方法的步骤为:
i.对第l尺度小波系数重构信号dl(t)进行希尔伯特变换,得到该尺度下的解析信号:其中,的希尔伯特变换;
ii.由解析信号zl(t)得到信号的相位角利用相位角的微分得到瞬时频率ωl(k),k=1,...,N/2l;
iii.利用瞬时频率wl(k)得到对应尺度l的时间尺度的大小sl=1/<wl(k)>,其中,<ωl(k)>为ωl(k)在该尺度的均值。
5.根据权利要求1所述的基于双树复小波变换的非平稳信号多重分形特征提取方法,其特征在于,所述得到广义hurst指数、各阶尺度指数、多分形奇异谱的方法是:
①对第l尺度沿信号的正反两个方向,利用时间尺度sl对波动成分进行无覆盖的分段,共得到2Ns段,每段记为εv(i),i=1,...,sl;
②计算每段的局部波动函数
③取q∈[-qlim,+qlim],计算q=0之外的各阶波动函数
对q=0,波动函数为
④对logFq(sl)和logsl进行最小二乘拟合,所得斜率即为广义hurst指数h(q),计算得到尺度指数:τ(q)=qh(q)-1;利用legendre变换,计算得到信号的奇异指数α和多分形奇异谱f(α)分别为α=h(q)+qh'(q)、f(α)=qα-τ(q),其中h'(q)是广义hurst指数h(q)的legendre变换函数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510744284.XA CN105426822B (zh) | 2015-11-05 | 2015-11-05 | 基于双树复小波变换的非平稳信号多重分形特征提取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510744284.XA CN105426822B (zh) | 2015-11-05 | 2015-11-05 | 基于双树复小波变换的非平稳信号多重分形特征提取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105426822A CN105426822A (zh) | 2016-03-23 |
CN105426822B true CN105426822B (zh) | 2018-09-04 |
Family
ID=55505022
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510744284.XA Active CN105426822B (zh) | 2015-11-05 | 2015-11-05 | 基于双树复小波变换的非平稳信号多重分形特征提取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105426822B (zh) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105615879B (zh) * | 2016-04-05 | 2018-07-17 | 陕西师范大学 | 基于多重分形消除趋势波动分析的脑电测谎方法 |
CN106092879B (zh) * | 2016-06-07 | 2019-07-12 | 西安向阳航天材料股份有限公司 | 基于振动响应信息的爆炸复合管结合状态检测方法 |
CN106770675B (zh) * | 2016-12-06 | 2019-05-21 | 郑州轻工业学院 | 基于声发射信号的金刚石压机顶锤裂纹在线检测方法 |
CN107411739A (zh) * | 2017-05-31 | 2017-12-01 | 南京邮电大学 | 基于双树复小波的脑电信号情绪识别特征提取方法 |
CN108181098A (zh) * | 2017-12-15 | 2018-06-19 | 天津金岸重工有限公司 | 一种门座式起重机低速重载部件故障特征提取方法 |
KR20210091737A (ko) * | 2018-11-09 | 2021-07-22 | 아규리 시스템스 엘티디. | 비정상 기계 성능의 자동화된 분석 |
CN110160791B (zh) * | 2019-06-27 | 2021-03-23 | 郑州轻工业学院 | 基于小波-谱峭度的感应电机轴承故障诊断系统及诊断方法 |
CN111256965B (zh) * | 2020-01-20 | 2022-03-11 | 郑州轻工业大学 | 多尺度信息融合的堆叠稀疏自编码旋转机械故障诊断方法 |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7450779B2 (en) * | 2004-05-21 | 2008-11-11 | Imaging Dynamics Company Ltd. | De-noising digital radiological images |
EP1788937A4 (en) * | 2004-09-16 | 2009-04-01 | Brainscope Co Inc | METHOD FOR THE ADAPTIVE, COMPLEX, WAVELET-BASED FILTERING OF EEG SIGNALS |
CN102411708A (zh) * | 2011-12-02 | 2012-04-11 | 湖南大学 | 一种融合双树复小波变换和离散小波变换的人脸识别方法 |
CN104200436B (zh) * | 2014-09-01 | 2017-01-25 | 西安电子科技大学 | 基于双树复小波变换的多光谱图像重构方法 |
-
2015
- 2015-11-05 CN CN201510744284.XA patent/CN105426822B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN105426822A (zh) | 2016-03-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105426822B (zh) | 基于双树复小波变换的非平稳信号多重分形特征提取方法 | |
Zhang et al. | A novel hybrid model based on VMD-WT and PCA-BP-RBF neural network for short-term wind speed forecasting | |
CN105572501B (zh) | 一种基于sst变换和ls-svm的电能质量扰动识别方法 | |
CN104572886B (zh) | 基于k线图表示的金融时间序列相似性查询方法 | |
CN103163050B (zh) | 一种基于电磁感应信号的润滑油系统金属磨粒检测方法 | |
CN106501602B (zh) | 一种基于滑窗频谱分离的基波参数测量方法 | |
Paul et al. | Wavelet based denoising technique for liquid level system | |
CN102288843A (zh) | 一种电能质量扰动信号检测方法 | |
CN106980044B (zh) | 一种适应风电接入的电力系统谐波电流估计方法 | |
CN110081967A (zh) | 基于谱图小波变换的机械振动信号阈值降噪方法 | |
Huda et al. | Power quality signals detection using S-transform | |
CN112268615B (zh) | 一种机电设备振动信号特征提取方法 | |
Tong et al. | Parameter estimation of FH signals based on STFT and music algorithm | |
Bai et al. | Application of Time‐Frequency Analysis in Rotating Machinery Fault Diagnosis | |
CN106788092A (zh) | 一种基于原子分解法的同步电机参数辨识方法 | |
CN106137184B (zh) | 基于小波变换的心电信号qrs波检测方法 | |
CN107491412A (zh) | 一种基于经验小波变换的用户用电负荷特征提取方法 | |
CN112798888B (zh) | 一种无人驾驶列车车载电气系统故障非侵入诊断方法 | |
CN108090270B (zh) | 一种基于形态学滤波和盲源分离的暂态振荡参数识别方法 | |
CN105606892A (zh) | 一种基于sst变换的电网谐波与间谐波分析方法 | |
CN103345739A (zh) | 一种基于纹理的高分辨率遥感影像建筑区指数计算方法 | |
Zheng et al. | Multi-gradient surface electromyography (SEMG) movement feature recognition based on wavelet packet analysis and support vector machine (SVM) | |
CN105572472A (zh) | 分布式电源环境的频率测量方法与系统 | |
Qiwei et al. | EMD algorithm based on bandwidth and the applicationon one economic data analysis | |
Zhang et al. | Identification inrush current and internal faults of transformer based on Hyperbolic S-transform |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |