CN104504181B - 一种基于稀疏复原的信号包络线提取方法 - Google Patents

Abstract

本发明公开了一种基于稀疏复原的信号包络线提取方法，其先找出待提取包络线的信号中的所有极大值点和所有极小值点，对应构成极大值点向量和极小值点向量；然后构建一个DCT基，从DCT基中提取出行号与每个极大值点的下标一致的每行元素构成一个矩阵，并从DCT基中提取出行号与每个极小值点的下标一致的每行元素构成一个矩阵；接着将极大值点向量作为观测向量、对应的矩阵作为感知矩阵获取上包络线，并将极小值点向量作为观测向量、对应的矩阵作为感知矩阵获取下包络线；最后根据上包络线和下包络线各自的平滑度，并结合DCT基的带宽的变化因子获取最佳上包络线和最佳下包络线；优点是不仅能够有效地提高包络线的精度，而且能够有效地抑制端点效应。

Description

一种基于稀疏复原的信号包络线提取方法

技术领域

本发明涉及一种信号处理技术，尤其是涉及一种基于稀疏复原的信号包络线提取方法。

背景技术

包络分析方法的核心是把调制在中、高频带的低频故障信息，解调到低频进行分析处理，从而提取出故障信息。在机械故障诊断过程中，由于受旋转机械的干扰信号和噪声的影响，提高了信噪比，因而在故障诊断，尤其是齿轮箱、滚动轴承故障诊断中，包络分析方法具有其他故障检测方法不可替代的作用，是目前诊断轴承和齿轮故障的最有效方法。

在信号处理中，目前常用的包络分析方法有：Hilbert(希尔伯特)变换、三次样条插值等。基于Hilbert变换提取信号包络线的方法可以有效地提取调制频率及具有一定的抗噪性，但是随着信号信噪比的变小，Hilbert变换所得的包络误差会逐渐变大，导致所提取的信号包络线不光滑，从而影响信号包络线提取的精度。而基于三次样条插值提取信号包络线的方法提取出的信号包络线不仅有很好的光滑度，而且当节点逐渐加密时，能够很好的提高信号包络线提取的精度，但是样条插值函数需要数据向量两端数据的一阶和二阶导数，而由数据曲线得不到所需要的端点处信息，因此包络线在端点会发生大的摆动，形成非常棘手的端点效应问题。

发明内容

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

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

①假定待提取包络线的信号为x，则将x以行向量的形式表示为x＝[x₁ x₂ … x_N-1x_N]，其中，在此符号“[ ]”为向量表示符号，N表示x的采样点数，x₁表示x中的第1个采样值，x₂表示x中的第2个采样值，x_N-1表示x中的第N-1个采样值，x_N表示x中的第N个采样值；

②找出x中的所有极大值点和所有极小值点，然后将从x中找出的所有极大值点按序排列构成一个x的极大值点向量，记为P_a，并将从x中找出的所有极小值点按序排列构成一个x的极小值点向量，记为P_b；

③根据所要构建的DCT基的阶数和用于改变所要构建的DCT基的带宽的变化因子，构建一个DCT基，记为Ψ，其中，Ψ为一个N阶方阵；

④从Ψ中提取出行号与x中的每个极大值点的下标一致的每行元素，然后将提取出的所有行按行号顺序排列构成一个维数为K₁×N的矩阵，记为H₁，其中，K₁表示x中的极大值点的总个数，1≤K₁<N；

同样，从Ψ中提取出行号与x中的每个极小值点的下标一致的每行元素，然后将提取出的所有行按行号顺序排列构成一个维数为K₂×N的矩阵，记为H₂，其中，K₂表示x中的极小值点的总个数，1≤K₂<N；

⑤将P_a作为观测向量，将H₁作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的上包络线，记为x_a；

同样，将P_b作为观测向量，将H₂作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的下包络线，记为x_b；

⑥构建一个维数为(N-1)×N的差分矩阵，记为D，D中第i行第i列的元素的值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素外的所有元素的值均为0；然后根据x_a和D获取x_a的平滑度，记为H_a，同样根据x_b和D获取x_b的平滑度，记为H_b；

⑦判断H_a是否小于min_H_a，如果是，则令min_H_a＝H_a，并将x_a作为x的最佳上包络线，然后执行步骤⑧，否则，直接执行步骤⑧，其中，min_H_a的初始值为无穷大；

