CN108169560A - 一种分段正弦拟合分解方法 - Google Patents

Abstract

本发明涉及一种分段正弦拟合分解方法，步骤一：初始化循环变量k、最大迭代次数Q和收敛阈值ε，令k＝1,r_k(t)＝x(t)；步骤二：寻找迭代信号r_k(t)所有过零点；步骤三：寻找相邻过零点之间局部最大值点和局部最小值点；步骤四：拟合中值曲线m(t)，从迭代信号r_k(t)中减去中值曲线m(t)，得到临时信号r₀(t)；步骤五：对临时信号r₀(t)进行分段正弦拟合，获得拟合函数s_k(t)；步骤六：令r_k+1(t)＝r_k(t)‑s_k(t)；步骤七：令k＝k+1；步骤八：如果k>Q或条件成立，则算法迭代过程结束，否则重复执行步骤二至步骤八，本发明不需要筛选过程来实现极值点个数和过零点数的约束条件，只需减除一次中值曲线，通过分段的正弦曲线进行拟合，有效降低计算量，瞬时频率曲线可直接从分解结果得到，简化了计算步骤。

Description

一种分段正弦拟合分解方法

技术领域

本发明涉及一种非平稳信号的分析方法，特别涉及一种分段正弦拟合分解方法。

背景技术

在自然界以及声、光、电、磁等工程领域存在大量非线性的、时变的、非平稳的信号，非平稳信号的主要特征是其时变性，其频率随时间而改变。以傅立叶变换为基础的经典分析方法是包括线性、高斯、周期、平稳假设条件的信号分析的最有力的工具，它能够反映信号整个时间域内的频率分布情况，但不能有效反映信号的时变频谱特性，无法获得许多关键信息，难以满足实际应用的需要。为了改善傅立叶方法的不足，短时傅立叶变换、小波分析、Wigner Ville分布、S变换等时频分析方法被用来分析时变信号，它们在不同程度上描述了非平稳信号的频率随时间变化的关系，但它们最终是基于傅立叶分析的，都属于全局分析方法的范畴，存在各自的局限性，不是受线性约束，就是受平稳性束缚，不能完全意义上处理非线性非平稳信号。例如短时傅立叶变换由于受Heisenberg测不准原理约束，不能同时获得较高的时频分辨率，Wigner Ville等二次型时频分布，存在交叉项干扰，固定核函数只适合特定信号。

Hilbert-Huang变换(HHT)是一种数据自适应的局域分析方法，它使用经验模态分解方法EMD(Empirical Mode Decomposition)，把复杂数据自适应分解成有限多个内蕴模态函数IMFS(Intrinsic Mode Functions)分量。它最突出的优点是分解过程不依赖特定的分解基，完全自适应数据本身，普适性好。另外，对IMF分量进行Hilbert变换，可以获得时间-频率-能量的三维谱图，它不受Heisenberg测不准原理制约，可以在时间和频率同时达到较高的分辨率，因此，HHT理论在分析非线性非平稳数据方面取得了优势性的进展。然而也存在一些局限：一是内蕴模态函数必须满足极值点和过零点数量相等或最多相差一个的条件，而极值点和过零点的数量是在一个时间窗口内统计的，上下包络的计算也是在窗口内进行的，因此对应的EMD算法并不是完全意义的局部波算法；二是EMD的筛选次数难以确定，筛选过程极值统计和包络拟合的计算量也比较大；三是Hilbert变换的谱分析方法通过积分变换计算解析信号，根据解析信号相位函数的导数计算的瞬时频率，然而积分的全局变换是一种滑动平均意义的处理，不能真正的反映信号瞬时的改变。

发明内容

为克服上述技术缺陷，本发明要解决的技术问题是提供一种基于数据自适应的分段正弦拟合分解方法。

为实现上述目的，本发明的技术方案是：

步骤一：初始化循环变量k、最大迭代次数Q和收敛阈值ε，令k＝1,r_k(t)＝x(t)，其中，x(t)为原信号，r_k(t)为迭代信号；

步骤二：寻找迭代信号r_k(t)的所有过零点，其中过零点对应时间序列满足：

T_Zero＝{tz_i|i∈[1,K],tz_i-1＜tz_i}

其中T_Zero代表迭代信号r_k(t)中过零点时间索引的集合,tz_i代表第i个过零点时间，K代表信号中过零点的个数；

过零点对应幅度序列满足：

V_Zero＝{az_i|az_i＝r_k(tz_i),i∈[1,K],tz_i∈T_zero}

其中r_k(tz_i)＝0代表过零点处的函数值为零，az_i代表第i个过零点的幅度；

步骤三：寻找相邻过零点之间的局部最大值点和局部最小值点，所述的局部最大值点是当相邻过零点之间的值都大于零时相邻过零点之间的最大值点，所述的局部最小值点是相邻过零点之间的值都小于零时相邻过零点之间的最小值点，所述的局部最大值点对应的时间索引序列满足：

