CN108491563B - 一种基于稀疏重构最优化算法的信号包络线提取方法 - Google Patents

Abstract

本发明公开了一种基于稀疏重构最优化算法的信号包络线提取方法，其采用了基于粒子群稀疏重构最优化算法，将改变DCT基的频带宽度的变化因子作为基于粒子群稀疏重构最优化算法中的每个粒子的属性，通过多次迭代获得上包络线全局最优位置和下包络线全局最优位置，对应作为用于改变粒子确定的上包络线获取过程所需的DCT基的频带宽度的变化因子和用于改变粒子确定的下包络线获取过程所需的DCT基的频带宽度的变化因子，进而得到最适合上包络线变化趋势的最佳DCT基和最适合下包络线变化趋势的最佳DCT基，最终自适应地提取出最佳上包络线和最佳下包络线，提高了上包络线和下包络线的提取精度，同时提高了上包络线和下包络线提取的抗噪性能。

Description

一种基于稀疏重构最优化算法的信号包络线提取方法

技术领域

本发明涉及一种信号处理技术，尤其是涉及一种基于稀疏重构最优化算法的信号包络线提取方法。

背景技术

机械故障诊断及桥梁振动检测领域中，许多振动信号包含着调制信号成分，不同故障产生不同的调制形式。对调制信号进行分析处理的一个重要内容就是解调，解调分析能够获得调制信号的包络曲线和包络解调谱，而调制信号的包络往往集中地携带着有用的信息，因此，在信号处理和故障诊断等诸多领域中，包络解调方法的研究一直是众多学者关注的热点。

目前常用的包络分析方法有：广义检波滤波、Hilbert变换、三次样条插值等。其中，广义检波滤波在包络线提取过程中可能会产生混频效应，因此要求选择合适的采样频率对信号进行采样；基于Hilbert变换提取信号包络线的方法抗噪性比较弱，当信噪比变小时，Hilbert变换所得的包络误差会随之增大；基于三次样条插值方法提取包络线会产生端点效应问题，进而会影响经验模态分解的精度。

发明内容

本发明所要解决的技术问题是提供一种基于稀疏重构最优化算法的信号包络线提取方法，其不仅能够有效地抑制端点效应问题，而且能够有效地提高包络线提取的精度与抗噪性。

本发明解决上述技术问题所采用的技术方案为：一种基于稀疏重构最优化算法的信号包络线提取方法，其特征在于包括以下步骤：

步骤1)：将所要提取的包络线信号记为x(t)，将x(t)表示为列向量形式，x(t)＝[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]^T；其中，t表示采样时间，t的单位为秒，0≤t≤3，N表示采样点的总个数，N为正整数，N≥250，x₁(t)、x₂(t)、x_N-1(t)、x_N(t)对应表示x(t)中的第1个采样值、第2个采样值、第N-1个采样值、第N个采样值，符号“[]”为向量表示符号，[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]^T为[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]的转置；

步骤2)：从x(t)中提取出所有极大值点和所有极小值点；然后将从x(t)中提取出的所有极大值点构成的列向量记为y_d(t)，

并将从x(t)中提取出的所有极小值点构成的列向量记为y_s(t)，

其中，x_d,1、x_d,2和

对应表示从x(t)中提取出的第1个极大值点、第2个极大值点和第K₁个极大值点，K₁表示从x(t)中提取出的极大值点的总个数，K₁为正整数，1≤K₁＜N，

为

的转置，x_s,1、x_s,2和

对应表示从x(t)中提取出的第1个极小值点、第2个极小值点和第K₂个极小值点，K₂表示从x(t)中提取出的极小值点的总个数，K₂为正整数，1≤K₂＜N，

为

的转置；

步骤3)：基于粒子群稀疏重构最优化算法，提取x(t)的上包络线和下包络线，具体过程为：

3)_1、令α表示粒子群的迭代次数，令α_max表示粒子群的最大迭代次数；设粒子群中有L个粒子；其中，α为正整数，α的初始值为1，α_max为正整数，30≤α_max≤200，L为正整数，1＜L＜50；

3)_2、在第α次迭代过程中，将粒子群中当前待处理的粒子定义为当前粒子；

3)_3、设当前粒子为粒子群中的第i个粒子，其中，i为正整数，1≤i≤L；

3)_4、针对当前粒子确定一个上包络线获取过程所需的DCT基和一个下包络线获取过程所需的DCT基，对应记为Ψ_d,i和Ψ_s,i；并将用于改变Ψ_d,i的频带宽度的变化因子记为m_d,i，将用于改变Ψ_s,i的频带宽度的变化因子记为m_s,i；其中，Ψ_d,i和Ψ_s,i均为N阶方阵，0≤m_d,i≤1，0≤m_s,i≤1，当α＝1时m_d,i的值为上包络线获取过程中当前粒子的位置的初始值、m_s,i的值为下包络线获取过程中当前粒子的位置的初始值，上包络线获取过程中当前粒子的位置的初始值和下包络线获取过程中当前粒子的位置的初始值均为当前粒子的位置的初始值；

3)_5、将Ψ_d,i中行号与x(t)中的所有极大值点的索引值对应的行进行标记，然后从Ψ_d,i中提取出标记的K₁行构建成一个新的维数为K₁×N的矩阵，记为H_d,i；并将Ψ_s,i中行号与x(t)中的所有极小值点的索引值对应的行进行标记，然后从Ψ_s,i中提取出标记的K₂行构建成一个新的维数为K₂×N的矩阵，记为H_s,i；

3)_6、将y_d(t)作为当前粒子在稀疏重构过程中的第一观测值，并将H_d,i作为当前粒子在稀疏重构过程中的第一感知矩阵；然后根据y_d(t)和H_d,i，并利用内点法，求解当前粒子在稀疏重构过程中的第一稀疏解，记为η_d,i；其中，η_d,i的维数为N×1；

同理，将y_s(t)作为当前粒子在稀疏重构过程中的第二观测值，并将H_s,i作为当前粒子在稀疏重构过程中的第二感知矩阵；然后根据y_s(t)和H_s,i，并利用内点法，求解当前粒子在稀疏重构过程中的第二稀疏解，记为η_s,i；其中，η_s,i的维数为N×1；

3)_7、计算当前粒子所对应的x(t)的上包络线，记为

并计算当前粒子所对应的x(t)的下包络线，记为

3)_8、构建一个维数为(N-1)×N的一阶微分矩阵，记为D，

接着根据D和

确定当前粒子所对应的用于体现

的变化程度的上包络线平稳度函数，记为E_d,i(t)，

并根据D和

确定当前粒子所对应的用于体现

的变化程度的下包络线平稳度函数，记为E_s,i(t)，

其中，符号“|| ||₂”为求矩阵的2-范数符号；

3)_9、计算当前粒子的上包络线适应度函数，记为f_d,i，