⑧判断H_b是否小于min_H_b，如果是，则令min_H_b＝H_b，并将x_b作为x的最佳下包络线，然后执行步骤⑨，否则，直接执行步骤⑨，其中，min_H_b的初始值为无穷大；

⑨判断m是否等于md，如果是，则分别输出x的最佳上包络线和x的最佳下包络线，否则，令m＝m+1，然后返回步骤③继续执行，其中，m表示用于改变所要构建的DCT基的带宽的变化因子，m的初始值为1，1≤m≤md，md表示设定的变化因子最大值，m＝m+1中的“＝”为赋值符号。

所述的步骤②中x的极大值点向量的获取过程为：

②-1a、对x求一阶差分，得到x的一阶差分向量，记为dx；

②-2a、根据dx中的每个元素的值，获取一个元素的值为1或0的新向量，记为A_x；对于dx中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小于或等于零，则将该元素的值置为0，其中，1≤j≤N；

②-3a、对A_x求一阶差分，得到A_x的一阶差分向量，记为B_x；

②-4a、在B_x中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构成一个位置向量，记为W1_x；再将W1_x中的每个元素的值加1，得到新的位置向量，记为W1_x'，W1_x'中的任一个元素的值为x中的一个极大值点的下标；

②-5a、根据W1_x'中的每个元素的值，在x中找出所有极大值点，对于W1_x'中的任一个元素，其值为x中的一个极大值点的下标；然后将从x中找出的所有极大值点按序排列构成一个x的极大值点向量，记为P_a；

所述的步骤②中x的极小值点向量的获取过程为：

②-1b、对x进行取反操作，得到x的反向量，记为然后对求一阶差分，得到的一阶差分向量，记为

②-2b、根据中的每个元素的值，获取一个元素的值为1或0的新向量，记为对于中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小于或等于零，则将该元素的值置为0，其中，1≤j≤N；

②-3b、对求一阶差分，得到的一阶差分向量，记为

②-4b、在中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构成一个位置向量，记为接着将中的每个元素的值加1，得到新的位置向量，记为再对进行取反操作，得到的反向量，记为中的任一个元素的值为x中的一个极小值点的下标；

②-5b、根据中的每个元素的值，在x中找出所有极小值点，对于中的任一个元素，其值为x中的一个极小值点的下标；然后将从x中找出的所有极小值点按序排列构成一个x的极小值点向量，记为P_b。

所述的步骤③的具体过程为：

③-1、假定所要构建的DCT基的阶数为N阶，并令m表示用于改变所要构建的DCT基的带宽的变化因子，其中，m的初始值为1，1≤m≤md，md表示设定的变化因子最大值；

③-2、构建一个N阶方阵，记为Ψ，将Ψ中第p行第q列的元素的值记为Ψ(p,q)，其中，p和q的初始值均为1，1≤p≤N,1≤q≤N，cos()为求余弦函数；

③-3、将步骤③-2构建的Ψ作为DCT基。

所述的步骤⑤中x的上包络线的获取过程为：

⑤-1a、令ka表示迭代次数，令P1表示维数为K₁的向量，令B1表示维数为K₁×K₁的矩阵，其中，ka的初始值为1，1≤ka≤K₁；

⑤-2a、计算H₁中的每一列与第ka-1次迭代的残差值r_ka-1之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H₁中的一列的列号作为P1中的第ka个元素的值，并将N个相关系数中的最大值对应的H₁中的一列作为B1中的第ka列；再在H₁中将N个相关系数中的最大值对应的H₁中的一列剔除；

其中，H₁中的第j列与第ka-1次迭代的残差值r_ka-1之间的相关系数为H₁中的第j列与第ka-1次迭代的残差值r_ka-1之间的内积，1≤j≤N，当ka＝1时取r_ka-1＝P_a；

⑤-3a、利用最小二乘法，计算r_ka-1在B1中的每一列上的稀疏映射系数，然后将所有稀疏映射系数按序构成稀疏映射系数序列，记为a_ka；

⑤-4a、根据B1和a_ka，计算第ka次迭代的残差值，记为r_ka，r_ka＝r_ka-1-B1×a_ka；

⑤-5a、判断ka是否等于K₁，如果是，则执行步骤⑤-6a，否则，令ka＝ka+1，然后返回步骤⑤-2a继续执行，其中，ka＝ka+1中的“＝”为赋值符号；

⑤-6a、根据P1和a_ka获取稀疏信号，记为共包含K₁个非零信号值和N-K₁个零值，中下标与P1中的第i个元素的值一致的元素的值为非零信号值且等于a_ka中的第i个元素的值，中下标与P1中的任一个元素的值不一致的元素的值等于0，其中，1≤i≤K1；