T_Max＝{tx_i|i∈[1,L],tx_i＜tx_i+1}

其中T_Max代表信号r_k(t)中局部最大值点时间索引的集合,tx_i代表第i个局部最大值点对应的时间，L代表局部极大值点个数；

所述的局部最大值点对应的幅度序列满足：

V_Max＝{ax_i|ax_i＝r_k(tx_i),i∈[1,L],tx_i∈T_Max}

其中r_k(tx_i)代表信号局部最大值时刻tx_i对应的值，ax_i代表第i个局部最大值点的幅度；

所述的局部最小值点对应的时间索引序列满足：

T_Min＝{tn_i|i∈[1,M],tn_i＜tn_i+1}

其中T_Min代表信号r_k(t)中局部最小值点时间索引的集合,tn_i代表第i个局部最小值点对应的时间，M代表局部最小值点的个数；

所述的局部最小值点对应的幅度序列满足：

V_Min＝{an_i|an_i＝r_k(tn_i),i∈[1,M],tn_i∈T_Min}

其中r_k(tn_i)代表信号局部最小值时刻tn_i对应的值，an_i代表第i个局部最小值点的幅度；

过零点个数、局部极大值点个数和局部极小值点个数满足：

|K-(M+L)|＝1or 0

将过零点时间、局部最大值点时间和局部最小值点时间放在一起，并按时间先后排序，获得时间序列，即：

T_All＝{t_i|i∈[1,M+L+K],t_i＜t_i+1}＝{T_Zero,T_Min,T_Max}

按照时间序列T_All，可得到幅度序列：

V_All＝{a_i|a_i＝r_k(t_i),i∈[1,M+L+K],t_i∈T_All}＝{V_Zero,V_Min,V_Max}

其中t_i代表T_All中第i个时间，为过零点时间或者局部最大值点时间或者局部最小值点时间，t₁或t_M+L+K对应过零点时间，a_i是对应T_All中第i个时间的信号幅度，对应时间序列T_All的瞬时相位序列满足：

其中为对应时刻t_i的瞬时相位,如果a₂>0,则否则

步骤四：拟合中值曲线m(t)，从迭代信号r_k(t)中减去中值曲线m(t)，得到临时信号r₀(t)，临时信号r₀(t)满足：r₀(t)＝r_k(t)-m(t)，其中中值曲线m(t)拟合过程满足：

第一类点：用直线连接相邻的局部极大值点和局部极小值点，取其中点；第二类点：用直线连接相邻的两个局部极大值点，过所述相邻的两个局部极大值点之间的局部极小值点做平行于y轴的直线并与连接相邻的两个局部极大值点的直线相交，取交点和局部极小值点连线的中点；第三类点：用直线连接相邻的两个局部极小值点，过所述的相邻两个局部极小值点之间的局部极大值点做平行于y轴的直线并与连接相邻的两个局部极小值点的直线相交，取交点和局部极大值点连线的中点；以所述第一、二、三类点为节点进行插值，拟合信号的中值曲线m(t)；

步骤五：对所述的临时信号r₀(t)进行分段正弦拟合，获得拟合函数s_k(t)，所述的拟合函数s_k(t)满足：

其中a_i∈V_All,t_i∈T_All,1≤i≤K,且k∈N是自然数,K是过零点个数，如果信号的左右端点不为零，则：

瞬时相位函数满足：

瞬时频率是相位函数的一阶导数，表示为：

对应的瞬时幅度函数表示为：

a_k(t)＝{|a_2i|,t∈[t_2i-1,t_2i+1),i∈[1,K]}

步骤六：令r_k+1(t)＝r_k(t)-s_k(t)

步骤七：令k＝k+1；

步骤八：如果k>Q或条件成立，则算法迭代过程结束，否则重复执行步骤二至步骤八，迭代过程结束后，原信号满足：

其中：Q为分解阶数，与最大迭代次数相同；s_k(t)为第k阶分段正弦拟合函数，r_k+1(t)为分解后的残差函数。

本发明还包括：

步骤四中所述插值方法包括线性插值法、多项式插值法、样条插值法。

本发明的有益效果：本发明提供了一种新的基于数据自适应的非平稳信号时频分析方法，定义了时间尺度、局部极值和中值曲线的新含义，并给出了新的中值曲线拟合的方法。通过发明的分段正弦拟合分解PSFD(Piecewise Sinusoidal FittingDecomposition)，能够将信号分解成有限多个分段正弦基函数PSBFs(PiecewiseSinusoidal(PWS)Basis Functions)和残差函数的集合。该算法采用新的时间尺度、局部极值和中值曲线的定义，可以不需要筛选过程，瞬时频率也可以直接给出，大大降低了计算成本，且PSBFs的能量随分解阶数由高到低分布。与EMD方法相比，不需要筛选过程来实现极值点个数和过零点数的约束条件，只需减除一次中值曲线，然后通过分段的正弦曲线进行拟合，有效降低了计算量，瞬时频率曲线可直接从分解结果得到，不需要通过Hilbert变换获取，简化了计算步骤，且PSBFs的能量随分解阶数由高到低分布，分解的时间尺度由大到小分布。