并计算当前粒子的下包络线适应度函数，记为

其中，λ为正则化参数，λ＞0，符号“|| ||₁”为求矩阵的1-范数符号；

3)_10、在第α次迭代过程中，将粒子群中下一个待处理的粒子作为当前粒子，然后返回步骤3)_3继续执行，直至粒子群中的所有粒子处理完毕；

3)_11、在第α次迭代过程中，根据L个上包络线适应度函数，从L个粒子中找出使上包络线平稳度函数最小的粒子，将找出的粒子的位置作为粒子群中的所有粒子的上包络线全局最优位置，记为g_d,best；并针对粒子群中的每个粒子，从第1次迭代过程至第α次迭代过程中找出使上包络线平稳度函数最小的迭代序号，将找出的迭代序号对应的迭代过程中该粒子的位置作为该粒子的上包络线局部最优位置，将粒子群中的第i个粒子的上包络线局部最优位置记为p_d,i,best；同样，根据L个下包络线适应度函数，从L个粒子中找出使下包络线平稳度函数最小的粒子，将找出的粒子的位置作为粒子群中的所有粒子的下包络线全局最优位置，记为g_s,best；并针对粒子群中的每个粒子，从第1次迭代过程至第α次迭代过程中找出使下包络线平稳度函数最小的迭代序号，将找出的迭代序号对应的迭代过程中该粒子的位置作为该粒子的下包络线局部最优位置，将粒子群中的第i个粒子的下包络线局部最优位置记为p_s,i,best；

3)_12、判断α＜α_max是否成立，如果α＜α_max成立，则更新上包络线获取过程中粒子群中的每个粒子的速度和位置、下包络线获取过程中粒子群中的每个粒子的速度和位置，将上包络线获取过程中粒子群中的第i个粒子更新后的速度记为v_d,i,back，v_d,i,back＝w×v_d,i,pre+c₁×r₁×(p_d,i,best,pre-p_d,i,pre)+c₂×r₂×(g_d,best,pre-p_d,i,best,pre)，将上包络线获取过程中粒子群中的第i个粒子更新后的位置记为p_d,i,back，p_d,i,back＝p_d,i,pre+v_d,i,pre，将下包络线获取过程中粒子群中的第i个粒子更新后的速度记为v_s,i,back，v_s,i,back＝w×v_s,i,pre+c₁×r₁×(p_s,i,best,pre-p_s,i,pre)+c₂×r₂×(g_s,best,pre-p_s,i,best,pre)，将下包络线获取过程中粒子群中的第i个粒子更新后的位置记为p_s,i,back，p_s,i,back＝p_s,i,pre+v_s,i,pre；然后将上包络线获取过程中粒子群中的每个粒子更新后的速度和位置对应作为下一次迭代过程中上包络线获取过程中该粒子的速度和位置，将下包络线获取过程中粒子群中的每个粒子更新后的速度和位置对应作为下一次迭代过程中下包络线获取过程中该粒子的速度和位置，接着令α＝α+1；再返回步骤3)_2继续执行；其中，w为惯性因子，w的取值范围为[0,1]，v_d,i,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的速度，c₁和c₂均为学习因子，r₁和r₂均为区间[0,1]内的随机数，p_d,i,best,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的上包络线局部最优位置，p_d,i,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的位置，g_d,best,pre表示更新前粒子群中的所有粒子的上包络线全局最优位置，v_s,i,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的速度，p_s,i,best,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的下包络线局部最优位置，p_s,i,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的位置，g_s,best,pre表示更新前粒子群中的所有粒子的下包络线全局最优位置，α＝α+1中的“＝”为赋值符号；

如果α＜α_max不成立，则结束迭代过程，然后将第α次迭代过程中获得的g_d,best赋值给用于改变一个粒子确定的上包络线获取过程所需的DCT基的频带宽度的变化因子，将第α次迭代过程中获得的g_s,best赋值给用于改变一个粒子确定的下包络线获取过程所需的DCT基的频带宽度的变化因子，再按照步骤3)_4至步骤3)_7的过程，以相同的方式执行一次，获得x(t)的最终上包络线即为最佳上包络线和x(t)的最终下包络线即为最佳下包络线。

所述的步骤2)中的y_d(t)的获取过程为：

步骤2_a1)：对x(t)求一阶导数，将得到的列向量记为d_x1；然后将d_x1中大于零的所有元素的索引值组成的位置向量记为a；

步骤2_a2)：对x(t)求二阶导数，将得到的列向量记为d_x2；然后将d_x2中小于零的所有元素的索引值组成的位置向量记为b；

步骤2_a3)：将a和b中相同的元素提取出来，且只选用从a中提取出的所有元素或只选用从b中提取出的所有元素；然后将选用的所有元素组成列向量，记为γ；再将γ中的每个元素加1，得到新的列向量，记为γ'，γ'中的每个元素即为x(t)中的一个极大值点的索引值；

步骤2_a4)：从x(t)中提取出索引值与γ'中的每个元素相同的采样值；然后将提取出的这些采样值按序排列组成y_d(t)；

所述的步骤2)中的y_s(t)的获取过程为：

步骤2_b1)：取x(t)中的每个采样值的相反数，将得到的新的列向量记为x'(t)；

步骤2_b2)：对x'(t)求一阶导数，将得到的列向量记为d'_x1；然后将d'_x1中大于零的所有元素的索引值组成的位置向量记为a'；

步骤2_b3)：对x'(t)求二阶导数，将得到的列向量记为d'_x2；然后将d'_x2中小于零的所有元素的索引值组成的位置向量记为b'；

步骤2_b4)：将a'和b'中相同的元素提取出来，且只选用从a'中提取出的所有元素或只选用从b'中提取出的所有元素；然后将选用的所有元素组成列向量，记为γ*；接着将γ*中的每个元素加1，得到新的列向量，记为γ*'；再取γ*'中的每个元素的相反数，将得到的新的列向量记为γ'*'，γ'*'中的每个元素即为x(t)中的一个极小值点的索引值；

步骤2_b5)：从x(t)中提取出索引值与γ'*'中的每个元素相同的采样值；然后将提取出的这些采样值按序排列组成y_s(t)。

所述的步骤3)_4中的Ψ_d,i的获取过程为：将Ψ_d,i中下标为(k,n)的元素记为

其中，k和n均为正整数，1≤k≤N，1≤n≤N，cos()为求余弦函数；

所述的步骤3)_4中的Ψ_s,i的获取过程为：将Ψ_s,i中下标为(k,n)的元素记为

所述的步骤3)_6中η_d,i的求解过程为：

3)_6_a1、令k1表示求解迭代次数，令η_d,i的初始值为[0,0,…,0]^T _N×1；其中，k1为正整数，k1的初始值为1；