⑤-7a、根据和Ψ，计算由P_a恢复出的信号，记为x_a，然后将x_a作为x的信号上包络线；

同样，将P_b作为观测向量，将H₂作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的下包络线，记为x_b；

所述的步骤⑤中x的下包络线的获取过程为：

⑤-1b、令kb表示迭代次数，令P2表示维数为K₂的向量，令B2表示维数为K₂×K₂的矩阵，其中，kb的初始值为1，1≤kb≤K₂；

⑤-2b、计算H₂中的每一列与第kb-1次迭代的残差值r_kb-1之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H₂中的一列的列号作为P2中的第kb个元素的值，并将N个相关系数中的最大值对应的H₂中的一列作为B2中的第kb列；再在H₂中将N个相关系数中的最大值对应的H₂中的一列剔除；

其中，H₂中的第j列与第kb-1次迭代的残差值r_kb-1之间的相关系数为H₂中的第j列与第kb-1次迭代的残差值r_kb-1之间的内积，1≤j≤N，当kb＝1时取r_kb-1＝P_b；

⑤-3b、利用最小二乘法，计算r_kb-1在B2中的每一列上的稀疏映射系数，然后将所有稀疏映射系数按序构成稀疏映射系数序列，记为a_kb；

⑤-4b、根据B2和a_kb，计算第kb次迭代的残差值，记为r_kb，r_kb＝r_kb-1-B2×a_kb；

⑤-5b、判断kb是否等于K₂，如果是，则执行步骤⑤-6b，否则，令kb＝kb+1，然后返回步骤⑤-2b继续执行，其中，kb＝kb+1中的“＝”为赋值符号；

⑤-6b、根据P2和a_kb获取稀疏信号，记为共包含K₂个非零信号值和N-K₂个零值，中下标与P2中的第i个元素的值一致的元素的值为非零信号值且等于a_kb中的第i个元素的值，中下标与P2中的任一个元素的值不一致的元素的值等于0，其中，1≤i≤K2；

⑤-7b、根据和Ψ，计算由P_b恢复出的信号，记为x_b，然后将x_b作为x的信号下包络线。

所述的步骤⑥的具体过程为：

⑥-1、构建一个维数为(N-1)×N的差分矩阵，记为D，D中第i行第i列的元素的值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素外的所有元素的值均为0，即

⑥-2、根据x_a和D计算的差分向量，记为dx_a，dx_a＝x_a×D；同样，根据x_b和D计算x_b的差分向量，记为dx_b，dx_b＝x_b×D；

⑥-3、计算dx_a中的所有元素的值的平方和，记为H_a，然后将H_a作为x_a的平滑度；同样，计算dx_b中的所有元素的值的平方和，记为H_b，然后将H_b作为x_b的平滑度。

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

1)本发明方法引入一个用于改变所要构建的DCT基的带宽的变换因子，可通过改变变换因子的值来构建不同的稀疏基即DCT基，这种方式可以利用正交匹配追踪算法自适应的提取出信号的最佳上包络线和最佳下包络线，因此本发明方法具有很好的自适应性，能够有效地提高包络线提取的精度。

2)本发明方法采用的稀疏复原方法是正交匹配追踪算法，分别以极大值点向量和极小值点向量作为正交匹配追踪算法的观测矩阵，相应求出信号的最佳上包络线和最佳下包络线，正交匹配追踪算法不仅在程序上易于实现，而且计算成本低、速度快。

3)本发明方法采用的是基于稀疏复原方法拟合出极大值点和极小值点的曲线，即将极大值点向量作为观测向量进行稀疏复原提取出信号的上包络线，并将极小值点向量作为观测向量进行稀疏复原提取出信号的下包络线，这就避免了三次样条插值拟合曲线所存在的端点效应问题，因此本发明方法可以有效的抑制端点效应，该特性在经验模态分解中可以得到很好的应用。

附图说明

图1为本发明方法的流程框图；

图2为采用本发明方法提取出的信号的最佳上包络线的示意图；

图3为采用本发明方法提取出的信号的最佳下包络线的示意图。

具体实施方式

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

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

