一种基于稀疏复原的信号包络线提取方法
【专利摘要】本发明公开了一种基于稀疏复原的信号包络线提取方法，其先找出待提取包络线的信号中的所有极大值点和所有极小值点，对应构成极大值点向量和极小值点向量；然后构建一个DCT基，从DCT基中提取出行号与每个极大值点的下标一致的每行元素构成一个矩阵，并从DCT基中提取出行号与每个极小值点的下标一致的每行元素构成一个矩阵；接着将极大值点向量作为观测向量、对应的矩阵作为感知矩阵获取上包络线，并将极小值点向量作为观测向量、对应的矩阵作为感知矩阵获取下包络线；最后根据上包络线和下包络线各自的平滑度，并结合DCT基的带宽的变化因子获取最佳上包络线和最佳下包络线；优点是不仅能够有效地提高包络线的精度，而且能够有效地抑制端点效应。
【专利说明】-种基于稀疏复原的信号包络线提取方法
【技术领域】
[0001] 本发明设及一种信号处理技术，尤其是设及一种基于稀疏复原的信号包络线提取 方法。
【背景技术】
[0002] 包络分析方法的核屯、是把调制在中、高频带的低频故障信息，解调到低频进行分 析处理，从而提取出故障信息。在机械故障诊断过程中，由于受旋转机械的干扰信号和噪声 的影响，提高了信噪比，因而在故障诊断，尤其是齿轮箱、滚动轴承故障诊断中，包络分析方 法具有其他故障检测方法不可替代的作用，是目前诊断轴承和齿轮故障的最有效方法。
[0003] 在信号处理中，目前常用的包络分析方法有出i化ert (希尔伯特）变换、=次样条 插值等。基于化化ed变换提取信号包络线的方法可W有效地提取调制频率及具有一定的 抗噪性，但是随着信号信噪比的变小，HUbed变换所得的包络误差会逐渐变大，导致所提 取的信号包络线不光滑，从而影响信号包络线提取的精度。而基于=次样条插值提取信号 包络线的方法提取出的信号包络线不仅有很好的光滑度，而且当节点逐渐加密时，能够很 好的提高信号包络线提取的精度，但是样条插值函数需要数据向量两端数据的一阶和二阶 导数，而由数据曲线得不到所需要的端点处信息，因此包络线在端点会发生大的摆动，形成 非常棘手的端点效应问题。
【发明内容】
[0004] 本发明所要解决的技术问题是提供一种基于稀疏复原的信号包络线提取方法，其 不仅能够有效地提高包络线的精度，而且能够有效地抑制端点效应。
[0005] 本发明解决上述技术问题所采用的技术方案为；一种基于稀疏复原的信号包络线 提取方法，其特征在于包括W下步骤：
[0006] ①假定待提取包络线的信号为X，则将X W行向量的形式表示为X= [Xi X2…Xw_i Xw]，其中，在此符号"□"为向量表示符号，N表示X的采样点数，Xi表示X中的第1个采样 值，X康示X中的第2个采样值，X 康示X中的第N-1个采样值，X W表示X中的第N个采 样值；
[0007] ②找出X中的所有极大值点和所有极小值点，然后将从X中找出的所有极大值点 按序排列构成一个X的极大值点向量，记为P。，并将从X中找出的所有极小值点按序排列构 成一个X的极小值点向量，记为Pb;
[000引⑨根据所要构建的DCT基的阶数和用于改变所要构建的DCT基的带宽的变化因 子，构建一个DCT基，记为W，其中，W为一个N阶方阵；
[0009] ④从W中提取出行号与X中的每个极大值点的下标一致的每行元素，然后将提取 出的所有行按行号顺序排列构成一个维数为KiXN的矩阵，记为Hi,其中，Ki表示X中的极 大值点的总个数，1《Ki<N ;
[0010] 同样，从W中提取出行号与X中的每个极小值点的下标一致的每行元素，然后将 提取出的所有行按行号顺序排列构成一个维数为KsXN的矩阵，记为&，其中，馬表示X中 的极小值点的总个数，1《K2<N ;
[00川⑥将P。作为观测向量，将H 1作为感知矩阵，利用正交匹配追踪算法恢复出的信号 即为X的上包路线，记为X。；
[0012] 同样，将Pb作为观测向量，将H2作为感知矩阵，利用正交匹配追踪算法恢复出的信 号即为X的下包路线，记为Xb;
[001引⑧构建一个维数为（N-1) XN的差分矩阵，记为D，D中第i行第i列的元素的值为 1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的 元素外的所有元素的值均为0 ;然后根据X。和D获取X。的平滑度，记为H。，同样根据Xb和D 获取Xb的平滑度，记为Hb;
[0014] ⑦判断H点否小于min_H如果是，则令min_&= H并将X。作为X的最佳上包 络线，然后执行步骤⑨，否则，直接执行步骤⑨，其中，min_H。的初始值为无穷大；
[001引⑨判断Hb是否小于min_Hb，如果是，则令min_Hb= Hb，并将Xb作为X的最佳下包 络线，然后执行步骤⑨，否则，直接执行步骤⑨，其中，min_Hb的初始值为无穷大；
[0016] ⑨判断m是否等于md，如果是，则分别输出X的最佳上包络线和X的最佳下包络 线，否则，令m = m+l，然后返回步骤⑨继续执行，其中，m表示用于改变所要构建的DCT基的 带宽的变化因子，m的初始值为1，1《m《md，md表示设定的变化因子最大值，m = m+1中 的"="为赋值符号。
[0017] 所述的步骤②中X的极大值点向量的获取过程为：
[001引②-la、对X求一阶差分，得到X的一阶差分向量，记为dx ;
[0019] ②-2a、根据dx中的每个元素的值，获取一个元素的值为1或0的新向量，记为A,; 对于dx中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值 小于或等于零，则将该元素的值置为0,其中，1《j《N ;
[0020] ②-3a、对A,求一阶差分，得到A ,的一阶差分向量，记为B
[0021] ②-4a、在By中找出值小于零的所有元素；然后将找出的所有元素的下标按序排 列构成一个位置向量，记为化再将化，中的每个元素的值加1，得到新的位置向量，记为 Wly'，Wly'中的任一个元素的值为X中的一个极大值点的下标；
[0022] ②-5a、根据化中的每个元素的值，在X中找出所有极大值点，对于化中的任 一个元素，其值为X中的一个极大值点的下标；然后将从X中找出的所有极大值点按序排列 构成一个X的极大值点向量，记为P。；
[0023] 所述的步骤②中X的极小值点向量的获取过程为：
[0024] ②-化、对X进行取反操作，得到X的反向量，记为X ;然后对X求一阶差分，得到i 的一阶差分向量，记为化；
[0025] ②-2b、根据植中的每个元素的值，获取一个元素的值为1或0的新向量，记为 As ;对于化中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素 的值小于或等于零，则将该元素的值置为0,其中，1《j《N ;
[0026] ②-3b、对A;求一阶差分，得到A;的一阶差分向量，记为B;;
[0027] ②-4b、在中找出值小于零的所有元素；然后将找出的所有元素的下标按序排 列构成一个位置向量，记为Wl;;接着将Wl,中的每个元素的值加1，得到新的位置向量，记 为WV;再对WV进行取反操作，得到WV的反向量，记为WV'，W1,"中的任一个元素 的值为X中的一个极小值点的下标；
[002引②-5b、根据WV'中的每个元素的值，在X中找出所有极小值点，对于WV'中的 任一个元素，其值为X中的一个极小值点的下标；然后将从X中找出的所有极小值点按序排 列构成一个X的极小值点向量，记为Pb。
[0029] 所述的步骤⑨的具体过程为：
[0030] ⑨-1、假定所要构建的DCT基的阶数为N阶，并令m表示用于改变所要构建的DCT 基的带宽的变化因子，其中，m的初始值为1，1《m《md，md表示设定的变化因子最大值；
[0031] ⑨-2、构建一个N阶方阵，记为W，将W中第P行第q列的元素的值记为 W (p，q)，^{P'q)=.y]2/N xc:〇s^(/，().5)x(<y-().5)]，其中，p 和 q 的初始值均为 1， 1《P《N，1《q《N，cos 0为求余弦函数；
[0032] ⑨-3、将步骤⑨-2构建的W作为DCT基。
[0033] 所述的步骤⑥中X的上包络线的获取过程为：
[0034] ⑥-la、令ka表示迭代次数，令P1表示维数为Ki的向量，令B1表示维数为K 1XKi 的矩阵，其中，ka的初始值为1，1《ka《Ki;
[00对⑥-2a、计算H冲的每一列与第ka-1次迭代的残差值r kw之间的相关系数，共得 到N个相关系数；然后将N个相关系数中的最大值对应的Hi中的一列的列号作为P1中的 第ka个元素的值，并将N个相关系数中的最大值对应的Hi中的一列作为B1中的第ka列； 再在Hi中将N个相关系数中的最大值对应的H 1中的一列剔除；
[0036] 其中，H冲的第j列与第ka-1次迭代的残差值Tka4之间的相关系数为H冲的第 j列与第ka-1次迭代的残差值Tkd之间的内积，1《j《N，当ka = 1时取r k,_i= P。；
[0037] ⑥-3a、利用最小二乘法，计算',。_在B1中的每一列上的稀疏映射系数，然后将所 有稀疏映射系数按序构成稀疏映射系数序列，记为
[00測⑥-4a、根据B1和aka，计算第ka次迭代的残差值，记为心町。=r ka-1-Bl X aka;
[0039] ⑥-5a、判断ka是否等于Ki，如果是，则执行步骤⑥-6a，否则，令ka = ka+1，然后 返回步骤⑥-2a继续执行，其中，ka = ka+1中的"="为赋值符号；
[0040] ⑥-Sa、根据P1和ak满取稀疏信号，记为Sp。，Sp。共包含Ki个非零信号值和N-K 1 个零值，Sp,中下标与P1中的第i个元素的值一致的元素的值为非零信号值且等于aka中 的第i个元素的值，Spu中下标与P1中的任一个元素的值不一致的元素的值等于0,其中， 1《i《K1 ;
[0041] ⑥-7a、根据Sp。和W，计算由P。恢复出的信号，记为X。，X。=Sp。X Y，然后将X。作 信号即为X的下包络线，记为Xb;
[0042] 所述的步骤⑥中X的下包络线的获取过程为：
[0043] ⑥-化、令化表示迭代次数，令P2表示维数为K2的向量，令B2表示维数为K2XK2 的矩阵，其中，化的初始值为1，1《化《馬；
[0044] ⑥-2b、计算H,中的每一列与第化-1次迭代的残差值r kb_i之间的相关系数，共得 到N个相关系数；然后将N个相关系数中的最大值对应的&中的一列的列号作为P2中的 第化个元素的值，并将N个相关系数中的最大值对应的中的一列作为B2中的第化列； 再在&中将N个相关系数中的最大值对应的H2中的一列剔除；
[0045] 其中，中的第j列与第化-1次迭代的残差值rkb_i之间的相关系数为中的第 j列与第化-1次迭代的残差值rkb_i之间的内积，1《j《N，当化=1时取r kb_i= P b;
[0046] ⑥-3b、利用最小二乘法，计算心在B2中的每一列上的稀疏映射系数，然后将所 有稀疏映射系数按序构成稀疏映射系数序列，记为akb;
[0047] ⑥-4b、根据B2和3化，计算第化次迭代的残差值，记为rkb，fkb二r kb-1-B2 X a化；
[0048] ⑥-5b、判断化是否等于馬，如果是，则执行步骤⑥-6b，否则，令化=化+1，然后 返回步骤⑥-化继续执行，其中，化=化+1中的"="为赋值符号；
[0049] ⑥-6b、根据P2和akb获取稀疏信号，记为Sp>，Sp>共包含馬个非零信号值和N-K 2 个零值，中下标与P2中的第i个元素的值一致的元素的值为非零信号值且等于akb中 的第i个元素的值，Sa中下标与P2中的任一个元素的值不一致的元素的值等于0,其中， 1《i《K2 ;
[0050] ⑥-7b、根据Spi和W，计算由Pb恢复出的信号，记为Xb，而=SaX'F，然后将Xb作 为X的信号下包络线。
[0化1] 所述的步骤⑧的具体过程为：
[005引⑧-1、构建一个维数为（N-1) XN的差分矩阵，记为D，D中第i行第i列的元素的 值为1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1 列的元素外的所有元素的值均为0,即
【权利要求】
1. 一种基于稀疏复原的信号包络线提取方法，其特征在于包括以下步骤： ① 假定待提取包络线的信号为X，则将X以行向量的形式表示为X= [X1X2…xN_ixN]， 其中，在此符号" □"为向量表示符号，N表示x的采样点数，X1表示x中的第1个采样值，x2 表示X中的第2个采样值，xN_i表示X中的第N-I个采样值，XN表示X中的第N个采样值； ② 找出X中的所有极大值点和所有极小值点，然后将从X中找出的所有极大值点按序 排列构成一个X的极大值点向量，记为Pa，并将从X中找出的所有极小值点按序排列构成一 个X的极小值点向量，记为Pb; ③ 根据所要构建的DCT基的阶数和用于改变所要构建的DCT基的带宽的变化因子，构 建一个DCT基，记为W，其中，W为一个N阶方阵； ④ 从W中提取出行号与X中的每个极大值点的下标一致的每行元素，然后将提取出的 所有行按行号顺序排列构成一个维数为K1XN的矩阵，记为H1，其中，K1表示X中的极大值 点的总个数，1彡Ki〈N; 同样，从W中提取出行号与X中的每个极小值点的下标一致的每行元素，然后将提取 出的所有行按行号顺序排列构成一个维数SK2XN的矩阵，记为H2，其中，K2表示X中的极 小值点的总个数，1彡K2〈N; ⑤ 将Pa作为观测向量，将H1作为感知矩阵，利用正交匹配追踪算法恢复出的信号即为 X的上包路线，记为xa; 同样，将Pb作为观测向量，将H2作为感知矩阵，利用正交匹配追踪算法恢复出的信号即 为X的下包路线，记为Xb; ⑥ 构建一个维数为（N-I)XN的差分矩阵，记为D，D中第i行第i列的元素的值为1，D 中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的元素 外的所有元素的值均为O;然后根据^和D获取X3的平滑度，记为Ha，同样根据&和D获 取xb的平滑度，记为Hb; ⑦ 判断札是否小于min_Ha，如果是，则令min_Ha=Ha，并将xa作为X的最佳上包络线， 然后执行步骤⑧，否则，直接执行步骤⑧，其中，min_Ha的初始值为无穷大； ⑧ 判断Hb是否小于min_Hb，如果是，则令min_Hb=Hb，并将xb作为X的最佳下包络线， 然后执行步骤⑨，否则，直接执行步骤⑨，其中，min_Hb的初始值为无穷大； ⑨ 判断m是否等于md，如果是，则分别输出X的最佳上包络线和X的最佳下包络线，否 贝1J，令m=m+1，然后返回步骤③继续执行，其中，m表示用于改变所要构建的DCT基的带宽 的变化因子，m的初始值为1，1彡m彡md，md表示设定的变化因子最大值，m=m+1中的 "="为赋值符号。
2. 根据权利要求1所述的一种基于稀疏复原的信号包络线提取方法，其特征在于所述 的步骤②中X的极大值点向量的获取过程为： ②-la、对X求一阶差分，得到X的一阶差分向量，记为dx; ②_2a、根据dx中的每个元素的值，获取一个元素的值为1或O的新向量，记为Ax;对 于dx中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小 于或等于零，则将该元素的值置为〇,其中，1彡j彡N; ②-3a、对Ax求一阶差分，得到A-阶差分向量，记为Bx; ②-4a、在匕中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构 成一个位置向量，记为Wlx;再将Wlx中的每个元素的值加I，得到新的位置向量，记为Wlx'，Wlx'中的任一个元素的值为X中的一个极大值点的下标； ②-5a、根据Wlx'中的每个元素的值，在X中找出所有极大值点，对于Wlx'中的任一个 元素，其值为X中的一个极大值点的下标；然后将从X中找出的所有极大值点按序排列构成 一个X的极大值点向量，记为Pa; 所述的步骤②中X的极小值点向量的获取过程为： ②-lb、对X进行取反操作，得到X的反向量，记为i;然后对i求一阶差分，得到X的一 阶差分向量，记为狀； ②-2b、根据dX中的每个元素的值，获取一个元素的值为1或O的新向量，记为Ai;对 于dS中的第j个元素，如果该元素的值大于零，则将该元素的值置为1，如果该元素的值小 于或等于零，则将该元素的值置为〇,其中，1彡j彡N; ②-3b 求一阶差分，得到Ai的一阶差分向量，记为Bi; ②-4b、在Bx中找出值小于零的所有元素；然后将找出的所有元素的下标按序排列构 成一个位置向量，记为Wli;接着将Wli中的每个元素的值加i，得到新的位置向量，记为 Wli 再对WV进行取反操作，得到Wli '的反向量，记为Wli "，Wli "中的任一个元素的 值为X中的一个极小值点的下标； ② -5b、根据Wli "中的每个元素的值，在X中找出所有极小值点，对于Wli "中的任一 个元素，其值为X中的一个极小值点的下标；然后将从X中找出的所有极小值点按序排列构 成一个X的极小值点向量，记为Pb。
3. 根据权利要求1或2所述的一种基于稀疏复原的信号包络线提取方法，其特征在于 所述的步骤③的具体过程为： ③ -1、假定所要构建的DCT基的阶数为N阶，并令m表示用于改变所要构建的DCT基的 带宽的变化因子，其中，m的初始值为1，I<m<md，md表示设定的变化因子最大值； ③-2、构建一个N阶方阵，记为叫，将屯中第p行第q列的元素的值记为
l<p<N，l<q<N，cos()为求余弦函数； ③-3、将步骤③-2构建的屯作为DCT基。
4. 根据权利要求3所述的一种基于稀疏复原的信号包络线提取方法，其特征在于所述 的步骤⑤中X的上包络线的获取过程为： ⑤-la、令ka表示迭代次数，令Pl表示维数为K1的向量，令Bl表示维数为KiXK1的矩 阵，其中，ka的初始值为1，I<ka<K1; ⑤-2a、计算氏中的每一列与第ka-1次迭代的残差值r之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H1中的一列的列号作为Pl中的第ka 个元素的值，并将N个相关系数中的最大值对应的H1中的一列作为Bl中的第ka列；再在 H1中将N个相关系数中的最大值对应的H:中的一列剔除； 其中，H1中的第j列与第ka-1次迭代的残差值rw之间的相关系数为Hi中的第j列 与第ka-1次迭代的残差值!Tkiri之间的内积，1彡j彡N，当ka= 1时取rJ5iri=Pa; ⑤_3a、利用最小二乘法，计算在BI中的每一列上的稀疏映射系数，然后将所有稀 疏映射系数按序构成稀疏映射系数序列，记为aka; ⑤-4a、根据BI和aka，计算第ka次迭代的残差值，记为rka，:Tka=rh-BIXaka; ⑤-5a、判断ka是否等于K1，如果是，则执行步骤⑤_6a，否则，令ka=ka+1，然后返回 步骤⑤_2a继续执行，其中，ka=ka+1中的"="为赋值符号； ⑤-6a、根据Pl和aka获取稀疏信号，记为SPa，SPu共包含&个非零信号值和n-kif 零值，Spii中下标与Pl中的第i个元素的值一致的元素的值为非零信号值且等于中的 第i个元素的值，\中下标与Pl中的任一个元素的值不一致的元素的值等于0,其中， 1彡i彡Kl; ⑤-7a、根据St和W，计算由Pa恢复出的信号，记为Xa，气=S11X中，然后将Xa作为X的信号上包络线； 同样，将Pb作为观测向量，将H2作为感知矩阵，利用正交匹配追踪算法恢复出的信号即 为X的下包络线，记为Xb; 所述的步骤⑤中X的下包络线的获取过程为： ⑤-lb、令kb表示迭代次数，令P2表示维数为K2的向量，令B2表示维数为K2XK2的矩 阵，其中，kb的初始值为1，I<kb<K2; ⑤-2b、计算4中的每一列与第kb-1次迭代的残差值r之间的相关系数，共得到N个相关系数；然后将N个相关系数中的最大值对应的H2中的一列的列号作为P2中的第kb 个元素的值，并将N个相关系数中的最大值对应的H2中的一列作为B2中的第kb列；再在 H2中将N个相关系数中的最大值对应的H2中的一列剔除； 其中，H2中的第j列与第kb-1次迭代的残差值r之间的相关系数为H2中的第j列 与第kb-1次迭代的残差值!Tklrt之间的内积，1彡j彡N，当kb= 1时取rH=Pb; ⑤_3b、利用最小二乘法，计算Fklrf在B2中的每一列上的稀疏映射系数，然后将所有稀 疏映射系数按序构成稀疏映射系数序列，记为akb; ⑤-4b、根据B2和akb，计算第kb次迭代的残差值，记为rkb，rkb=rJdrt-BZXal5b; ⑤-5b、判断kb是否等于K2，如果是，则执行步骤⑤-6b，否则，令kb=kb+1，然后返回 步骤⑤_2b继续执行，其中，kb=kb+1中的"="为赋值符号； ⑤-6b、根据P2和akb获取稀疏信号，记为SP6，SP6共包含1(2个非零信号值和N-K2个 零值，Sn中下标与P2中的第i个元素的值一致的元素的值为非零信号值且等于akb中的 第i个元素的值，sP6中下标与P2中的任一个元素的值不一致的元素的值等于0,其中， 1彡i彡K2 ; ⑤-7b、根据Sp,和W，计算由Pb恢复出的信号，记为xb，弋=SP6X Y，然后将Xb作为x的 信号下包络线。
5.根据权利要求4所述的一种基于稀疏复原的信号包络线提取方法，其特征在于所述 的步骤⑥的具体过程为： ⑥-1、构建一个维数为（N-I)XN的差分矩阵，记为D，D中第i行第i列的元素的值为 1，D中第i行第i+1列的元素的值为-1，D中除第i行第i列的元素和第i行第i+1列的
⑥-2、根据xjPD计算的差分向量，记为dxa，dxa=XaXD;同样，根据xjPD计算Xb 的差分向量，记为dxb，dxb=xbXD; ⑥-3、计算dxa中的所有元素的值的平方和，记为Ha，然后将Ha作为xa的平滑度；同样， 计算dxb中的所有元素的值的平方和，记为Hb，然后将Hb作为Xb的平滑度。
【文档编号】G06F17/50GK104504181SQ201410751425
【公开日】2015年4月8日 申请日期:2014年12月10日 优先权日:2014年12月10日
【发明者】徐静妹, 叶庆卫, 周宇, 王晓东 申请人:宁波大学