3)_6_b1、构造第k1次求解迭代过程中的惩罚函数Α(η_d,i,u_d,i)，Α(η_d,i,u_d,i)的表达式为：

并将Α(η_d,i,u_d,i)的惩罚因子记为τ，将Α(η_d,i,u_d,i)的衰减系数记为ρ；其中，q为正整数，1≤q≤N，u_d,i(q)表示所有元素均大于零且维数为N×1的列向量u_d,i中的第q个元素，且满足-u_d,i(q)≤η_d,i(q)≤+u_d,i(q)，u_d,i(q)的初始值为1，η_d,i(q)表示η_d,i中的第q个元素，τ的初始值大于0，ρ的值为区间(0,1)之间的数；

3)_6_c1、在Α(η_d,i,u_d,i)基础上构造第k1次求解迭代过程中的无约束子问题T(η_d,i,u_d,i)，T(η_d,i,u_d,i)的表达式为：

其中，符号“|| ||₂”为求矩阵的2-范数符号，λ为正则化参数，λ＞0；

3)_6_d1、利用共轭梯度法求解T(η_d,i,u_d,i)，得到η_d,i在第k1次求解迭代后的解值；

3)_6_e1、在第k1次求解迭代过程中，判断τ×Α(η_d,i,u_d,i)≤ε是否满足，如果满足，则停止求解迭代过程，并将η_d,i在第k1次求解迭代后的解值作为η_d,i的最终值；如果不满足，则将η_d,i在第k1次求解迭代后的解值作为下一次求解迭代过程中的值，然后令k1＝k1+1，并令τ＝ρ×τ，再返回步骤3)_6_b1继续执行；其中，k1＝k1+1和τ＝ρ×τ中的“＝”为赋值符号，ε为误差终止值，0≤ε≤1；

所述的步骤3)_6中η_s,i的求解过程为：

3)_6_a2、令k2表示求解迭代次数，令η_s,i的初始值为[0,0,…,0]^T _N×1；其中，k2为正整数，k2的初始值为1；

3)_6_b2、构造第k2次求解迭代过程中的惩罚函数Α(η_s,i,u_s,i)，Α(η_s,i,u_s,i)的表达式为：

并将Α(η_s,i,u_s,i)的惩罚因子记为τ'，将Α(η_s,i,u_s,i)的衰减系数记为ρ'；其中，q为正整数，1≤q≤N，u_s,i(q)表示所有元素均大于零且维数为N×1的列向量u_s,i中的第q个元素，且满足-u_s,i(q)≤η_s,i(q)≤+u_s,i(q)，u_s,i(q)的初始值为1，η_s,i(q)表示η_s,i中的第q个元素，τ'的初始值大于0，ρ'的值为区间(0,1)之间的数；

3)_6_c2、在Α(η_s,i,u_s,i)基础上构造第k2次求解迭代过程中的无约束子问题T(η_s,i,u_s,i)，T(η_s,i,u_s,i)的表达式为：

其中，符号“|| ||₂”为求矩阵的2-范数符号，λ为正则化参数，λ＞0；

3)_6_d2、利用共轭梯度法求解T(η_s,i,u_s,i)，得到η_s,i在第k2次求解迭代后的解值；

3)_6_e2、在第k2次求解迭代过程中，判断τ'×Α(η_s,i,u_s,i)≤ε是否满足，如果满足，则停止求解迭代过程，并将η_s,i在第k2次求解迭代后的解值作为η_s,i的最终值；如果不满足，则将η_s,i在第k2次求解迭代后的解值作为下一次求解迭代过程中的值，然后令k2＝k2+1，并令τ'＝ρ'×τ'，再返回步骤3)_6_b2继续执行；其中，k2＝k2+1和τ'＝ρ'×τ'中的“＝”为赋值符号，ε为误差终止值，0≤ε≤1。

与现有技术相比，本发明的优点在于：

1)本发明方法采用了基于粒子群稀疏重构最优化算法，将改变DCT基的频带宽度的变化因子作为基于粒子群稀疏重构最优化算法中的每一个粒子的属性，通过多次迭代获得上包络线全局最优位置和下包络线全局最优位置，对应作为用于改变每个粒子确定的上包络线获取过程所需的DCT基的频带宽度的变化因子和用于改变每个粒子确定的下包络线获取过程所需的DCT基的频带宽度的变化因子，进而得到最适合上包络线变化趋势的最佳DCT基和最适合下包络线变化趋势的最佳DCT基，最终自适应地提取出最佳上包络线和最佳下包络线，本发明方法通过基于粒子群稀疏重构最优化算法提高了上包络线和下包络线的提取精度，同时提高了上包络线和下包络线提取的抗噪性能。

2)本发明方法中的稀疏重构过程采用的是内点法，以极值点组成的向量作为稀疏复原中的观测值，然后寻找到最优DCT基复原出包络线信号，内点法的引入将包络线的提取过程变得简单易操作，并且解决了端点效应问题，同时也提高了包络线提取的抗噪性能。

附图说明

图1为本发明方法的总体实现框图；

图2为所要提取的上包络线信号、采用本发明方法自适应地从未含噪信号中提取出的最佳上包络线、采用现有的三次样条插值法从未含噪信号中提取出的上包络线的对比示意图；

图3为所要提取的上包络线信号、采用本发明方法自适应地从信噪比为15dB含噪信号中提取出的最佳上包络线、采用现有的三次样条插值法从信噪比为15dB含噪信号中提取出的上包络线的对比示意图；

图4为所要提取的上包络线信号、采用本发明方法自适应地从信噪比为10dB含噪信号中提取出的最佳上包络线、采用现有的三次样条插值法从信噪比为10dB含噪信号中提取出的上包络线的对比示意图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种基于稀疏重构最优化算法的信号包络线提取方法，其流程框图如图1所示，其包括以下步骤：

步骤1)：将所要提取的包络线信号记为x(t)，将x(t)表示为列向量形式，x(t)＝[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]^T；其中，t表示采样时间，t的单位为秒，0≤t≤3，如取t＝1.5，N表示采样点的总个数，N为正整数，N≥250，一般N的取值要适当，N的取值过小的话会影响包络线的复原，N的取值过大的话会增加计算复杂度，因此在本实施例中仿真信号的N取值为500，x₁(t)、x₂(t)、x_N-1(t)、x_N(t)对应表示x(t)中的第1个采样值、第2个采样值、第N-1个采样值、第N个采样值，符号“[]”为向量表示符号，[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]^T为[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]的转置。

步骤2)：从x(t)中提取出所有极大值点和所有极小值点；然后将从x(t)中提取出的所有极大值点构成的列向量记为y_d(t)，

并将从x(t)中提取出的所有极小值点构成的列向量记为y_s(t)，

其中，x_d,1、x_d,2和