附图说明

图1为准周期信号的时间尺度和极值；

图2为中值曲线拟合图；

图3为分段正弦拟合分解流程图；

图4为风噪声信号的分解结果；

图5为PWSFs的瞬时频率曲线图；

图6为PWSFs随阶数的能量分布曲线；

图7为时频功率谱分布；

图8为加权平均时频功率谱分布。

具体实施方式

下面结合附图对技术方案的实施作进一步的详细描述：

(一)实施方式一：

步骤一：初始化循环变量k、最大迭代次数Q和收敛阈值ε，令k＝1,r_k(t)＝x(t)，其中，x(t)为原信号，r_k(t)为迭代信号；

步骤二：寻找迭代信号r_k(t)的所有过零点，其中过零点对应时间序列满足：

T_Zero＝{tz_i|i∈[1,K],tz_i-1＜tz_i} (1)

其中T_Zero代表迭代信号r_k(t)中过零点时间索引的集合，tz_i代表第i个过零点时间，K代表信号中过零点的个数；

过零点对应幅度序列满足：

V_Zero＝{az_i|az_i＝r_k(tz_i),i∈[1,K],tz_i∈T_Zero} (2)

其中r_k(tz_i)＝0代表过零点处的函数值为零，az_i代表第i个过零点的幅度；

步骤三：寻找相邻过零点之间的局部最大值点和局部最小值点，局部最大值点是当相邻过零点之间的值都大于零时相邻过零点之间的最大值点，局部最小值点是相邻过零点之间的值都小于零时相邻过零点之间的最小值点，局部最大值点对应的时间索引序列满足：

T_Max＝{tx_i|i∈[1,L],tx_i＜tx_i+1} (3)

其中T_Max代表迭代信号r_k(t)中局部最大值点时间索引的集合,tx_i代表第i个局部最大值点对应的时间，L代表局部极大值点个数；

局部最大值点对应的幅度序列满足：

V_Max＝{ax_i|ax_i＝r_k(tx_i),i∈[1,L],tx_i∈T_Max} (4)

其中r_k(tx_i)代表信号局部最大值时刻tx_i对应的值，ax_i代表第i个局部最大值点的幅度；

局部最小值点对应的时间索引序列满足：

T_Min＝{tn_i|i∈[1,M],tn_i＜tn_i+1} (5)

其中T_Min代表信号r_k(t)中局部最小值点时间索引的集合,tn_i代表第i个局部最小值点对应的时间，M代表局部最小值点的个数；

局部最小值点对应的幅度序列满足：

V_Min＝{an_i|an_i＝r_k(tn_i),i∈[1,M],tn_i∈T_Min} (6)

其中r_k(tn_i)代表信号局部最小值时刻tn_i对应的值，an_i代表第i个局部最小值点的幅度；

过零点个数、局部极大值点个数和局部极小值点个数满足：

|K-(M+L)|＝1or 0 (7)

将过零点时间、局部最大值点时间和局部最小值点时间放在一起，并按时间先后排序，获得时间序列，即：

T_All＝{t_i|i∈[1,M+L+K],t_i＜t_i+1}＝{T_Zero,T_Min,T_Max} (8)

按照时间序列T_All，可得到幅度序列：

V_All＝{a_i|a_i＝r_k(t_i),i∈[1,M+L+K],t_i∈T_All}＝{V_Zero,V_Min,V_Max} (9)

其中t_i代表T_All中第i个时间，为过零点时间或者局部最大值点时间或者局部最小值点时间，t₁或t_M+L+K对应过零点时间，a_i是对应T_All中第i个时间的信号幅度，对应时间序列T_All的瞬时相位序列满足：

其中为对应时刻t_i的瞬时相位,如果a₂>0,则否则

步骤四：拟合中值曲线m(t)，从迭代信号r_k(t)中减去中值曲线m(t)，得到临时信号r₀(t)，临时信号r₀(t)满足：r₀(t)＝r_k(t)-m(t)，其中中值曲线m(t)拟合过程满足：

第一类点：用直线连接相邻的局部极大值点和局部极小值点，取其中点；第二类点：用直线连接相邻的两个局部极大值点，过所述相邻的两个局部极大值点之间的局部极小值点做平行于y轴的直线并与连接相邻的两个局部极大值点的直线相交，取交点和局部极小值点连线的中点；第三类点：用直线连接相邻的两个局部极小值点，过相邻两个局部极小值点之间的局部极大值点做平行于y轴的直线并与连接相邻的两个局部极小值点的直线相交，取交点和局部极大值点连线的中点；以所述第一、二、三类点为节点进行插值，拟合信号的中值曲线m(t)；