①假定待提取包络线的信号为x，则将x以行向量的形式表示为x＝[x₁ x₂ … x_N-1x_N]，其中，在此符号“[ ]”为向量表示符号，N表示x的采样点数，N的具体值一般根据处理的信号而定，通常采样点数越多越好，但是采样点数越多则计算复杂度会增加，因此在实际处理过程中一般都折衷考虑，如本实施例在仿真时取N＝256，x₁表示x中的第1个采样值，x₂表示x中的第2个采样值，x_N-1表示x中的第N-1个采样值，x_N表示x中的第N个采样值。

②找出x中的所有极大值点和所有极小值点，然后将从x中找出的所有极大值点按序排列构成一个x的极大值点向量，记为P_a，并将从x中找出的所有极小值点按序排列构成一个x的极小值点向量，记为P_b。

在本实施例中，步骤②中x的极大值点向量的获取过程为：

②-1a、对x求一阶差分，得到x的一阶差分向量，记为dx。在此，对x求一阶差分采用现有技术。

②-2a、根据dx中的每个元素的值，获取一个元素的值为1或0的新向量，记为A_x；对于dx中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小于或等于零，则将该元素的值置为0，其中，1≤j≤N。

②-3a、对A_x求一阶差分，得到A_x的一阶差分向量，记为B_x。

②-4a、在B_x中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构成一个位置向量，记为W1_x；再将W1_x中的每个元素的值加1，得到新的位置向量，记为W1_x'，W1_x'中的任一个元素的值为x中的一个极大值点的下标。

②-5a、根据W1_x'中的每个元素的值，在x中找出所有极大值点，对于W1_x'中的任一个元素，其值为x中的一个极大值点的下标；然后将从x中找出的所有极大值点按序排列构成一个x的极大值点向量，记为P_a。

在本实施例中，步骤②中x的极小值点向量的获取过程与x的极大值点向量的获取过程大致相同，只是先将x取反，然后按照步骤②-1a至步骤②-4a的过程，以相同的方式获取新的位置向量，接着对新的位置向量取反，再按照步骤②-5a的过程，以相同的方式获取x的极大值点向量，即步骤②中x的极小值点向量的获取过程为：

②-1b、对x进行取反操作，得到x的反向量，记为然后对求一阶差分，得到的一阶差分向量，记为在此，对x进行取反操作，即对x中的每个元素取其相反数，如2的相反数为-2。

②-2b、根据中的每个元素的值，获取一个元素的值为1或0的新向量，记为对于中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小于或等于零，则将该元素的值置为0，其中，1≤j≤N。

②-3b、对求一阶差分，得到的一阶差分向量，记为

②-4b、在中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构成一个位置向量，记为接着将中的每个元素的值加1，得到新的位置向量，记为再对进行取反操作，得到的反向量，记为中的任一个元素的值为x中的一个极小值点的下标。在此，对进行取反操作，即对中的每个元素取其相反数。

②-5b、根据中的每个元素的值，在x中找出所有极小值点，对于中的任一个元素，其值为x中的一个极小值点的下标；然后将从x中找出的所有极小值点按序排列构成一个x的极小值点向量，记为P_b。

③根据所要构建的DCT基的阶数和用于改变所要构建的DCT基的带宽的变化因子，构建一个DCT基，记为Ψ，其中，Ψ为一个N阶方阵。

在此具体实施例中，步骤③的具体过程为：

③-1、假定所要构建的DCT基的阶数为N阶，并令m表示用于改变所要构建的DCT基的带宽的变化因子，其中，m的初始值为1，1≤m≤md，md表示设定的变化因子最大值，在本实施例中取md＝20，该值是在本发明方法的基础上通过大量实验获得的实验结果。

③-2、构建一个N阶方阵，记为Ψ，将Ψ中第p行第q列的元素的值记为Ψ(p,q)，其中，p和q的初始值均为1，1≤p≤N,1≤q≤N，cos()为求余弦函数。

③-3、将步骤③-2构建的Ψ作为DCT基。

④从Ψ中提取出行号与x中的每个极大值点的下标一致的每行元素，然后将提取出的所有行按行号顺序排列构成一个维数为K₁×N的矩阵，记为H₁，其中，K₁表示x中的极大值点的总个数，亦表示提取出的行的总行数，1≤K₁<N。

同样，从Ψ中提取出行号与x中的每个极小值点的下标一致的每行元素，然后将提取出的所有行按行号顺序排列构成一个维数为K₂×N的矩阵，记为H₂，其中，K₂表示x中的极小值点的总个数，亦表示提取出的行的总行数，1≤K₂<N。

⑤将P_a作为观测向量，将H₁作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的上包络线，记为x_a。

在此具体实施例中，步骤⑤中x的上包络线的获取过程为：