对应表示从x(t)中提取出的第1个极大值点、第2个极大值点和第K₁个极大值点，K₁表示从x(t)中提取出的极大值点的总个数，K₁为正整数，1≤K₁＜N，

为

的转置，x_s,1、x_s,2和

对应表示从x(t)中提取出的第1个极小值点、第2个极小值点和第K₂个极小值点，K₂表示从x(t)中提取出的极小值点的总个数，K₂为正整数，1≤K₂＜N，

为

的转置。

在本实施例中，步骤2)中的y_d(t)的获取过程为：

步骤2_a1)：对x(t)求一阶导数，将得到的列向量记为d_x1；然后将d_x1中大于零的所有元素的索引值组成的位置向量记为a。

步骤2_a2)：对x(t)求二阶导数，将得到的列向量记为d_x2；然后将d_x2中小于零的所有元素的索引值组成的位置向量记为b。

步骤2_a3)：将a和b中相同的元素提取出来，且只选用从a中提取出的所有元素或只选用从b中提取出的所有元素；然后将选用的所有元素组成列向量，记为γ；再将γ中的每个元素加1，得到新的列向量，记为γ'，γ'中的每个元素即为x(t)中的一个极大值点的索引值。

步骤2_a4)：从x(t)中提取出索引值与γ'中的每个元素相同的采样值；然后将提取出的这些采样值按序排列组成y_d(t)。

在本实施例中，步骤2)中的y_s(t)的获取过程为：

步骤2_b1)：取x(t)中的每个采样值的相反数，将得到的新的列向量记为x'(t)。

步骤2_b2)：对x'(t)求一阶导数，将得到的列向量记为d'_x1；然后将d'_x1中大于零的所有元素的索引值组成的位置向量记为a'。

步骤2_b3)：对x'(t)求二阶导数，将得到的列向量记为d'_x2；然后将d'_x2中小于零的所有元素的索引值组成的位置向量记为b'。

步骤2_b4)：将a'和b'中相同的元素提取出来，且只选用从a'中提取出的所有元素或只选用从b'中提取出的所有元素；然后将选用的所有元素组成列向量，记为γ*；接着将γ*中的每个元素加1，得到新的列向量，记为γ*'；再取γ*'中的每个元素的相反数，将得到的新的列向量记为γ'*'，γ'*'中的每个元素即为x(t)中的一个极小值点的索引值。

步骤2_b5)：从x(t)中提取出索引值与γ'*'中的每个元素相同的采样值；然后将提取出的这些采样值按序排列组成y_s(t)。

步骤3)：基于粒子群稀疏重构最优化算法，提取x(t)的上包络线和下包络线，具体过程为：

3)_1、令α表示粒子群的迭代次数，令α_max表示粒子群的最大迭代次数；设粒子群中有L个粒子；其中，α为正整数，α的初始值为1，α_max为正整数，30≤α_max≤200，L为正整数，1＜L＜50，一般情况下L的取值不宜太大，太大会降低搜索速度，因此在本实施例中取L＝20。

3)_2、在第α次迭代过程中，将粒子群中当前待处理的粒子定义为当前粒子。

3)_3、设当前粒子为粒子群中的第i个粒子，其中，i为正整数，1≤i≤L。

3)_4、针对当前粒子确定一个上包络线获取过程所需的DCT基和一个下包络线获取过程所需的DCT基，对应记为Ψ_d,i和Ψ_s,i；并将用于改变Ψ_d,i的频带宽度的变化因子记为m_d,i，将用于改变Ψ_s,i的频带宽度的变化因子记为m_s,i；其中，Ψ_d,i和Ψ_s,i均为N阶方阵，0≤m_d,i≤1，0≤m_s,i≤1，当α＝1时m_d,i的值为上包络线获取过程中当前粒子的位置的初始值、m_s,i的值为下包络线获取过程中当前粒子的位置的初始值，上包络线获取过程中当前粒子的位置的初始值和下包络线获取过程中当前粒子的位置的初始值均为当前粒子的位置的初始值。

在本实施例中，步骤3)_4中的Ψ_d,i的获取过程为：将Ψ_d,i中下标为(k,n)的元素记为

其中，k和n均为正整数，1≤k≤N，1≤n≤N，cos()为求余弦函数；

在本实施例中，步骤3)_4中的Ψ_s,i的获取过程为：将Ψ_s,i中下标为(k,n)的元素记为

3)_5、将Ψ_d,i中行号与x(t)中的所有极大值点的索引值对应的行进行标记，然后从Ψ_d,i中提取出标记的K₁行构建成一个新的维数为K₁×N的矩阵，记为H_d,i；并将Ψ_s,i中行号与x(t)中的所有极小值点的索引值对应的行进行标记，然后从Ψ_s,i中提取出标记的K₂行构建成一个新的维数为K₂×N的矩阵，记为H_s,i。

3)_6、将y_d(t)作为当前粒子在稀疏重构过程中的第一观测值，并将H_d,i作为当前粒子在稀疏重构过程中的第一感知矩阵；然后根据y_d(t)和H_d,i，并利用内点法，求解当前粒子在稀疏重构过程中的第一稀疏解，记为η_d,i；其中，η_d,i的维数为N×1。

在本实施例中，步骤3)_6中η_d,i的求解过程为：

3)_6_a1、令k1表示求解迭代次数，令η_d,i的初始值为[0,0,...,0]^T _N×1；其中，k1为正整数，k1的初始值为1。

3)_6_b1、构造第k1次求解迭代过程中的惩罚函数Α(η_d,i,u_d,i)，Α(η_d,i,u_d,i)的表达式为：

并将Α(η_d,i,u_d,i)的惩罚因子记为τ，将Α(η_d,i,u_d,i)的衰减系数记为ρ；其中，q为正整数，1≤q≤N，u_d,i(q)表示所有元素均大于零且维数为N×1的列向量u_d,i中的第q个元素，且满足-u_d,i(q)≤η_d,i(q)≤+u_d,i(q)，u_d,i(q)的初始值为1，即u_d,i的初始值中的所有元素均为1，u_d,i的初始值为[1,1,...,1]^T _N×1，η_d,i(q)表示η_d,i中的第q个元素，τ的初始值大于0，在本实施例中取τ的初始值为1，ρ的值为区间(0,1)之间的数，在本实施例中取ρ的值为0.5。

3)_6_c1、在Α(η_d,i,u_d,i)基础上构造第k1次求解迭代过程中的无约束子问题T(η_d,i,u_d,i)，T(η_d,i,u_d,i)的表达式为：

其中，符号“|| ||₂”为求矩阵的2-范数符号，λ为正则化参数，λ＞0，在本实施例中取λ＝0.1。

3)_6_d1、利用共轭梯度法求解T(η_d,i,u_d,i)，得到η_d,i在第k1次求解迭代后的解值。