步骤五：对临时信号r₀(t)进行分段正弦拟合，获得拟合函数s_k(t)，拟合函数s_k(t)满足：

其中a_i∈V_All,t_i∈T_All,1≤i≤K,且k∈N是自然数,K是过零点个数，如果信号的左右端点不为零，则：

瞬时相位函数满足：

瞬时频率是相位函数的一阶导数，表示为：

对应的瞬时幅度函数表示为：

a_k(t)＝{|a_2i|,t∈[t_2i-1,t_2i+1),i∈[1,K]} (15)

步骤六：令r_k+1(t)＝r_k(t)-s_k(t)

步骤七：令k＝k+1；

步骤八：如果k>Q或条件成立，则算法迭代过程结束，否则重复执行步骤二至步骤八，迭代过程结束后，原信号满足：

其中：Q为分解阶数，与最大迭代次数相同；s_k(t)为第k阶分段正弦拟合函数，r_k+1(t)为分解后的残差函数。

步骤四中的插值方法包括线性插值法、多项式插值法、样条插值法。

(二)实施方式二

实施实例是一段风噪声信号，如图4最上面的子图x(t)所示，它是待分解的原始信号，用x(t)表示。图2所示为从待分解信号x(t)中截取的一段作为实施实例进行说明，则对x(t)进行分段正弦拟合分解的算法流程如图3所示，具体如下：

步骤一：初始化循环变量k、最大迭代次数Q和收敛阈值ε，令k＝1,r_k(t)＝x(t)，其中，x(t)为原信号，r_k(t)为迭代信号；

步骤二：寻找迭代信号r_k(t)的所有过零点，其中过零点对应时间序列满足：

T_Zero＝{tz_i|i∈[1,K],tz_i-1＜tz_i}

其中T_Zero代表迭代信号r_k(t)中过零点时间索引的集合，tz_i代表第i个过零点时间，K代表信号中过零点的个数；

过零点对应幅度序列满足：

V_Zero＝{az_i|az_i＝r_k(tz_i),i∈[1,K],tz_i∈T_Zero}

其中r_k(tz_i)＝0代表过零点处的函数值为零，az_i代表第i个过零点的幅度；

如图2所示的线1为待分解信号波形，Z_2i-1,Z_2i,Z_2i+1,Z_2i+2,Z_2i+3是待分解信号对应的过零点，它的横坐标就是对应的过零点时间，纵坐标为零。

步骤三：寻找相邻过零点之间的局部最大值点和局部最小值点，局部最大值点是当相邻过零点之间的值都大于零时相邻过零点之间的最大值点，局部最小值点是相邻过零点之间的值都小于零时相邻过零点之间的最小值点，局部最大值点对应的时间索引序列满足：

T_Max＝{tx_i|i∈[1,L],tx_i＜tx_i+1}

其中T_Max代表迭代信号r_k(t)中局部最大值点时间索引的集合,tx_i代表第i个局部最大值点对应的时间，L代表局部极大值点个数；

局部最大值点对应的幅度序列满足：

V_Max＝{ax_i|ax_i＝r_k(tx_i),i∈[1,L],tx_i∈T_Max}

其中r_k(tx_i)代表信号局部最大值时刻tx_i对应的值，ax_i代表第i个局部最大值点的幅度；

如图2所示的P_i-1,P_i,P_i+1对应实例曲线的局部最大值点，他们的横坐标分别对应局部最大值时间索引，纵坐标对应局部最大值点的幅度。

局部最小值点对应的时间索引序列满足：

T_Min＝{tn_i|i∈[1,M],tn_i＜tn_i+1}

其中T_Min代表信号r_k(t)中局部最小值点时间索引的集合,tn_i代表第i个局部最小值点对应的时间，M代表局部最小值点的个数；

局部最小值点对应的幅度序列满足：

V_Min＝{an_i|an_i＝r_k(tn_i),i∈[1,M],tn_i∈T_Min}

其中r_k(tn_i)代表信号局部最小值时刻tn_i对应的值，an_i代表第i个局部最小值点的幅度；

如图2所示的V_i-1,V_i,V_i+1分别对应实例曲线的局部最小值点，他们的横坐标分别对应局部最小值时间索引，纵坐标对应局部最小值点的幅度。

过零点个数、局部极大值点个数和局部极小值点个数满足：

|K-(M+L)|＝1or 0

如图2所示，该段信号有4个局部最大值点，3个局部最小值点，7个过零点，满足上式的要求。

将过零点时间、局部最大值点时间和局部最小值点时间放在一起，并按时间先后排序，获得时间序列，即：