⑤-1a、令ka表示迭代次数，令P1表示维数为K₁的向量，令B1表示维数为K₁×K₁的矩阵，其中，ka的初始值为1，1≤ka≤K₁。

⑤-2a、计算H₁中的每一列与第ka-1次迭代的残差值r_ka-1之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H₁中的一列的列号作为P1中的第ka个元素的值，并将N个相关系数中的最大值对应的H₁中的一列作为B1中的第ka列；再在H₁中将N个相关系数中的最大值对应的H₁中的一列剔除。

其中，H₁中的第j列与第ka-1次迭代的残差值r_ka-1之间的相关系数为H₁中的第j列与第ka-1次迭代的残差值r_ka-1之间的内积，1≤j≤N，当ka＝1时取r_ka-1＝P_a。

⑤-3a、利用最小二乘法，计算r_ka-1在B1中的每一列上的稀疏映射系数，然后将所有稀疏映射系数按序构成稀疏映射系数序列，记为a_ka。

⑤-4a、根据B1和a_ka，计算第ka次迭代的残差值，记为r_ka，r_ka＝r_ka-1-B1×a_ka。

⑤-5a、判断ka是否等于K₁，如果是，则执行步骤⑤-6a，否则，令ka＝ka+1，然后返回步骤⑤-2a继续执行，其中，ka＝ka+1中的“＝”为赋值符号。

⑤-6a、根据P1和a_ka获取稀疏信号，记为共包含K₁个非零信号值和N-K₁个零值，中下标与P1中的第i个元素的值一致的元素的值为非零信号值且等于a_ka中的第i个元素的值，中下标与P1中的任一个元素的值不一致的元素的值等于0，其中，1≤i≤K1。

当ka等于K₁时，P1和a_ka中均包含有K₁元素，假设P1＝[P1₁,P1₂,...,P1_i,...,P1_K1]，a_ka＝[a_ka1,a_ka2,...,a_kai,...,a_kaK1]，则中下标为P1₁的元素的值为非零信号值且等于a_ka1，中下标为P1₂的元素的值为非零信号值且等于a_ka2，中下标为P1_i的元素的值为非零信号值且等于a_kai，中下标为P1_K1的元素的值为非零信号值且等于a_kaK1，中的其他元素的值为0。

⑤-7a、根据和Ψ，计算由P_a恢复出的信号，记为x_a，然后将x_a作为x的信号上包络线。

同样，将P_b作为观测向量，将H₂作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的下包络线，记为x_b。

在此具体实施例中，步骤⑤中x的下包络线的获取过程为：

⑤-1b、令kb表示迭代次数，令P2表示维数为K₂的向量，令B2表示维数为K₂×K₂的矩阵，其中，kb的初始值为1，1≤kb≤K₂。

⑤-2b、计算H₂中的每一列与第kb-1次迭代的残差值r_kb-1之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H₂中的一列的列号作为P2中的第kb个元素的值，并将N个相关系数中的最大值对应的H₂中的一列作为B2中的第kb列；再在H₂中将N个相关系数中的最大值对应的H₂中的一列剔除。

其中，H₂中的第j列与第kb-1次迭代的残差值r_kb-1之间的相关系数为H₂中的第j列与第kb-1次迭代的残差值r_kb-1之间的内积，1≤j≤N，当kb＝1时取r_kb-1＝P_b。

⑤-3b、利用最小二乘法，计算r_kb-1在B2中的每一列上的稀疏映射系数，然后将所有稀疏映射系数按序构成稀疏映射系数序列，记为a_kb。

⑤-4b、根据B2和a_kb，计算第kb次迭代的残差值，记为r_kb，r_kb＝r_kb-1-B2×a_kb。

⑤-5b、判断kb是否等于K₂，如果是，则执行步骤⑤-6b，否则，令kb＝kb+1，然后返回步骤⑤-2b继续执行，其中，kb＝kb+1中的“＝”为赋值符号。

⑤-6b、根据P2和a_kb获取稀疏信号，记为共包含K₂个非零信号值和N-K₂个零值，中下标与P2中的第i个元素的值一致的元素的值为非零信号值且等于a_kb中的第i个元素的值，中下标与P2中的任一个元素的值不一致的元素的值等于0，其中，1≤i≤K2。