3)_6_e1、在第k1次求解迭代过程中，判断τ×Α(η_d,i,u_d,i)≤ε是否满足，如果满足，则停止求解迭代过程，并将η_d,i在第k1次求解迭代后的解值作为η_d,i的最终值；如果不满足，则将η_d,i在第k1次求解迭代后的解值作为下一次求解迭代过程中的值，然后令k1＝k1+1，并令τ＝ρ×τ，再返回步骤3)_6_b1继续执行；其中，k1＝k1+1和τ＝ρ×τ中的“＝”为赋值符号，ε为误差终止值，0≤ε≤1，在本实施例中取ε＝10^-6。

同理，将y_s(t)作为当前粒子在稀疏重构过程中的第二观测值，并将H_s,i作为当前粒子在稀疏重构过程中的第二感知矩阵；然后根据y_s(t)和H_s,i，并利用内点法，求解当前粒子在稀疏重构过程中的第二稀疏解，记为η_s,i；其中，η_s,i的维数为N×1。

在本实施例中，步骤3)_6中η_s,i的求解过程为：

3)_6_a2、令k2表示求解迭代次数，令η_s,i的初始值为[0,0,...,0]^T _N×1；其中，k2为正整数，k2的初始值为1。

3)_6_b2、构造第k2次求解迭代过程中的惩罚函数Α(η_s,i,u_s,i)，Α(η_s,i,u_s,i)的表达式为：

并将Α(η_s,i,u_s,i)的惩罚因子记为τ'，将Α(η_s,i,u_s,i)的衰减系数记为ρ'；其中，q为正整数，1≤q≤N，u_s,i(q)表示所有元素均大于零且维数为N×1的列向量u_s,i中的第q个元素，且满足-u_s,i(q)≤η_s,i(q)≤+u_s,i(q)，u_s,i(q)的初始值为1，即u_s,i的初始值中的所有元素均为1，u_s,i的初始值为[1,1,...,1]^T _N×1，η_s,i(q)表示η_s,i中的第q个元素，τ'的初始值大于0，在本实施例中取τ'的初始值为1，ρ'的值为区间(0,1)之间的数，在本实施例中取ρ'的值为0.5。

3)_6_c2、在Α(η_s,i,u_s,i)基础上构造第k2次求解迭代过程中的无约束子问题T(η_s,i,u_s,i)，T(η_s,i,u_s,i)的表达式为：

其中，符号“|| ||₂”为求矩阵的2-范数符号，λ为正则化参数，λ＞0，在本实施例中取λ＝0.1。

3)_6_d2、利用共轭梯度法求解T(η_s,i,u_s,i)，得到η_s,i在第k2次求解迭代后的解值。

3)_6_e2、在第k2次求解迭代过程中，判断τ'×Α(η_s,i,u_s,i)≤ε是否满足，如果满足，则停止求解迭代过程，并将η_s,i在第k2次求解迭代后的解值作为η_s,i的最终值；如果不满足，则将η_s,i在第k2次求解迭代后的解值作为下一次求解迭代过程中的值，然后令k2＝k2+1，并令τ'＝ρ'×τ'，再返回步骤3)_6_b2继续执行；其中，k2＝k2+1和τ'＝ρ'×τ'中的“＝”为赋值符号，ε为误差终止值，0≤ε≤1，在本实施例中取ε＝10^-6。

3)_7、计算当前粒子所对应的x(t)的上包络线，记为

并计算当前粒子所对应的x(t)的下包络线，记为

3)_8、构建一个维数为(N-1)×N的一阶微分矩阵，记为D，

接着根据D和

确定当前粒子所对应的用于体现

的变化程度的上包络线平稳度函数，记为E_d,i(t)，

并根据D和

确定当前粒子所对应的用于体现

的变化程度的下包络线平稳度函数，记为E_s,i(t)，

其中，D中第ω行第ω列的元素的值为1，第ω行第ω+1列的元素的值为-1，其余行和列的所有元素均为0，符号“|| ||₂”为求矩阵的2-范数符号。

3)_9、计算当前粒子的上包络线适应度函数，记为f_d,i，

式中

保证了上包络线通过所有的极大值点，式中λ||η_d,i||₁确保了所求的η_d,i是最稀疏的，E_d,i(t)用于使得所求上包络线最光滑；并计算当前粒子的下包络线适应度函数，记为f_s,i，

式中

保证了下包络线通过所有的极小值点，式中λ||η_s,i||₁确保了所求的η_s,i是最稀疏的，E_s,i(t)用于使得所求下包络线最光滑；其中，λ为正则化参数，λ＞0，在本实施例中取λ＝0.1，符号“|| ||₁”为求矩阵的1-范数符号。

3)_10、在第α次迭代过程中，将粒子群中下一个待处理的粒子作为当前粒子，然后返回步骤3)_3继续执行，直至粒子群中的所有粒子处理完毕。

3)_11、在第α次迭代过程中，根据L个上包络线适应度函数，从L个粒子中找出使上包络线平稳度函数最小的粒子，将找出的粒子的位置作为粒子群中的所有粒子的上包络线全局最优位置，记为g_d,best；并针对粒子群中的每个粒子，从第1次迭代过程至第α次迭代过程中找出使上包络线平稳度函数最小的迭代序号，将找出的迭代序号对应的迭代过程中该粒子的位置作为该粒子的上包络线局部最优位置，将粒子群中的第i个粒子的上包络线局部最优位置记为p_d,i,best；同样，根据L个下包络线适应度函数，从L个粒子中找出使下包络线平稳度函数最小的粒子，将找出的粒子的位置作为粒子群中的所有粒子的下包络线全局最优位置，记为g_s,best；并针对粒子群中的每个粒子，从第1次迭代过程至第α次迭代过程中找出使下包络线平稳度函数最小的迭代序号，将找出的迭代序号对应的迭代过程中该粒子的位置作为该粒子的下包络线局部最优位置，将粒子群中的第i个粒子的下包络线局部最优位置记为p_s,i,best。

