一种基于多尺度线调频基的稀疏信号分解方法✳

彭富强,于德介,刘 坚

(湖南大学汽车车身先进设计制造国家重点实验室,湖南长沙 410082)

引 言

近年来,国内外学者通常采用短时 Fourier变换、小波变换、Gabor变换等时频分析方法来获得信号的时频分布,但这些时频分析方法均存在一定局限。短时Fourier变换中谱窗函数的引用会降低局部谱的分辩率[1,2]。小波变换使用平行分割的等面积时频窗[3],只适用于分析具有固定比例带宽的非平稳信号,然而实际情况下信号频率成分通常是呈曲线变化的。Gabor变换等二次型时频分析则由于交叉干扰项限制了它对多分量信号分析的应用。EMD方法在理论上还存在诸多问题有待进一步研究,如EMD方法中的过包络、欠包络、模态混淆和端点效应等问题[4～6]。为此亟待新的时频分析方法来解决现有时频分析方法的一些固有缺陷。

Mallat等人于 1993年提出了稀疏信号分解的概念[7],在过完备库上根据信号本身的变换特点,自适应的选择基函数对信号进行投影分解,得到信号的稀疏表示,其过程称为信号的稀疏分解。 Emmanuel JCandès等人于 2007年提出了线调频小波路径追踪算法[8],该算法通过对线调频小波图中的线调频小波原子进行连接,自适应的获得频率呈曲线变化的信号分量[9～12]。 Emmanuel JCandès等人提出的方法目前仅应用在地震引力波的分析中,且只能对单分量信号进行分析。文章将稀疏信号分解与线调频小波路径追踪算法进行结合,提出一种基于多尺度线调频基的稀疏信号分解方法。通过数学推导,发现了投影系数中包含了分解信号的幅值和初始相位信息,并进一步利用投影系数和基函数推导出了分解信号的数学形式,得到分量信号的稀疏表示,成功地将稀疏信号方法的引进来,解决了线调频小波路径追踪算法只适应用于对单分量信号进行分解的局限。该方法在单次分解中采用线调频小波路径追踪算法获得与分析信号具有最大相关系数的信号分量,然后采用稀疏信号分解的概念,将该分量从分析信号中除去,再进一步分解,这样逐次将各信号分量分解出来。该分解方法可以自适应地根据信号频率的变化复杂度采用相应尺度的线调频基函数进行投影分解,在频率变化复杂的地方采用时间支撑区较短的线调频基函数进行投影分解,反之则采用时间支撑区较长的线调频基函数进行分解,所以该方法非常适合于分解频率呈曲线变化的多分量非平稳信号。为了验证该方法的有效性,将该方法与EMD方法进行了对比,结果表明该方法对多分量非平稳信号的分解效果优于 EMD方法,适用于多分量非平稳信号的分解。

1 多尺度线调频基函数库

根据信号分析理论,任意信号f(t)可以展开为一组基函数的线性组合[13],即

如果该组基函数为正交基,则可用内积计算它们的展开系数 ,即

式中an的大小反映了f(t)与基函数的相似程度[13]。本文所用的基函数库为多尺度的线性调频基函数

要求‖ha,b,I‖ = 1,则

式中D为本文定义的基函数库;ha,b,I(t)为多尺度线调频基函数;I为动态时间支撑区,I=[kN2-j～(k+ 1)N2-j],j为分析尺度系数,j= 0,1,…,log2(N-7),N为采样长度,k= 0,1,…,2j-1;Ka,b,I为归一化系数,使得‖ha,b,I‖ =1;a为频率偏置系数,取采样频率的一半;2b为频率斜率,根据信号频率变化快慢设定。a,b与尺度系数j有关,根据采样定理,a+2bt应该小于fs/2;1I(t)为矩形窗函数 ,当t∈I时为 1,当t∉I时为 0。

式(3)定义的多尺度线性调频基函数在动态时间支撑区内的瞬时频率为a+2bt,适合在小的动态时间支撑区内逐段拟合频率呈曲线变化的频率成分。动态时间支持区长度太小将影响分析结果,误差较大,所以动态时间支持区长度必须大于32,尺度系数j取0,1,… ,log2(N-7)。本文方法要求采样长度为2的整数次方。通过计算可获得每个支撑区I内的最大投影系数和对应的线调频基函数。线调频基函数的多尺度特性使得它具有了动态匹配分析信号的特性,而不像传统的稀疏信号分解一样,在整个信号分析时间段内采用单一尺度基函数对信号进行分解,基函数中包含的频率斜率信息则使得其适合分析频率呈曲线变化的非平稳信号,克服了小波变换存在的具有多尺度特性,但不具有频率斜率信息的问题。