当kb等于K₂时，P2和a_kb中均包含有K₂元素，假设P2＝[P2₁,P2₂,...,P2_i,...,P2_K2]，a_kb＝[a_kb1,a_kb2,...,a_kbi,...,a_kbK2]，则中下标为P2₁的元素的值为非零信号值且等于a_kb1，中下标为P2₂的元素的值为非零信号值且等于a_kb2，中下标为P2_i的元素的值为非零信号值且等于a_kbi，中下标为P2_K2的元素的值为非零信号值且等于a_kbK2，中的其他元素的值为0。

⑤-7b、根据和Ψ，计算由P_b恢复出的信号，记为x_b，然后将x_b作为x的信号下包络线。

⑥构建一个维数为(N-1)×N的差分矩阵，记为D，D中第i行第i列的元素的值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素外的所有元素的值均为0；然后根据x_a和D获取x_a的平滑度，记为H_a，同样根据x_b和D获取x_b的平滑度，记为H_b。

在此具体实施例中，步骤⑥的具体过程为：

⑥-1、构建一个维数为(N-1)×N的差分矩阵，记为D，D中第i行第i列的元素的值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素外的所有元素的值均为0，即

⑥-2、根据x_a和D计算x_a的差分向量，记为dx_a，dx_a＝x_a×D；同样，根据x_b和D计算x_b的差分向量，记为dx_b，dx_b＝x_b×D。

⑥-3、计算dx_a中的所有元素的值的平方和，记为H_a，然后将H_a作为x_a的平滑度；同样，计算dx_b中的所有元素的值的平方和，记为H_b，然后将H_b作为x_b的平滑度。

⑦判断H_a是否小于min_H_a，如果是，则令min_H_a＝H_a，并将x_a作为x的最佳上包络线，然后执行步骤⑧，否则，直接执行步骤⑧，其中，min_H_a的初始值为无穷大。

⑧判断H_b是否小于min_H_b，如果是，则令min_H_b＝H_b，并将x_b作为x的最佳下包络线，然后执行步骤⑨，否则，直接执行步骤⑨，其中，min_H_b的初始值为无穷大。

⑨判断m是否等于md，如果是，则分别输出x的最佳上包络线和x的最佳下包络线，否则，令m＝m+1，然后返回步骤③继续执行，其中，m表示用于改变所要构建的DCT基的带宽的变化因子，m的初始值为1，1≤m≤md，md表示设定的变化因子最大值，在本实施例中取md＝20，该值是在本发明方法的基础上通过大量实验获得的实验结果，m＝m+1中的“＝”为赋值符号。

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

假设一维信号x(t)＝0.8exp(-1.5t)cos(2πf₁t+0.1)+0.9exp(-2t)cos(2πf₂t+0.2)，并假设采样点数N为256，信号采样频率为f_s＝800Hz，式中f₁＝50Hz，f₂＝100Hz。利用本发明方法分别提取一维信号x(t)的最佳上包络线和最佳下包络线，图2给出了采用本发明方法提取出的信号x(t)的最佳上包络线的示意图，图3给出了采用本发明方法提取出的信号x(t)的最佳下包络线的示意图。图2中的m为用于表示改变DCT基的带宽的变化因子，dd表示的是采用本发明方法拟合出的曲线的平滑程度。由于包络线具有变化缓慢的特性，因此dd越小越好。图2中的稀疏复原即为本发明方法，从图2中可以看出利用本发明方法可以准确地提取出信号的最佳上包络线，并且端点处的包络线也能够准确提取出来，由此证明了本发明方法的可行性和有效性。图3中的稀疏复原即为本发明方法，从图3中可以看出利用本发明方法也能准确地提取出信号的最佳下包络线，再次验证了本发明方法的可行性和有效性。

Claims (4)

1.一种基于稀疏复原的信号包络线提取方法，其特征在于包括以下步骤：

①假定待提取包络线的信号为x，则将x以行向量的形式表示为x＝[x₁ x₂ … x_N-1 x_N]，其中，在此符号“[]”为向量表示符号，N表示x的采样点数，x₁表示x中的第1个采样值，x₂表示x中的第2个采样值，x_N-1表示x中的第N-1个采样值，x_N表示x中的第N个采样值；

②找出x中的所有极大值点和所有极小值点，然后将从x中找出的所有极大值点按序排列构成一个x的极大值点向量，记为P_a，并将从x中找出的所有极小值点按序排列构成一个x的极小值点向量，记为P_b；

③根据所要构建的DCT基的阶数和用于改变所要构建的DCT基的带宽的变化因子，构建一个DCT基，记为Ψ，其中，Ψ为一个N阶方阵；