3)_12、判断α＜α_max是否成立，如果α＜α_max成立，则更新上包络线获取过程中粒子群中的每个粒子的速度和位置、下包络线获取过程中粒子群中的每个粒子的速度和位置，将上包络线获取过程中粒子群中的第i个粒子更新后的速度记为v_d,i,back，v_d,i,back＝w×v_d,i,pre+c₁×r₁×(p_d,i,best,pre-p_d,i,pre)+c₂×r₂×(g_d,best,pre-p_d,i,best,pre)，将上包络线获取过程中粒子群中的第i个粒子更新后的位置记为p_d,i,back，p_d,i,back＝p_d,i,pre+v_d,i,pre，将下包络线获取过程中粒子群中的第i个粒子更新后的速度记为v_s,i,back，v_s,i,back＝w×v_s,i,pre+c₁×r₁×(p_s,i,best,pre-p_s,i,pre)+c₂×r₂×(g_s,best,pre-p_s,i,best,pre)，将下包络线获取过程中粒子群中的第i个粒子更新后的位置记为p_s,i,back，p_s,i,back＝p_s,i,pre+v_s,i,pre；然后将上包络线获取过程中粒子群中的每个粒子更新后的速度和位置对应作为下一次迭代过程中上包络线获取过程中该粒子的速度和位置，将下包络线获取过程中粒子群中的每个粒子更新后的速度和位置对应作为下一次迭代过程中下包络线获取过程中该粒子的速度和位置，接着令α＝α+1；再返回步骤3)_2继续执行；其中，w为惯性因子，w的取值范围为[0,1]，在本实施例中取w＝0.5，v_d,i,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的速度，c₁和c₂均为学习因子，在本实施例中取c₁＝1和c₂＝2，r₁和r₂均为区间[0,1]内的随机数，p_d,i,best,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的上包络线局部最优位置，p_d,i,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的位置，g_d,best,pre表示更新前粒子群中的所有粒子的上包络线全局最优位置，v_s,i,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的速度，p_s,i,best,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的下包络线局部最优位置，p_s,i,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的位置，g_s,best,pre表示更新前粒子群中的所有粒子的下包络线全局最优位置，α＝α+1中的“＝”为赋值符号。

如果α＜α_max不成立，则结束迭代过程，然后将第α次迭代过程中获得的g_d,best赋值给用于改变一个粒子确定的上包络线获取过程所需的DCT基的频带宽度的变化因子，将第α次迭代过程中获得的g_s,best赋值给用于改变一个粒子确定的下包络线获取过程所需的DCT基的频带宽度的变化因子，再按照步骤3)_4至步骤3)_7的过程，以相同的方式执行一次，获得x(t)的最终上包络线即为最佳上包络线和x(t)的最终下包络线即为最佳下包络线。

为了进一步说明本发明方法的可行性和有效性，对本发明方法进行仿真实验。

在本实施例中仿真信号为：x(t)＝1.8e^-3tcos(2πf₁t+1.3)+0.6e^-3tcos(2πf₂t+0.7)+0.1e^-2tcos(2πf₃t+0.7)，f₁、f₂和f₃均为模态频率，分别取f₁为0.8Hz、f₂为3Hz、f₃为10Hz，采样时间t的取值范围为[0,3]，单位为秒，采样点数N为500个。

为了突出本发明方法的优越性，结合基于粒子群稀疏重构最优化算法对上面的仿真信号采用内点法复原上包络线，同时与三次样条差值算法对比。图2给出了所要提取的上包络线信号(原始信号)、采用本发明方法自适应地从未含噪信号中提取出的最佳上包络线(粒子群稀疏优化)、采用现有的三次样条插值法从未含噪信号中提取出的上包络线(三次样条插值)的对比示意图。图2中dd值为定义的用于体现上包络线变化程度的平稳度函数值E_d(t)，m为改变DCT基的频带宽度的变化因子。从图2中可以看出，基于粒子群稀疏重构最优化算法复原上包络线明显地抑制了端点效应问题，而且自适应地找到了改变DCT基的频带宽度的变化因子，从而自适应地寻找到一条符合上包络线缓慢变化特性的最佳曲线，其中上包络线的变化程度可以从dd值大小看出由于上包络线具有变化缓慢的特性，因此dd值越小越好，本发明方法使得dd＝0.02。

为了检验本发明方法对信号处理的抗噪性能，仍然采用上述仿真信号x(t)，然后给x(t)添加不同信噪比的高斯白噪声n(t)，得到加噪信号x₁(t)，然后对x₁(t)用本发明方法提取出上包络线。图3给出了所要提取的上包络线信号、采用本发明方法自适应地从信噪比为15dB含噪信号中提取出的最佳上包络线、采用现有的三次样条插值法从信噪比为15dB含噪信号中提取出的上包络线的对比示意图；图4给出了所要提取的上包络线信号、采用本发明方法自适应地从信噪比为10dB含噪信号中提取出的最佳上包络线、采用现有的三次样条插值法从信噪比为10dB含噪信号中提取出的上包络线的对比示意图。图3和图4中dd值为定义的用于体现上包络线变化程度的平稳度函数值E_d(t)，m为改变DCT基的频带宽度的变化因子，图3中添加的高斯白噪声n(t)的信噪比为15dB，图4中添加的高斯白噪声n(t)的信噪比为10dB，信噪比的数值越小，表明噪声越大。从图3和图4所示的仿真结果中可以看出，在含有噪声的情况下本发明方法仍可以精确地提取出上包络信号，同样可以解决端点效应问题，并且随着噪声强度的加大，提取上包络线的精度没有明显降低，其误差仍可以控制在一定的范围之内，从而说明了本发明方法的有效性与实用性。

Claims (4)

1.一种基于稀疏重构最优化算法的信号包络线提取方法，其特征在于包括以下步骤：

步骤1)：将所要提取的包络线信号记为x(t)，将x(t)表示为列向量形式，x(t)＝[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]^T；其中，t表示采样时间，t的单位为秒，0≤t≤3，N表示采样点的总个数，N为正整数，N≥250，x₁(t)、x₂(t)、x_N-1(t)、x_N(t)对应表示x(t)中的第1个采样值、第2个采样值、第N-1个采样值、第N个采样值，符号“[]”为向量表示符号，[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]^T为[x₁(t),x₂(t),...,x_N-1(t),x_N(t)]的转置；

步骤2)：从x(t)中提取出所有极大值点和所有极小值点；然后将从x(t)中提取出的所有极大值点构成的列向量记为y_d(t)，

并将从x(t)中提取出的所有极小值点构成的列向量记为y_s(t)，

其中，x_d,1、x_d,2和

对应表示从x(t)中提取出的第1个极大值点、第2个极大值点和第K₁个极大值点，K₁表示从x(t)中提取出的极大值点的总个数，K₁为正整数，1≤K₁＜N，

为

的转置，x_s,1、x_s,2和

对应表示从x(t)中提取出的第1个极小值点、第2个极小值点和第K₂个极小值点，K₂表示从x(t)中提取出的极小值点的总个数，K₂为正整数，1≤K₂＜N，

为

的转置；

步骤3)：基于粒子群稀疏重构最优化算法，提取x(t)的上包络线和下包络线，具体过程为：

3)_1、令α表示粒子群的迭代次数，令α_max表示粒子群的最大迭代次数；设粒子群中有L个粒子；其中，α为正整数，α的初始值为1，α_max为正整数，30≤α_max≤200，L为正整数，1＜L＜50；

3)_2、在第α次迭代过程中，将粒子群中当前待处理的粒子定义为当前粒子；

3)_3、设当前粒子为粒子群中的第i个粒子，其中，i为正整数，1≤i≤L；