2 基于多尺度线调频基的稀疏信号分解

2.1 最大投影系数求解与分解信号定义

根据稀疏信号分解理论,基于多尺度线调频基的稀疏信号分解方法通过从基函数库D中挑选一组基函数来对f(t)进行分解,f(t)在该组基函数的动态时间支撑区I上具有最大的投影系数,且该组基函数的支撑区应该覆盖整个分析信号,不重叠。支撑区I内的最大投影系数UI计算公式为

假设分析信号为

有

周期函数在一个周期内的积分为0,所以

当θ(t)-at-bt2=0时上式成立。即投影系数UI包含了分解信号的幅值信息和初始相位信息,对应的基函数则包含了分解信号的频偏信息和频率斜率信息。根据以上推导文章定义cI(t)为在动态时间支撑区I内最大投影系数代表的分解信号,则

在动态时间支撑区I内

式中rI(t)为分解的残余信号,因为rI(t)与cI(t)对应的基函数正交[7],所以

6)关联度：γi=1/N×∑εi(k)，式中：γi为比较数列Xi对参考数列X0的关联度；N为评价指标个数。

投影系数的最大化保证了残余信号能量的最小。通过式(6)～ (8)可以计算出每个动态时间支撑区内的最大投影系数和相应的基函数。从式(6)可以看出,当多个信号分量幅值相同,但频率成分不同时将会具有几乎相同的投影系数,为了在此种情况下区分不同频率成分的信号,文章方法将同时保留投影系数相差不大的基函数参与分解。

2.2 动态时间支撑区连接算法

根据线调频小波路径追踪算法,要求找到一种动态时间支撑区连接方法,在该连接方法下满足在整个分析时间内使分解信号的总能量最大,即

式中n代表第n次分解,∏n覆盖整个分析时间段,不能重叠。∏n的连接方法应保证在本次分解中分解信号总能量最大,线调频小波路径追踪算法提出的连接算法如下[8]:

(1)初始化:d(i)为第i个动态分析段之前分解信号的总能量,pre(i)为连接到第i个动态时间支撑区的前置动态时间支撑区序号,e(i)为第i个动态时间支撑区最大投影系数对应的信号能量,初始化时,置d(i)=0,pre(i)=0;

(2)对于动态时间支撑区集合{Ii,i∈Z}中的每一个元素Ii,查找出与其相邻的所有下一个动态时间支撑区集合{Ij}:稀疏信号分解方法根据投影系数大小进行分解,投影系数大的先分解,小的后分解,但当多个分量具有相同的幅值时,其具有相同的投影系数,出现交叉分解的现象,为解决等振幅分解的问题,文章在支撑区连接算法中引入了保留系数W,即保留投影系数相差不大的基函数参与分解,如果

有

W取1,只有当分解信号频率出现交叉解调的现象时,才将W的值逐次减小后重新分解,减小的步长取 0.1,直到分解信号光滑连续。

为了保证分段调频信号频率变化的连续性和光滑性,支撑区连接算法还应遵循以下原则:

2.优先选择相邻动态时间支撑区频率变化较平坦的频率成分。

2.3 稀疏信号分解

通过上述搜寻算法可以保证获得总信号能量最大的分解信号,参数pre(i)中包含了连接路径,原始信号与分解信号的关系,即

式中rn为第n次分解的残余信号,cn为分解信号。基于多尺度线调频基的稀疏信号分解方法是一个逐次分解的过程,本次分解的残余信号可以进一步地分解,直到残余信号能量小于一定的阈值便停止分解。随着n的增大,残余信号rn的能量迅速衰减到0。分解流程如图 1所示。

图1 分解流程图

3 仿真分析与噪声分析

3.1 多分量非平稳仿真信号分析

为了说明本文方法的有效性,在 MATLAB上对仿真信号进行了分析。

仿真信号为

式中 e(t)为高斯白噪声,信噪比为 5 d B,采样时长5.115 s,采样频率为 200 Hz,采样点数为 1 024点。仿真信号包含 3个频率分量,其中后两个频率分量的幅值相等,每个分量的频率都呈正弦变化。图2为原始信号图形。

图2 原始仿真信号