所述的步骤③的具体过程为：

③-1、假定所要构建的DCT基的阶数为N阶，并令m表示用于改变所要构建的DCT基的带宽的变化因子，其中，m的初始值为1，1≤m≤md，md表示设定的变化因子最大值；

③-2、构建一个N阶方阵，记为Ψ，将Ψ中第p行第q列的元素的值记为Ψ(p,q)，其中，p和q的初始值均为1，1≤p≤N,1≤q≤N，cos()为求余弦函数；

③-3、将步骤③-2构建的Ψ作为DCT基；

④从Ψ中提取出行号与x中的每个极大值点的下标一致的每行元素，然后将提取出的所有行按行号顺序排列构成一个维数为K₁×N的矩阵，记为H₁，其中，K₁表示x中的极大值点的总个数，1≤K₁<N；

同样，从Ψ中提取出行号与x中的每个极小值点的下标一致的每行元素，然后将提取出的所有行按行号顺序排列构成一个维数为K₂×N的矩阵，记为H₂，其中，K₂表示x中的极小值点的总个数，1≤K₂<N；

⑤将P_a作为观测向量，将H₁作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的上包络线，记为x_a；

同样，将P_b作为观测向量，将H₂作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的下包络线，记为x_b；

⑥构建一个维数为(N-1)×N的差分矩阵，记为D，D中第i行第i列的元素的值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素外的所有元素的值均为0；然后根据x_a和D获取x_a的平滑度，记为H_a，同样根据x_b和D获取x_b的平滑度，记为H_b；

⑦判断H_a是否小于min_H_a，如果是，则令min_H_a＝H_a，并将x_a作为x的最佳上包络线，然后执行步骤⑧，否则，直接执行步骤⑧，其中，min_H_a的初始值为无穷大；

⑧判断H_b是否小于min_H_b，如果是，则令min_H_b＝H_b，并将x_b作为x的最佳下包络线，然后执行步骤⑨，否则，直接执行步骤⑨，其中，min_H_b的初始值为无穷大；

⑨判断m是否等于md，如果是，则分别输出x的最佳上包络线和x的最佳下包络线，否则，令m＝m+1，然后返回步骤③继续执行，其中，m表示用于改变所要构建的DCT基的带宽的变化因子，m的初始值为1，1≤m≤md，md表示设定的变化因子最大值，m＝m+1中的“＝”为赋值符号。

2.根据权利要求1所述的一种基于稀疏复原的信号包络线提取方法，其特征在于所述的步骤②中x的极大值点向量的获取过程为：

②-1a、对x求一阶差分，得到x的一阶差分向量，记为dx；

②-2a、根据dx中的每个元素的值，获取一个元素的值为1或0的新向量，记为A_x；对于dx中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小于或等于零，则将该元素的值置为0，其中，1≤j≤N；

②-3a、对A_x求一阶差分，得到A_x的一阶差分向量，记为B_x；

②-4a、在B_x中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构成一个位置向量，记为W1_x；再将W1_x中的每个元素的值加1，得到新的位置向量，记为W1_x'，W1_x'中的任一个元素的值为x中的一个极大值点的下标；

②-5a、根据W1_x'中的每个元素的值，在x中找出所有极大值点，对于W1_x'中的任一个元素，其值为x中的一个极大值点的下标；然后将从x中找出的所有极大值点按序排列构成一个x的极大值点向量，记为P_a；

所述的步骤②中x的极小值点向量的获取过程为：

②-1b、对x进行取反操作，得到x的反向量，记为然后对求一阶差分，得到的一阶差分向量，记为

②-2b、根据中的每个元素的值，获取一个元素的值为1或0的新向量，记为对于中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小于或等于零，则将该元素的值置为0，其中，1≤j≤N；

②-3b、对求一阶差分，得到的一阶差分向量，记为

②-4b、在中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构成一个位置向量，记为接着将中的每个元素的值加1，得到新的位置向量，记为再对进行取反操作，得到的反向量，记为中的任一个元素的值为x中的一个极小值点的下标；

②-5b、根据中的每个元素的值，在x中找出所有极小值点，对于中的任一个元素，其值为x中的一个极小值点的下标；然后将从x中找出的所有极小值点按序排列构成一个x的极小值点向量，记为P_b。

3.根据权利要求1所述的一种基于稀疏复原的信号包络线提取方法，其特征在于所述的步骤⑤中x的上包络线的获取过程为：