T_All＝{t_i|i∈[1,M+L+K],t_i＜t_i+1}＝{T_Zero,T_Min,T_Max}

按照时间序列T_All，可得到幅度序列：

V_All＝{a_i|a_i＝r_k(t_i),i∈[1,M+L+K],t_i∈T_All}＝{V_Zero,V_Min,V_Max}

其中t_i代表T_All中第i个时间，为过零点时间或者局部最大值点时间或者局部最小值点时间，t₁或t_M+L+K对应过零点时间，a_i是对应T_All中第i个时间的信号幅度，对应时间序列T_All的瞬时相位序列满足：

其中为对应时刻t_i的瞬时相位,如果a₂>0,则否则

如图2所示，将V_i-1,V_i,V_i+1,P_i-1,P_i,P_i+1,Z_2i-1,Z_2i,Z_2i+1,Z_2i+2,Z_2i+3对应的点按照时间顺序先后排序，各点相位按照初始相位，每隔一个点增加0.5π。

步骤四：拟合中值曲线m(t)，从迭代信号r_k(t)中减去中值曲线m(t)，得到临时信号r₀(t)，临时信号r₀(t)满足：r₀(t)＝r_k(t)-m(t)，其中中值曲线m(t)拟合过程满足：

第一类点：用直线连接相邻的局部极大值点和局部极小值点，取其中点；第二类点：用直线连接相邻的两个局部极大值点，过所述相邻的两个局部极大值点之间的局部极小值点做平行于y轴的直线并与连接相邻的两个局部极大值点的直线相交，取交点和局部极小值点连线的中点；第三类点：用直线连接相邻的两个局部极小值点，过所述的相邻两个局部极小值点之间的局部极大值点做平行于y轴的直线并与连接相邻的两个局部极小值点的直线相交，取交点和局部极大值点连线的中点；以所述第一、二、三类点为节点进行插值，拟合信号的中值曲线m(t)；

如图2所示，用线段连接相邻的局部极大值和极小值点，即线段P_i-1V_i-1，V_i-1P_i，P_iV_i，V_kP_k+1，分别过局部极大值点和局部极小值点做平行y轴的直线，与对应相邻局部极小值或极大值的连线相交，得到线段V_i-1C_2(i-1)+1,P_iC_2(i-1)+2,V_iC_2i+1,P_i+1C_2i-+2,取这些线段的中点M_3(i-1)+1,M_3(i-1)+2,M_3(i-1)+3,M_3i+1,M_3i+2,M_3i+3,M_3(i+1)+1…作为基准点进行插值，图2所示的线2为插值后的中值曲线。

步骤五：对所述的临时信号r₀(t)进行分段正弦拟合，获得拟合函数s_k(t)，所述的拟合函数s_k(t)满足：

其中a_i∈V_All,t_i∈T_All,1≤i≤K,且k∈N是自然数,K是过零点个数。如图2实例所示第一段准周期(点P_i-1Z_2i-1Z_2iP对应的段)的分段正弦函数可写为：

瞬时相位函数满足：

如图2实例所示第一个段准周期的瞬时相位为：

瞬时频率是相位函数的一阶导数，表示为：

如图2实例所示第一个段准周期的瞬时频率为：

对应的瞬时幅度函数表示为：

a_k(t)＝{|a_2i|,t∈[t_2i-1,t_2i+1),i∈[1,K]}

如图2实例所示第一个段准周期对应的瞬时幅度为：

步骤六：令r_k+1(t)＝r_k(t)-s_k(t)

步骤七：令k＝k+1；

步骤八：如果k>Q或条件成立，则算法迭代过程结束，否则重复执行步骤二至步骤八，迭代过程结束后，原信号满足：

其中：Q为分解阶数，与最大迭代次数相同；s_k(t)为第k阶分段正弦拟合函数，r_k+1(t)为分解后的残差函数。

步骤四中的插值方法包括线性插值法、多项式插值法、样条插值法。

如图4所示为风噪声分段正弦拟合分解的结果，停止参数残差能量与原信号能量之比小于0.01。如图5所示为1-3阶分段正弦拟合函数的瞬时频率随时间变化的曲线。用公式(17)可计算各阶的能量分布情况，如图6所示为随阶数变化的对数能量分布曲线，能量随阶数降低。

结合瞬时频率和瞬时幅度公式(14-15)，可给出信号的二维时频功率谱，如图7所示，横纵坐标分别是时间和频率，灰度是对应频率的瞬时功率。

根据公式(18)用瞬时幅度对同一时刻的频率进行加权平均，可得加权平均瞬时频率曲线，根据公式19计算平均瞬时幅度，则平均时频分布可在二维时频图中画出，如图8所示，可观察主要信号频率随时间的变化情况。

(三)实施方式三

时间尺度是信号分析的重要参数，单分量三角函数的时间尺度就是它的周期，但它只定义了时间尺度的全局均值，无论在幅度和频率上都脱离了随时间变化的事实。本发明采用如下的定义:

定义1：局部极值包括局部最大值和局部最小值。局部极值是相邻过零点之间的最大值或最小值(绝对值最大的值)。换句话说，如果相邻过零点之间的值都大于零，则局部最大值就是过零点之间的最大值。如果相邻过零点之间的值都小于零，则局部最小值就是过零点之间的最小值。局部最大值和局部最小值交替出现。

定义2：用过零点和局部极值点之间的时间间隔描述局部时间尺度,在类谐波的一个准周期内包括四段：过零点到局部极大值点的时间间隔τ₁，局部极大值点到过零点的时间间隔它τ₂，过零点到局部极小值τ₃，局部极小值到过零点τ₄。

定义3：用直线连接相邻的时间尺度内的极大值和极小值，取其中点，相邻的两个极大值区间一定有一个极小值区间，用直线连接这两个极大值，过极小值做平行于y轴的直线并与极大直线相交，取交点和极小值连线的中点。相邻的极小值区间也一定有一个极大值区间，同样方法获得一个中点。以这些中点为节点进行插值，拟合信号的中值曲线m(t)。

对于给定信号x(t),t∈[0,T]，其中T代表信号时间长度,对应的过零点时间序列表示为：

T_Zero＝{tz_i|i∈[1,K],tz_i-1＜tz_i} (1)

其中T_Zero代表信号x(t)中过零点时间索引的集合,tz_i代表第i个过零点时间，K代表信号中过零点的个数。局部最大值时间序列为：

T_Max＝{tx_i|i∈[1,L],tx_i＜tx_i+1} (2)

其中T_Max代表信号x(t)中局部最大值时间索引的集合,tx_i代表第i个局部最大值对应的时间，L代表局部极大值个数，这里的局部最大值是定义1中的含义。局部最小值时间序列为：

T_Min＝{tn_i|i∈[1,M],tn_i＜tn_i+1} (3)

其中T_Min代表信号x(t)中局部最小值时间索引的集合,tn_i代表第i个局部最小值对应的时间，M代表局部最小值的个数，这里的局部最小值是定义1中的含义。过零点序列为：

V_Zero＝{az_i|az_i＝x(tz_i),i∈[1,K],tz_i∈T_zero} (4)

其中x(tz_i)＝0代表过零点处的函数值为零。局部最大值序列为：

V_Max＝{ax_i|ax_i＝x(tx_i),i∈[1,L],tx_i∈T_Max} (5)

其中x(tx_i)代表信号局部最大值时刻tx_i对应的值。

局部最小值序列为：

V_Min＝{an_i|an_i＝x(tn_i),i∈[1,M],tn_i∈T_Min} (6)

上述序列，过零点个数和极值点个数满足如下条件：

|K-(M+L)|＝1or 0 (7)

将过零点时间、局部最大值时间和局部最小值时间放在一起，并按时间先后排序，即：

T_All＝{t_i|i∈[1,M+L+K],t_i＜t_i+1}＝{T_Zero,T_Min,T_Max} (8)

按照时间序列T，可得到幅度序列：

V_All＝{a_i|a_i＝x(t_i),i∈[1,M+L+K],t_i∈T_All}＝{V_Zero,V_Min,V_Max} (9)

其中t_i代表T_All中第i个时间，可能是过零点时间，也可能时极值时间，t₁或t_M+L+K对应零点时间。a_i是对应T_All中第i个时间的信号幅度。对应时间序列T_All的瞬时相位序列表示为：

其中为对应时刻t_i的瞬时相位,如果a₂>0,则否则

定义4：按照定义1和定义2，将信号的一个准周期分段，每段之间的相位差为0.5π，用分段正弦函数进行拟合，其表达式如下：

其中a_k∈V_All,t_k∈T_All,1≤k≤K,且k∈N是自然数,K是过零点个数。如果信号的左右端点不为零，则

对应的瞬时相位函数可以表示为：

瞬时频率是相位函数的一阶导数，表示为：

对应的瞬时幅度函数表示为：

a(t)＝a_2k,t∈[t_2k-1,t_2k+1) (15)

信号左右端点附近的瞬时相位、瞬时频率和瞬时幅度可根据公式(12)，采用类似的方法提取。

算法分解过程如下:

1)初始化最大迭代次数Q和收敛阈值ε，令k＝1,r₀(t)＝x(t)；

2)对信号r₀(t)重复如下操作:

a)寻找信号r₀(t)的所有过零点，对应的过零点时间序列用公式(1)表示。按照定义1，寻找相邻过零点之间的局部最大值和局部最小值点，用公式(2)和公式(3)表示对应的时间索引序列，用公式(5)和公式(6)表示极值点的幅度值；

b)按照定义3拟合中值曲线m(t)，并从信号减去中值曲线，得到剩余信号