取尺度系数0～ 5。分析点数太少将影响分析结果,误差较大,最小分析点数必须大于 32。根据采样率,信号的最大分析频率不超过100 Hz,所以分析频偏范围取 0～ 100 Hz,搜寻频偏分辨率为 1 Hz,信号的频率斜率范围为-473～ 473 Hz/s,所以搜寻频率斜率范围取-500～500 Hz/s,搜寻分辨率为1 Hz/s。第一次分解时,W取 1可以分解出幅值为 1的分量,没有出现等幅值交叉解调的现象,而后两个分量因为具有相同的幅值,所以W取1时出现了交叉解调的现象,如图3所示。图中虚线为仿真信号的频率变化曲线,实线为分解信号频率曲线,因此将W值减小取 0.9,则可以顺利分解出后两个等幅信号分量。根据信噪比计算,噪声信号能量占信号总能量的38.69%,所以当残余信号能量小于仿真信号能量的38.69%时停止分解。采用本文方法,经过3次分解,残余信号能量只占分析信号能量的24.45%,算法停止分解。

图3 交叉解调现象

图4 仿真信号与分解信号时频图

图4为仿真信号与分解信号时频图。图中虚线为仿真信号的频率变化曲线,实线为分解出的分解信号频率曲线,分解顺序如箭头所示。可以看到经过3次分解可以解调出仿真信号的3个信号分量,分解信号很好地匹配了仿真信号,频率拟合较光滑,没有出现同幅值分量的交叉解调现象。从图4的曲线1可以看到,在频率变化较平坦的时间段,该方法采用了动态时间支撑区较大的基函数进行投影,而在频率变化复杂的时间段则会自动采用动态时间支撑区较小的基函数进行投影。分解出的 3个信号分量如图5,6和7所示,其中第一个分量幅值近似为1,后两个分量幅值近似为 0.5,较好地匹配了原始信号分量。

图5 分解信号分量 1波形图cos{2×π× 3× [10t+ 3sin(2×π× 0.25×t)]}

图6 分解信号分量 2波形图0.5cos{2×π× 4× [10t+ 3sin(2×π× 0.25× t)]}

图7 分解信号分量 3波形图0.5cos{2×π× 2× [10t+ 3sin(2×π× 0.25× t)]}

将该仿真信号通过 EMD方法进行分解,分解出10个 IMF分量,对各个 IMF分量进行 Hilbert变换,获得瞬时频率,前两个 IMF分量的时频图如图 8和 9所示,虚线为仿真信号的频率变化曲线,实线为分解出的IMF分量频率变化曲线,从图中可以看出,EMD方法分解效果较差,没有很好地将 3个信号分量从噪声中提取出来。图 10为其余 IMF分量,依然没有分解出仿真信号的频率成分,对比图 4与图8～10,从分解的效果看,本文提出的方法明显优于 EMD分解方法,从而验证了本文方法的有效性。

图8 EMD分解第一个 IM F分量时频图

图9 EMD分解第二个 IMF分量时频图

图10 第 3～ 10个 IMF分量时频图

3.2 噪声分析

为验证文章方法的抗噪性能,对在不同信噪比下的加噪线调频信号和纯噪声信号进行分析,分别将其对基函数库中元素进行投影,获得最大投影系数。

加噪线调频信号为

采样率为200 Hz,采样点数256。将不同信噪比下的加噪线调频信号和纯噪声信号对基函数投影,其最大投影系数如图 11所示。从图中可以看出,当信噪比大于-8 d B时加噪线调频信号的最大投影系数相比纯噪声最大投影系数占优势,仍能很好地从噪声信号中提取线调频信号,证明了文章方法具有很好的抗噪性能,其最大投影系数对应基函数的频偏和频率斜率如表 1所示,表中参数与仿真线调频信号参数很好地吻合,证明了文章方法的有效性。

图11 仿真线调频信号与纯噪声信号最大投影系数

4 齿轮箱振动信号分析

为验证文章方法的有效性,对多级齿轮箱振动信号进行分解。因为试验条件有限,现对一正常单级传动齿轮箱在不同转速下测量两组数据,将两组数据叠加模拟多级齿轮箱振动信号,齿轮箱输入轴齿轮齿数为 55,输出轴齿轮齿数为 75。当齿轮箱正常时,通常多级齿轮箱振动信号以各齿轮副啮合频率分量为主,采集齿轮箱振动加速度信号,采样频率为4 096 Hz,采样时长为1.999 8 s。叠加振动信号时域波形图如图12所示。