3)_4、针对当前粒子确定一个上包络线获取过程所需的DCT基和一个下包络线获取过程所需的DCT基，对应记为Ψ_d,i和Ψ_s,i；并将用于改变Ψ_d,i的频带宽度的变化因子记为m_d,i，将用于改变Ψ_s,i的频带宽度的变化因子记为m_s,i；其中，Ψ_d,i和Ψ_s,i均为N阶方阵，0≤m_d,i≤1，0≤m_s,i≤1，当α＝1时m_d,i的值为上包络线获取过程中当前粒子的位置的初始值、m_s,i的值为下包络线获取过程中当前粒子的位置的初始值，上包络线获取过程中当前粒子的位置的初始值和下包络线获取过程中当前粒子的位置的初始值均为当前粒子的位置的初始值；

3)_5、将Ψ_d,i中行号与x(t)中的所有极大值点的索引值对应的行进行标记，然后从Ψ_d,i中提取出标记的K₁行构建成一个新的维数为K₁×N的矩阵，记为H_d,i；并将Ψ_s,i中行号与x(t)中的所有极小值点的索引值对应的行进行标记，然后从Ψ_s,i中提取出标记的K₂行构建成一个新的维数为K₂×N的矩阵，记为H_s,i；

3)_6、将y_d(t)作为当前粒子在稀疏重构过程中的第一观测值，并将H_d,i作为当前粒子在稀疏重构过程中的第一感知矩阵；然后根据y_d(t)和H_d,i，并利用内点法，求解当前粒子在稀疏重构过程中的第一稀疏解，记为η_d,i；其中，η_d,i的维数为N×1；

同理，将y_s(t)作为当前粒子在稀疏重构过程中的第二观测值，并将H_s,i作为当前粒子在稀疏重构过程中的第二感知矩阵；然后根据y_s(t)和H_s,i，并利用内点法，求解当前粒子在稀疏重构过程中的第二稀疏解，记为η_s,i；其中，η_s,i的维数为N×1；

3)_7、计算当前粒子所对应的x(t)的上包络线，记为

并计算当前粒子所对应的x(t)的下包络线，记为

3)_8、构建一个维数为(N-1)×N的一阶微分矩阵，记为D，

接着根据D和

确定当前粒子所对应的用于体现

的变化程度的上包络线平稳度函数，记为E_d,i(t)，

并根据D和

确定当前粒子所对应的用于体现

的变化程度的下包络线平稳度函数，记为E_s,i(t)，

其中，符号“|| ||₂”为求矩阵的2-范数符号；

3)_9、计算当前粒子的上包络线适应度函数，记为f_d,i，

并计算当前粒子的下包络线适应度函数，记为f_s,i，

其中，λ为正则化参数，λ＞0，符号“|| ||₁”为求矩阵的1-范数符号；

3)_10、在第α次迭代过程中，将粒子群中下一个待处理的粒子作为当前粒子，然后返回步骤3)_3继续执行，直至粒子群中的所有粒子处理完毕；

3)_11、在第α次迭代过程中，根据L个上包络线适应度函数，从L个粒子中找出使上包络线平稳度函数最小的粒子，将找出的粒子的位置作为粒子群中的所有粒子的上包络线全局最优位置，记为g_d,best；并针对粒子群中的每个粒子，从第1次迭代过程至第α次迭代过程中找出使上包络线平稳度函数最小的迭代序号，将找出的迭代序号对应的迭代过程中该粒子的位置作为该粒子的上包络线局部最优位置，将粒子群中的第i个粒子的上包络线局部最优位置记为p_d,i,best；同样，根据L个下包络线适应度函数，从L个粒子中找出使下包络线平稳度函数最小的粒子，将找出的粒子的位置作为粒子群中的所有粒子的下包络线全局最优位置，记为g_s,best；并针对粒子群中的每个粒子，从第1次迭代过程至第α次迭代过程中找出使下包络线平稳度函数最小的迭代序号，将找出的迭代序号对应的迭代过程中该粒子的位置作为该粒子的下包络线局部最优位置，将粒子群中的第i个粒子的下包络线局部最优位置记为p_s,i,best；

3)_12、判断α＜α_max是否成立，如果α＜α_max成立，则更新上包络线获取过程中粒子群中的每个粒子的速度和位置、下包络线获取过程中粒子群中的每个粒子的速度和位置，将上包络线获取过程中粒子群中的第i个粒子更新后的速度记为v_d,i,back，v_d,i,back＝w×v_d,i,pre+c₁×r₁×(p_d,i,best,pre-p_d,i,pre)+c₂×r₂×(g_d,best,pre-p_d,i,best,pre)，将上包络线获取过程中粒子群中的第i个粒子更新后的位置记为p_d,i,back，p_d,i,back＝p_d,i,pre+v_d,i,pre，将下包络线获取过程中粒子群中的第i个粒子更新后的速度记为v_s,i,back，v_s,i,back＝w×v_s,i,pre+c₁×r₁×(p_s,i,best,pre-p_s,i,pre)+c₂×r₂×(g_s,best,pre-p_s,i,best,pre)，将下包络线获取过程中粒子群中的第i个粒子更新后的位置记为p_s,i,back，p_s,i,back＝p_s,i,pre+v_s,i,pre；然后将上包络线获取过程中粒子群中的每个粒子更新后的速度和位置对应作为下一次迭代过程中上包络线获取过程中该粒子的速度和位置，将下包络线获取过程中粒子群中的每个粒子更新后的速度和位置对应作为下一次迭代过程中下包络线获取过程中该粒子的速度和位置，接着令α＝α+1；再返回步骤3)_2继续执行；其中，w为惯性因子，w的取值范围为[0,1]，v_d,i,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的速度，c₁和c₂均为学习因子，r₁和r₂均为区间[0,1]内的随机数，p_d,i,best,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的上包络线局部最优位置，p_d,i,pre表示上包络线获取过程中粒子群中的第i个粒子更新前的位置，g_d,best,pre表示更新前粒子群中的所有粒子的上包络线全局最优位置，v_s,i,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的速度，p_s,i,best,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的下包络线局部最优位置，p_s,i,pre表示下包络线获取过程中粒子群中的第i个粒子更新前的位置，g_s,best,pre表示更新前粒子群中的所有粒子的下包络线全局最优位置，α＝α+1中的“＝”为赋值符号；

如果α＜α_max不成立，则结束迭代过程，然后将第α次迭代过程中获得的g_d,best赋值给用于改变一个粒子确定的上包络线获取过程所需的DCT基的频带宽度的变化因子，将第α次迭代过程中获得的g_s,best赋值给用于改变一个粒子确定的下包络线获取过程所需的DCT基的频带宽度的变化因子，再按照步骤3)_4至步骤3)_7的过程，以相同的方式执行一次，获得x(t)的最终上包络线即为最佳上包络线和x(t)的最终下包络线即为最佳下包络线。