r₀(t)＝r₀(t)-m(t)；

c)对b)中剩余信号r₀(t)用定义4的方法进行分段正弦拟合，获得拟合函数s_k(t)；

d)从c)中的剩余信号r₀(t)中减去拟合函数s_k(t)，即r₀(t)＝r₀(t)-s_k(t)；

e)令k＝k+1；

f)如果k>Q或条件成立，则算法迭代过程结束。

经过迭代处理，原信号可以分解为多个分段正弦拟合函数和残差函数的和，即：

其中，Q是分解阶数，s_k(t)为第k阶分段正弦拟合函数。由公式(13)～(15)可计算每一阶正弦拟合函数的瞬时相位函数、瞬时幅度和瞬时频率函数。

(四)实施方式四

如图1所示，实施实例中信号是一个准周期信号，根据定义1，B和D分别是相邻过零点之间的局部极大值和局部极小值，而且在相邻过零点之间有且只有一个极值(局部极值中绝对值最大的极值)。根据定义2一个准周期可以分成四段，时间尺度分别为：τ₁＝t₂-t₁,τ₂＝t₃–t₂,τ₃＝t₄-t₃,τ₄＝t₅–t₄。

如图2所示，实施实例是一段风噪声信号，图中，P_k-1,P_k,P_k+1和V_k-1,V_k,V_k+1分别是定义1中所定义的局部极大值和局部极小值。根据定义3，用线段连接相邻的局部极大值和极小值点，即线段P_k-1V_k-1，V_k-1P_k，P_kV_k，V_kP_k+1，分别过局部极大值点和局部极小值点做平行y轴的直线，与对应相邻局部极小值或极大值的连线相交，得到线段V_k-1C_2(k-1)+1,P_kC_2(k-1)+2,V_kC_2k+1,P_k+1C_2k-+2,取这些线段的中点M_3(k-1)+1,M_3(k-1)+2,M_3(k-1)+3,M_3k+1,M_3k+2,M_3k+3,M_3(k+1)+1…作为基准点进行插值，插值方法可根据情况选择，例如线性插值法、多项式插值法或样条插值法等。

如图3所示，分段正弦拟合分解的算法流程图，分解步骤如下：

1)初始化最大迭代次数Q和收敛阈值ε，令k＝1,r₀(t)＝x(t)；

2)对信号r₀(t)重复如下操作:

a)寻找信号r₀(t)的所有过零点，对应的过零点时间序列用公式(1)表示。按照定义1，寻找相邻过零点之间的局部最大值和局部最小值点，用公式(2)和公式(3)表示对应的时间索引序列，用公式(5)和公式(6)表示极值点的幅度值；

b)按照定义3拟合中值曲线m(t)，并从信号减去中值曲线，得到剩余信号

r₀(t)＝r₀(t)-m(t)；

c)对b)中剩余信号r₀(t)用定义4的方法进行分段正弦拟合，获得拟合函数s_k(t)；

d)从c)中的剩余信号r₀(t)中减去拟合函数s_k(t)，即r₀(t)＝r₀(t)-s_k(t)；

e)令k＝k+1；

f)如果k>Q或条件成立，则算法迭代过程结束。

经过迭代处理，原信号可以分解为多个分段正弦拟合函数和残差函数的和，即：

其中分段正弦拟合部分，根据定义4，按照公式(10)给出对应极值点相位，结合公式(13)，图2中第一段准周期的分段正弦函数可写为：

对应的瞬时相位为：

对应的瞬时频率为：

对应的瞬时幅度为：

如图4所示为风噪声分段正弦拟合分解的结果，停止参数残差能量与原信号能量之比小于0.01。如图5所示为1-3阶分段正弦拟合函数的瞬时频率随时间变化的曲线。用公式(17)可计算各阶的能量分布情况，如图6所示为随阶数变化的对数能量分布曲线，能量随阶数降低。

结合瞬时频率和瞬时幅度公式(14-15)，可给出信号的二维时频功率谱，如图7所示，，横纵坐标分别是时间和频率，灰度是对应频率的瞬时功率。

根据公式(18)用瞬时幅度对同一时刻的频率进行加权平均，可得加权平均瞬时频率曲线，根据公式19计算平均瞬时幅度，则平均时频分布可在二维时频图中画出，如图8所示，可观察主要信号频率随时间的变化情况。

Claims (2)

1.一种分段正弦拟合分解方法，其特征在于：包括以下步骤：

步骤一：初始化循环变量k、最大迭代次数Q和收敛阈值ε，令k＝1,r_k(t)＝x(t)，其中，x(t)为原信号，r_k(t)为迭代信号；

步骤二：寻找迭代信号r_k(t)的所有过零点，其中过零点对应时间序列满足：

T_Zero＝{tz_i|i∈[1,K],tz_i-1＜tz_i}