表1 最大投影系数对应基函数频偏与频率斜率参数

图12 齿轮箱叠加振动信号时域波形图

对模拟多级齿轮箱振动信号进行基于多尺度线调频基的稀疏信号分解。根据采样率,分析频偏范围取 0～ 2 048 Hz,搜寻频偏分辨率为 1 Hz,由于齿轮箱的加减速需要一定时间,转速不可能在短时间内过度剧烈变化,所以啮合频率斜率不可能超过-100～ 100 Hz/s,因此尺度系数 0～ 4,搜寻频率斜率范围取-100～100 Hz/s,搜寻分辨率为1 Hz/s。图 13中虚线为测量所得的2个啮合频率分量,经 2次分解,得到 2个信号分量如图13中实线所示,分解顺序如箭头所示,文章方法很好地从振动信号中提取了 2个齿轮啮合频率分量,验证了该方法对多分量非平稳信号的分解的有效性。

图13 分解信号与啮合频率分量时频图

5 结 论

本文提出了一种基于多尺度线调频基的稀疏信号分解方法,根据最大投影系数和线调频基函数导出了分解信号的理论形式,证明了最大投影系数包含了分解信号的幅值和相位信息,并对时间支撑区的连接算法进行了优化,引入了保留系数W和偏差阈值X,解决了稀疏信号分解中的等振幅交叉解调问题。基于多尺度线调频基的稀疏信号分解方法具有柔性的匹配性能,良好的时频聚集性,正交分解等特点,克服了以往时频分析方法的诸多固有缺陷,非常适合于多分量非平稳信号的分解。对仿真信号的分析结果表明,与 EMD方法比较,基于多尺度线调频基的稀疏信号分解方法具有更高的频率匹配特性,是一种新的有效的时频分析方法。但该方法目前还不够完善,还存在很多值得研究的问题,如新的多尺度划分方法,端点拟合精度的提高,计算方法的优化等等,有待近一步深入研究。

[1] Bizarro J P S,Figueiredo A C A.Time-frequency analysis of fusion plasma signals beyond the shorttime Fourier transform paradigm:An overview[J].Fusion Engineering and Design,2008,83(2-3):350—353.

[2] 张贤达,保铮.非平稳信号分析与处理 [M].北京:国防工业出版社,1998:20— 26.

[3] Benko U,JuˇricˊicD.Frequency analysis of noisy shorttime stationary signals using filter-diagonalization[J].Signal Processing,2008,88(7):1 733—1 746.

[4] Qin S R,Zhong Y M.A new envelope algorithm of Hilbert-Huang transform[J].Mechanical Systems and Signal Processing,2006,20(8):1 941—1 952.

[5] Huang N E,Wu M-L C,Long S R.A confidence limit for the empirical mode decomposition and Hilbert spectral analysis[J].Proc.R.Soc.Lond.A,2003,459:2 317—2 345.

[6] Deng Yongjun,Wang Wei.Boundary processing technique in EMD method and Hilbert transform[J].Chinese Science Bulletin,2001,46(11):257—263.

[7] Mallat S,Zhang Z.Matching pursuit with time-frequency dictionaries[J].Signal Processing,1993,41(12):3 397—3 415.

[8] Candès E J,Charlton P R,Helgason H.Detecting highly oscillatory signals by chirplet path pursuit[J].Applied and Computational Harmonic Analysis,2008,24(1):14—40.

[9] Spanos P D,Giaralis A,Politis N P.Time-f requency representation of earthquakeaccelerograms and inelastic structural response records using the adaptive chirplet decomposition and empirical mode decomposition[J].Soil Dynamics and Earthquake Engineering,2007,27(7):675—689.

[10]Dugnol B,Fernández C,Galiano G,et al.On a chirplet transform-based method applied to separating and counting wolf howls[J].Signal Processing,2008,88(7):1 817—1 826.

[11]Jie Cui,Wong W.Investigation of short-term changes in visual evoked potentials with windowed adaptive chirplet transform[J].Biomedical Engineering,2008,55(4):1 449—1 454.

[12]Greenberg J M,Zhisong Wang,Jian Li.New approaches for chirplet approximation[J].Signal Processing,2007,55(2):734—741.

[13]Qian S,Chen D.Joint Time-Frequency Analysis[M].Englewood Cliffs,Prentice Hall,1996:124—130.