2.根据权利要求1所述的一种基于稀疏重构最优化算法的信号包络线提取方法，其特征在于所述的步骤2)中的y_d(t)的获取过程为：

步骤2_a1)：对x(t)求一阶导数，将得到的列向量记为d_x1；然后将d_x1中大于零的所有元素的索引值组成的位置向量记为a；

步骤2_a2)：对x(t)求二阶导数，将得到的列向量记为d_x2；然后将d_x2中小于零的所有元素的索引值组成的位置向量记为b；

步骤2_a3)：将a和b中相同的元素提取出来，且只选用从a中提取出的所有元素或只选用从b中提取出的所有元素；然后将选用的所有元素组成列向量，记为γ；再将γ中的每个元素加1，得到新的列向量，记为γ'，γ'中的每个元素即为x(t)中的一个极大值点的索引值；

步骤2_a4)：从x(t)中提取出索引值与γ'中的每个元素相同的采样值；然后将提取出的这些采样值按序排列组成y_d(t)；

所述的步骤2)中的y_s(t)的获取过程为：

步骤2_b1)：取x(t)中的每个采样值的相反数，将得到的新的列向量记为x'(t)；

步骤2_b2)：对x'(t)求一阶导数，将得到的列向量记为d'_x1；然后将d'_x1中大于零的所有元素的索引值组成的位置向量记为a'；

步骤2_b3)：对x'(t)求二阶导数，将得到的列向量记为d'_x2；然后将d'_x2中小于零的所有元素的索引值组成的位置向量记为b'；

步骤2_b4)：将a'和b'中相同的元素提取出来，且只选用从a'中提取出的所有元素或只选用从b'中提取出的所有元素；然后将选用的所有元素组成列向量，记为γ*；接着将γ*中的每个元素加1，得到新的列向量，记为γ*'；再取γ*'中的每个元素的相反数，将得到的新的列向量记为γ'*'，γ'*'中的每个元素即为x(t)中的一个极小值点的索引值；

步骤2_b5)：从x(t)中提取出索引值与γ'*'中的每个元素相同的采样值；然后将提取出的这些采样值按序排列组成y_s(t)。

3.根据权利要求1所述的一种基于稀疏重构最优化算法的信号包络线提取方法，其特征在于所述的步骤3)_4中的Ψ_d,i的获取过程为：将Ψ_d,i中下标为(k,n)的元素记为Ψ_d,i(k,n)，

其中，k和n均为正整数，1≤k≤N，1≤n≤N，cos()为求余弦函数；

所述的步骤3)_4中的Ψ_s,i的获取过程为：将Ψ_s,i中下标为(k,n)的元素记为Ψ_s,i(k,n)，

4.根据权利要求1所述的一种基于稀疏重构最优化算法的信号包络线提取方法，其特征在于所述的步骤3)_6中η_d,i的求解过程为：

3)_6_a1、令k1表示求解迭代次数，令η_d,i的初始值为[0,0,...,0]^T _N×1；其中，k1为正整数，k1的初始值为1；

3)_6_b1、构造第k1次求解迭代过程中的惩罚函数Α(η_d,i,u_d,i)，Α(η_d,i,u_d,i)的表达式为：

并将Α(η_d,i,u_d,i)的惩罚因子记为τ，将Α(η_d,i,u_d,i)的衰减系数记为ρ；其中，q为正整数，1≤q≤N，u_d,i(q)表示所有元素均大于零且维数为N×1的列向量u_d,i中的第q个元素，且满足-u_d,i(q)≤η_d,i(q)≤+u_d,i(q)，u_d,i(q)的初始值为1，η_d,i(q)表示η_d,i中的第q个元素，τ的初始值大于0，ρ的值为区间(0,1)之间的数；

3)_6_c1、在Α(η_d,i,u_d,i)基础上构造第k1次求解迭代过程中的无约束子问题T(η_d,i,u_d,i)，T(η_d,i,u_d,i)的表达式为：

其中，符号“|| ||₂”为求矩阵的2-范数符号，λ为正则化参数，λ＞0；

3)_6_d1、利用共轭梯度法求解T(η_d,i,u_d,i)，得到η_d,i在第k1次求解迭代后的解值；

3)_6_e1、在第k1次求解迭代过程中，判断τ×Α(η_d,i,u_d,i)≤ε是否满足，如果满足，则停止求解迭代过程，并将η_d,i在第k1次求解迭代后的解值作为η_d,i的最终值；如果不满足，则将η_d,i在第k1次求解迭代后的解值作为下一次求解迭代过程中的值，然后令k1＝k1+1，并令τ＝ρ×τ，再返回步骤3)_6_b1继续执行；其中，k1＝k1+1和τ＝ρ×τ中的“＝”为赋值符号，ε为误差终止值，0≤ε≤1；

所述的步骤3)_6中η_s,i的求解过程为：

3)_6_a2、令k2表示求解迭代次数，令η_s,i的初始值为[0,0,...,0]^T _N×1；其中，k2为正整数，k2的初始值为1；

3)_6_b2、构造第k2次求解迭代过程中的惩罚函数Α(η_s,i,u_s,i)，Α(η_s,i,u_s,i)的表达式为：

并将Α(η_s,i,u_s,i)的惩罚因子记为τ'，将Α(η_s,i,u_s,i)的衰减系数记为ρ'；其中，q为正整数，1≤q≤N，u_s,i(q)表示所有元素均大于零且维数为N×1的列向量u_s,i中的第q个元素，且满足-u_s,i(q)≤η_s,i(q)≤+u_s,i(q)，u_s,i(q)的初始值为1，η_s,i(q)表示η_s,i中的第q个元素，τ'的初始值大于0，ρ'的值为区间(0,1)之间的数；

3)_6_c2、在Α(η_s,i,u_s,i)基础上构造第k2次求解迭代过程中的无约束子问题T(η_s,i,u_s,i)，T(η_s,i,u_s,i)的表达式为：

其中，符号“|| ||₂”为求矩阵的2-范数符号，λ为正则化参数，λ＞0；

3)_6_d2、利用共轭梯度法求解T(η_s,i,u_s,i)，得到η_s,i在第k2次求解迭代后的解值；

3)_6_e2、在第k2次求解迭代过程中，判断τ'×Α(η_s,i,u_s,i)≤ε是否满足，如果满足，则停止求解迭代过程，并将η_s,i在第k2次求解迭代后的解值作为η_s,i的最终值；如果不满足，则将η_s,i在第k2次求解迭代后的解值作为下一次求解迭代过程中的值，然后令k2＝k2+1，并令τ'＝ρ'×τ'，再返回步骤3)_6_b2继续执行；其中，k2＝k2+1和τ'＝ρ'×τ'中的“＝”为赋值符号，ε为误差终止值，0≤ε≤1。

Images