其中T_Zero代表迭代信号r_k(t)中过零点时间索引的集合,tz_i代表第i个过零点时间，K代表信号中过零点的个数；

过零点对应幅度序列满足：

V_Zero＝{az_i|az_i＝r_k(tz_i),i∈[1,K],tz_i∈T_zero}

其中r_k(tz_i)＝0代表过零点处的函数值为零，az_i代表第i个过零点的幅度；

步骤三：寻找相邻过零点之间的局部最大值点和局部最小值点，所述的局部最大值点是当相邻过零点之间的值都大于零时相邻过零点之间的最大值点，所述的局部最小值点是相邻过零点之间的值都小于零时相邻过零点之间的最小值点，所述的局部最大值点对应的时间索引序列满足：

T_Max＝{tx_i|i∈[1,L],tx_i＜tx_i+1}

其中T_Max代表信号r_k(t)中局部最大值点时间索引的集合,tx_i代表第i个局部最大值点对应的时间，L代表局部极大值点个数；

所述的局部最大值点对应的幅度序列满足：

V_Max＝{ax_i|ax_i＝r_k(tx_i),i∈[1,L],tx_i∈T_Max}

其中r_k(tx_i)代表信号局部最大值时刻tx_i对应的值，ax_i代表第i个局部最大值点的幅度；

所述的局部最小值点对应的时间索引序列满足：

T_Min＝{tn_i|i∈[1,M],tn_i＜tn_i+1}

其中T_Min代表信号r_k(t)中局部最小值点时间索引的集合,tn_i代表第i个局部最小值点对应的时间，M代表局部最小值点的个数；

所述的局部最小值点对应的幅度序列满足：

V_Min＝{an_i|an_i＝r_k(tn_i),i∈[1,M],tn_i∈T_Min}

其中r_k(tn_i)代表信号局部最小值时刻tn_i对应的值，an_i代表第i个局部最小值点的幅度；

过零点个数、局部极大值点个数和局部极小值点个数满足：

|K-(M+L)|＝1or 0

将过零点时间、局部最大值点时间和局部最小值点时间放在一起，并按时间先后排序，获得时间序列，即：

T_All＝{t_i|i∈[1,M+L+K],t_i＜t_i+1}＝{T_Zero,T_Min,T_Max}

按照时间序列T_All，可得到幅度序列：

V_All＝{a_i|a_i＝r_k(t_i),i∈[1,M+L+K],t_i∈T_All}＝{V_Zero,V_Min,V_Max}

其中t_i代表T_All中第i个时间，为过零点时间或者局部最大值点时间或者局部最小值点时间，t₁或t_M+L+K对应过零点时间，a_i是对应T_All中第i个时间的信号幅度，对应时间序列T_All的瞬时相位序列满足：

其中为对应时刻t_i的瞬时相位,如果a₂>0,则否则

步骤四：拟合中值曲线m(t)，从迭代信号r_k(t)中减去中值曲线m(t)，得到临时信号r₀(t)，临时信号r₀(t)满足：r₀(t)＝r_k(t)-m(t)，其中中值曲线m(t)拟合过程满足：

第一类点：用直线连接相邻的局部极大值点和局部极小值点，取其中点；第二类点：用直线连接相邻的两个局部极大值点，过所述相邻的两个局部极大值点之间的局部极小值点做平行于y轴的直线并与连接相邻的两个局部极大值点的直线相交，取交点和局部极小值点连线的中点；第三类点：用直线连接相邻的两个局部极小值点，过所述的相邻两个局部极小值点之间的局部极大值点做平行于y轴的直线并与连接相邻的两个局部极小值点的直线相交，取交点和局部极大值点连线的中点；以所述第一、二、三类点为节点进行插值，拟合信号的中值曲线m(t)；

步骤五：对所述的临时信号r₀(t)进行分段正弦拟合，获得拟合函数s_k(t)，所述的拟合函数s_k(t)满足：

其中a_i∈V_All,t_i∈T_All,1≤i≤K,且k∈N是自然数,K是过零点个数，如果信号的左右端点不为零，则：

瞬时相位函数满足：

瞬时频率是相位函数的一阶导数，表示为：

对应的瞬时幅度函数表示为：

a_k(t)＝{|a_2i|,t∈[t_2i-1,t_2i+1),i∈[1,K]}

步骤六：令r_k+1(t)＝r_k(t)-s_k(t)

步骤七：令k＝k+1；

步骤八：如果k>Q或条件成立，则算法迭代过程结束，否则重复执行步骤二至步骤八，迭代过程结束后，原信号满足：

其中：Q为分解阶数，与最大迭代次数相同；s_k(t)为第k阶分段正弦拟合函数，r_k+1(t)为分解后的残差函数。

2.根据权利要求1所述的一种分段正弦拟合分解方法，其特征在于：步骤四中所述插值方法包括线性插值法、多项式插值法、样条插值法。