⑤-1a、令ka表示迭代次数，令P1表示维数为K₁的向量，令B1表示维数为K₁×K₁的矩阵，其中，ka的初始值为1，1≤ka≤K₁；

⑤-2a、计算H₁中的每一列与第ka-1次迭代的残差值r_ka-1之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H₁中的一列的列号作为P1中的第ka个元素的值，并将N个相关系数中的最大值对应的H₁中的一列作为B1中的第ka列；再在H₁中将N个相关系数中的最大值对应的H₁中的一列剔除；

其中，H₁中的第j列与第ka-1次迭代的残差值r_ka-1之间的相关系数为H₁中的第j列与第ka-1次迭代的残差值r_ka-1之间的内积，1≤j≤N，当ka＝1时取r_ka-1＝P_a；

⑤-3a、利用最小二乘法，计算r_ka-1在B1中的每一列上的稀疏映射系数，然后将所有稀疏映射系数按序构成稀疏映射系数序列，记为a_ka；

⑤-4a、根据B1和a_ka，计算第ka次迭代的残差值，记为r_ka，r_ka＝r_ka-1-B1×a_ka；

⑤-5a、判断ka是否等于K₁，如果是，则执行步骤⑤-6a，否则，令ka＝ka+1，然后返回步骤⑤-2a继续执行，其中，ka＝ka+1中的“＝”为赋值符号；

⑤-6a、根据P1和a_ka获取稀疏信号，记为共包含K₁个非零信号值和N-K₁个零值，中下标与P1中的第i个元素的值一致的元素的值为非零信号值且等于a_ka中的第i个元素的值，中下标与P1中的任一个元素的值不一致的元素的值等于0，其中，1≤i≤K1；

⑤-7a、根据和Ψ，计算由P_a恢复出的信号，记为x_a，然后将x_a作为x的信号上包络线；

同样，将P_b作为观测向量，将H₂作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为x的下包络线，记为x_b；

所述的步骤⑤中x的下包络线的获取过程为：

⑤-1b、令kb表示迭代次数，令P2表示维数为K₂的向量，令B2表示维数为K₂×K₂的矩阵，其中，kb的初始值为1，1≤kb≤K₂；

⑤-2b、计算H₂中的每一列与第kb-1次迭代的残差值r_kb-1之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H₂中的一列的列号作为P2中的第kb个元素的值，并将N个相关系数中的最大值对应的H₂中的一列作为B2中的第kb列；再在H₂中将N个相关系数中的最大值对应的H₂中的一列剔除；

其中，H₂中的第j列与第kb-1次迭代的残差值r_kb-1之间的相关系数为H₂中的第j列与第kb-1次迭代的残差值r_kb-1之间的内积，1≤j≤N，当kb＝1时取r_kb-1＝P_b；

⑤-3b、利用最小二乘法，计算r_kb-1在B2中的每一列上的稀疏映射系数，然后将所有稀疏映射系数按序构成稀疏映射系数序列，记为a_kb；

⑤-4b、根据B2和a_kb，计算第kb次迭代的残差值，记为r_kb，r_kb＝r_kb-1-B2×a_kb；

⑤-5b、判断kb是否等于K₂，如果是，则执行步骤⑤-6b，否则，令kb＝kb+1，然后返回步骤⑤-2b继续执行，其中，kb＝kb+1中的“＝”为赋值符号；

⑤-6b、根据P2和a_kb获取稀疏信号，记为共包含K₂个非零信号值和N-K₂个零值，中下标与P2中的第i个元素的值一致的元素的值为非零信号值且等于a_kb中的第i个元素的值，中下标与P2中的任一个元素的值不一致的元素的值等于0，其中，1≤i≤K2；

⑤-7b、根据和Ψ，计算由P_b恢复出的信号，记为x_b，然后将x_b作为x的信号下包络线。

4.根据权利要求3所述的一种基于稀疏复原的信号包络线提取方法，其特征在于所述的步骤⑥的具体过程为：

⑥-1、构建一个维数为(N-1)×N的差分矩阵，记为D，D中第i行第i列的元素的值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素外的所有元素的值均为0，即

⑥-2、根据x_a和D计算的差分向量，记为dx_a，dx_a＝x_a×D；同样，根据x_b和D计算x_b的差分向量，记为dx_b，dx_b＝x_b×D；

⑥-3、计算dx_a中的所有元素的值的平方和，记为H_a，然后将H_a作为x_a的平滑度；同样，计算dx_b中的所有元素的值的平方和，记为H_b，然后将H_b作为x_b的平滑度。