多成分信号的时频谱获取方法及装置
技术领域
本发明涉及信号处理技术领域,尤其涉及一种多成分信号的时频谱获取方法及装置。
背景技术
时频谱是分析和处理非平稳信号的有力工具,它将信号表示为时间和频率的联合函数,能够清楚描述信号频率随时间变化的分布关系。时频分析思想肇始于二十世纪四十年代,1946年提出的Gabor变换(Gabor D.Theory of communication[J].J.IEE,1946,93:429-457)奠定了在时间和频率联合域内分析信号的理论基础。1947年Potter等提出了简单而实用的短时傅里叶变换(STFT)。然而,由于不确定原理的限制,STFT不能兼顾频率和时间分辨率的需求。
1932年物理学家Wigner在量子力学中提出了著名的Wigner分布(Wigner E P.Onthe Quantum Correction for Thermodynamic Equilibrium[J].Physical Review,1931,40(40):749-759),1948年,J.Ville将Wigner分布引入到信号处理领域,从而发展成为后来最具有代表性的一种时频分析方法,即Wigner-Ville分布(Wigner-Ville Distribution,WVD)。Wigner-Ville分布是一种二次型时频表示方法,它能够满足大部分所希望的数学性质,如实值性、对称性、能量守恒、时频边缘特性、时频移位特性等,是描述信号时频分布的有力工具(Claasen T A C M.The Wigner distribution-A tool for time-frequencysignal analysis[J].Philips J Res,1980,35(4-5):276-300)。虽然Wigner-Ville分布的时频聚集度较高,但是Wigner-Ville分布不能保证非负性,而且对于多分量,会产生严重的交叉干扰,使其谱分布难以解释,严重限制了它的广泛应用。
WVD的性质表明,交叉项是实的,混杂于自项成分之间,且其幅度是自项成分的两倍。许多学者围绕如何克服交叉项干扰的弱点,进行了大量研究,设计出理想的核函数来消除交叉项的影响,如伪Wigner-Ville分布、平滑伪Wigner-Ville分布、Choi-William分布(CWD)、锥形核分布(CKD)等,这些分布统称为Cohen类时频分布(邹红星,戴琼海,李衍达,等.不含交叉项干扰且具有WVD聚集度的时频分布之不存在性[J].中国科学:,2001,31(4):348-354)。对于多分量信号,邹红星等证明,既能消除交叉项的干扰,又能保持时频聚集度最优的双线性方法是不存在的(Cohen L.Time-frequency distributions-a review[J].Proceedings of the IEEE,1989,77(7):941-981)。由此可见,Cohen类时频分布虽然能减弱交叉项的干扰,但同时必定降低Wigner-Ville分布的时频聚集度。1976年Kodera K提出时频重排的方法,通过对信号进行重排以提高时频聚集度(Kodera K,Villedary C D,Gendrin R.A new method for the numerical analysis of non-stationary signals[J].Physics of the Earth&Planetary Interiors,1976,12(2):142-150)。这种方法虽然能够提供更高的时频聚集度,但是并没有完全消除交叉项,同时存在重排振荡。
由Wigner-Ville分布的性质可知,对于单分量信号,WVD具有最优的时频聚集度,而且没有交叉项的干扰。也就是说,交叉项干扰只发生在多分量信号身上。因此,将信号分解和WVD相结合,即先将信号分解成多个单分量信号之和,然后求各单分量信号WVD之和,能得到没有交叉项的性能更优的时频谱。传统的信号分解方法是预先给定一组基,然后把信号表示成它们的线性组合,组合系数就是信号在基上的投影(Flandrin P.Time-frequency/time Scale Analysis[J].Academic Press Inc San Diego Ca,1999)。如熟知的傅立叶分析采用L2(R)空间的正交谐波基,短时傅立叶变换采用局域L2()的正交谐波基,小波变换采用一组小波基。Pachori等(Pachori R B,Sircar P.A new technique toreduce cross terms in the Wigner distribution[J].Digital Signal Processing,2007,17(2):466-474)提出采用Fourier-Besse展开先将多分量信号分解成单分量信号,再对各分量信号分别计算WVD后合并,以达到去除交叉项的目的。Mirela等(Bianu M,IsarA.The reduction of interference terms in the time-frequency plane[C]//International Symposium on Signals,Circuits and Systems.2003:461-464vol.2)通过对多分量信号进行Gabor展开,并对各分量信号做WVD后再合成的方法避免WVD的交叉项。对于复杂信号,采用给定的基去逼近信号有可能将一个连续的时频谱变成多个间断的时频谱之和,降低信号的物理意义。
Hilbert-Huang变换(HHT)是近些年发展起来的处理非线性非平稳信号的自适应时频分析方法。它先对信号进行经验模态分解(Empirical Mode Decomposition,EMD)分解,然后借助Hilbert变换引入瞬时频率,获得信号在时频平面上的能量分布,即Hilbert谱(Huang,N.E.,Z.Shen,and S.R.Long,M.C.Wu,E.H.Shih,Q.Zheng,C.C.Tung,and H.H.Liu,1998:The empirical mode decomposition method and the Hilbert spectrum fornon-stationary time series analysis,Proc.Roy.Soc.London,A454,903-995)(Huang,N.E.,Z.Shen,R.S.Long,1999:A new view of nonlinear water waves –the Hilbertspectrum,Ann.Rev.Fluid Mech.,31,417-457)。EMD是HHT的核心,可以在不需要知道任何先验知识的情况下,依据输入信号自身的特点,自适应地将非线性非平稳信号分解成若干个具有不同特征时间尺度的内蕴模态函数(Intrinsic mode function,IMF)之和。EMD在机械故障诊断(Yang B.,C.S.Suh,2004:Interpretation of crack-induced rotor non-linear response using instantaneous frequency,Mechanical Systems and SignalProcessing 18(3):491–513)、模态辨识、生物医学(Echeverria J.C.,J.A.Crowea,Woolfson,et.al.2001:Application of empirical mode decomposition to heart ratevariability analysis,Medical&Biological Engineering&Computing,39 471–479)、图像处理等领域取得了成功应用。HHT方法利用Hilbert变换和差分方法计算IMF的瞬时频率,当IMF残留有噪声,或存在很小的干扰时,或有多个成分叠加时,计算得到的瞬时频率可能失真。Xu et.al.(Xu C.H.,J.F.Liu,G.M.Chen,J.Xie,2010:Application of EMD and WVDto feature extraction from vibration signal of reciprocating pump waves,Journal of China University of Petroleum,34(3):99-103)和Su L.et.al.(Su L.,H.P.Nan,X.Y.Yu,L.H.Wu,J.Wang,2012:Analysis of hydro turbine vibration signalsbased on empirical mode decomposition and Wigner-Ville distribution,ActaJournal of Hydroelectric Engineering,31(2):240:245)分别提出将EMD和EEMD方法与WVD结合计算信号的时频分布,即抑制了WVD对于多分量非平稳信号存在固有交叉项干扰的缺陷,又发挥了WVD描述信号的时变特征性能。
EMD方法的一个主要缺陷是模态混叠,表示在一个IMF中包含差异极大的特征时间尺度,或是相近的特征时间尺度分布在不同的IMF中(Huang,N.E.,Z.Shen,and S.R.Long,M.C.Wu,E.H.Shih,Q.Zheng,C.C.Tung,and H.H.Liu,1998:The empirical modedecomposition method and the Hilbert spectrum for non-stationary time seriesanalysis,Proc.Roy.Soc.London,A454,903-995)(Wu,Z.and N.E.Huang,2004:A study ofthe characteristics of white noise using the empirical mode decompositionmethod,Proc.Roy.Soc.London,A.460,1597-1611)。若一个IMF分量包含多个特征时间尺度,则不满足单分量的条件,计算WVD将发生交叉项。
文献(Peng Z K,Tse P W,Chu F L.An improved Hilbert–Huang transform andits application in vibration signal analysis[J].Journal of Sound&Vibration,2005,286(1–2):187-205)的分析结果表明,EMD还有两个缺点可能导致交叉项:其一是第1个IMF的频段比较宽,除了包含高频成分外,同时也包含相邻的多个低频成分;其二是低能量成分信号可能附加在相邻的高能量IMF中。与模态混叠一样,上述缺点破坏了IMF是单分量信号的条件。在计算各IMF分量信号的WVD时,会产生令人烦恼的交叉项。
发明内容
本发明提供一种多成分信号的时频谱获取方法及装置,以去除多成分信号的Wigner-Ville时频谱分布中的交叉项干扰。
本发明实施例提供一种多成分信号的时频谱获取方法,包括:将多成分信号分离为多个子带信号;计算所述子带信号的Wigner-Ville分布;计算各所述子带信号的Wigner-Ville分布的和,得到所述多成分信号的时频谱。
一个实施例中,将多成分信号分离为多个子带信号,包括:计算所述多成分信号的Wigner-Ville分布及其模的频率边缘;利用所述多成分信号的Wigner-Ville分布及其模的频率边缘确定所述多成分信号的Wigner-Ville分布中交叉项的位置;基于所述多成分信号的Wigner-Ville分布中交叉项的位置,通过低通滤波将所述多成分信号分离为所述多个子带信号。
一个实施例中,利用所述多成分信号的Wigner-Ville分布及其模的频率边缘确定所述多成分信号的Wigner-Ville分布中交叉项的位置,包括:确定所述多成分信号的Wigner-Ville分布的模的频率边缘的极大值点;在所述极大值点处,计算所述多成分信号的Wigner-Ville分布的频率边缘和所述多成分信号的Wigner-Ville分布的模的频率边缘的差值;判断所述差值的模是否大于设定阈值,若是,将所述极大值点作为所述多成分信号的Wigner-Ville分布中交叉项的位置。
一个实施例中,基于所述多成分信号的Wigner-Ville分布中交叉项的位置,通过低通滤波将所述多成分信号分离为所述多个子带信号,包括:以所述多成分信号的Wigner-Ville分布中交叉项的位置为截止频率,对所述多成分信号进行低通滤波,得到第一子带信号,其中,所述多个子带信号包括所述第一子带信号。
一个实施例中,基于所述多成分信号的Wigner-Ville分布中交叉项的位置,通过低通滤波将所述多成分信号分离为所述多个子带信号,还包括:利用所述多成分信号减去所述第一子带信号,得到剩余信号;利用所述剩余信号的Wigner-Ville分布及其模的频率边缘确定所述剩余信号的Wigner-Ville分布中交叉项的位置;以所述剩余信号的Wigner-Ville分布中交叉项的位置为截止频率,对所述剩余信号进行低通滤波,得到第二子带信号,其中,所述多个子带信号还包括所述第二子带信号。
一个实施例中,将多成分信号分离为多个子带信号之前,还包括:利用窗极值经验模态分解对原始信号进行分解得到内蕴模态函数IMF分量信号;根据所述IMF分量信号得到所述多成分信号。
一个实施例中,根据所述IMF分量信号得到所述多成分信号,包括:对所述IMF分量信号进行小波包滤波,得到所述多成分信号。
本发明实施例还提供一种多成分信号的时频谱获取装置,包括:子带分离单元,用于:将多成分信号分离为多个子带信号;WVD计算单元,用于:计算所述子带信号的Wigner-Ville分布;时频谱计算单元,用于:计算各所述子带信号的Wigner-Ville分布的和,得到所述多成分信号的时频谱。
本发明实施例还提供一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现上述各实施例所述方法的步骤。
本发明实施例还提供一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述程序时实现上述各实施例所述方法的步骤。
本发明实施例的多成分信号的时频谱获取方法、装置、存储介质及计算机设备,通过先将包含多分量成分的多成分信号分离为仅包含但分量成分的子带信号,再对各子带信号的Wigner-Ville分布求和得到该多成分信号的时频谱,能够消除直接计算多成分信号的Wigner-Ville分布(时频分布)中产生的交叉项的干扰,从而能够解决谱分布难以分析的问题,并能够提高多成分信号时频分析的精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。在附图中:
图1是本发明一实施例中泄漏电流信号及其利用WE-EMD方法得到的IMF分量信号的曲线图。
图2(a)和图2(b)分别是图1中第1个IMF分量信号和第2个IMF分量信号的功率谱。
图3(a)和图3(b)分别是图1中第1个IMF分量信号的WVD和Hilbert谱。
图4是本发明一实施例的多成分信号的时频谱获取方法的流程示意图。
图5是本发明一实施例中将多成分信号分离为多个子带信号的方法流程示意图。
图6是本发明一实施例中利用多成分信号的Wigner-Ville分布及其模的频率边缘确定交叉项位置的方法流程示意图。
图7是本发明一实施例中基于多成分信号的Wigner-Ville分布中交叉项的位置通过低通滤波将多成分信号分离为子带信号的方法流程示意图。
图8是本发明另一实施例的多成分信号的时频谱获取方法的流程示意图。
图9是本发明一实施例中多成分信号的时频谱获取方法的流程示意图。
图10是利用本发明一实施例的方法得到的WVD及其模的频率边缘的曲线图。
图11是利用本发明一实施例的方法得到的图1中的第1个IMF分量及其子带分离信号的曲线图。
图12(a)和图12(b)分别是图1中第1个IMF分量信号的子带分量信号的WVD之和和全部IMF分量信号的子带分量信号的WVD之和。
图13(a)至图13(f)依次是图1中泄漏电流信号的WVD、伪WVD、重排WVD、连续小波变换、基于EMD的WVD及基于EEMD的WVD的时频谱。
图14是本发明一实施例的多成分信号的时频谱获取装置的结构示意图。
图15是本发明一实施例的计算机设备的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白,下面结合附图对本发明实施例做进一步详细说明。在此,本发明的示意性实施例及其说明用于解释本发明,但并不作为对本发明的限定。
首先,以利用WE-EMD方法(中国发明专利申请,公开号:CN 106844293 A)得到的IMF分量信号为例,说明多成分信号在计算WVD时会产生交叉项的现象。
图1是本发明一实施例中实测的泄漏电流信号及其利用WE-EMD方法得到的IMF分量信号的曲线图。如图1所示,从上之下,第1个是原始的泄漏电流信号,第2个~第6个是IMF分量信号,第7个是趋势项。图2(a)和图2(b)分别是图1中第1个IMF分量信号和第2个IMF分量信号的功率谱。从图2(b)可以看出,第2个IMF分量信号是频率成分为50Hz的单分量信号,是泄漏电流的基频成分。与之相比,如图2(a)所示,第1个IMF分量信号成分要复杂得多,既包括能量最大的150Hz的3倍基频成分,同时也包括与之相邻的50Hz基频成分;此外,能量较低的250Hz的5倍基频成分附加在相邻的能量较高的150Hz IMF中,这里刚好又是第1个IMF。图3(a)和图3(b)分别是图1中第1个IMF分量信号的WVD和Hilbert谱。从图3(a)可以看出,在150Hz和50Hz自项中间产生了100Hz的交叉项,在150Hz和250Hz自项中间产生了200Hz的交叉项,50Hz和250Hz自项之间产生的150Hz交叉项叠加在自项150Hz上,由于能量较低,显示不明显。从图3(b)可以看出,对于多分量信号,Hilbert的瞬时频率在几个频率之间跳跃,失去了物理意义。
由此可见,窗极值经验模态分解(WE-EMD)中碰会到的两个技术问题,它们破坏了IMF是单分量信号的条件,在计算各IMF分量信号的WVD时,会产生交叉项:(1)第1个IMF的频段比较宽,除了包含高频成分外,同时也包含相邻的1个或多个低频成分;(2)低能量成分信号可能附加在相邻的高能量IMF中。实际上,对于包含多分量成分的信号,即多成分信号,在计算Wigner-Ville分布(WVD)时均易出现交叉项。
为了解决多分量信号(例如,多分量的IMF分量信号)的Wigner-Ville分布中的交叉项干扰问题,本发明实施例提供了一种多成分信号的时频谱获取方法。
图4是本发明一实施例的多成分信号的时频谱获取方法的流程示意图。如图4所示,本实施例的多成分信号的时频谱获取方法,可包括:
步骤S110:将多成分信号分离为多个子带信号;
步骤S120:计算所述子带信号的Wigner-Ville分布;
步骤S130:计算各所述子带信号的Wigner-Ville分布的和,得到所述多成分信号的时频谱。
在上述步骤S110中,可以采用多种不同方法将多成分信号分离为多个子带信号之和,例如,离散小波、滤波器组和该多成分信号的Wigner-Ville分布模的频率边缘方法等。
本实施例中,通过先将包含多分量成分的多成分信号分离为仅包含但分量成分的子带信号,再对各子带信号的Wigner-Ville分布求和得到该多成分信号的时频谱,能够消除直接计算多成分信号的Wigner-Ville分布(时频分布)中产生的交叉项的干扰,从而能够解决谱分布难以分析的问题,并能够提高多成分信号时频分析的精度。
由WVD的频率边缘性质可知,WVD的自项是恒正的,交叉项是振荡的,WVD沿着时间轴的积分等于该信号在频率ω处的瞬时能量。于是发明人推论,若频率ω在整个时间内没有出现交叉项,则在频率ω处,WVD和其模沿着时间轴的积分相等,而且等于它的瞬时能量。因此,发明人发现可以根据多成分信号的WVD和多成分信号的WVD模的频率边缘的差别分辨交叉项的位置,即在哪些频率上出现了交叉项,并通过依次低通滤波的方式将多成分信号转换成单分量的子带信号。
图5是本发明一实施例中将多成分信号分离为多个子带信号的方法流程示意图。如图5所示,在上述步骤S110中,将多成分信号分离为多个子带信号的方法,可包括:
步骤S111:计算所述多成分信号的Wigner-Ville分布及其模的频率边缘;
步骤S112:利用所述多成分信号的Wigner-Ville分布及其模的频率边缘确定所述多成分信号的Wigner-Ville分布中交叉项的位置;
步骤S113:基于所述多成分信号的Wigner-Ville分布中交叉项的位置,通过低通滤波将所述多成分信号分离为所述多个子带信号。
在上述步骤S111中,本领域技术人员了解,频率边缘可指对信号的时频分布在时间维度求积分或求和,所得到的关于频率且类似于功率谱的结果。
在上述步骤S112中,由于发明人发现,在交叉项位置处,多成分信号的Wigner-Ville分布的频率边缘和多成分信号的Wigner-Ville分布的模的频率边缘存在较大差异,所以利用多成分信号的Wigner-Ville分布及其模的频率边缘能够确定该多成分信号的Wigner-Ville分布中交叉项的位置。
在上述步骤S113中,可以通过依次低通滤波,依次得到多成分信号中的各个单分量的子带信号。
本实施例中,首先利用多成分信号的Wigner-Ville分布及其模的频率边缘确定多成分信号的Wigner-Ville分布中交叉项的位置,再基于交叉项的位置进行低通滤波,不仅能够完全消除交叉项的干扰,而且不会降低Wigner-Ville分布的时频聚集度。
图6是本发明一实施例中利用多成分信号的Wigner-Ville分布及其模的频率边缘确定交叉项位置的方法流程示意图。如图6所示,在上述步骤S112中,利用所述多成分信号的Wigner-Ville分布及其模的频率边缘确定所述多成分信号的Wigner-Ville分布中交叉项的位置的方法,可包括:
步骤S1121:确定所述多成分信号的Wigner-Ville分布的模的频率边缘的极大值点;
步骤S1122:在所述极大值点处,计算所述多成分信号的Wigner-Ville分布的频率边缘和所述多成分信号的Wigner-Ville分布的模的频率边缘的差值;
步骤S1123:判断所述差值的模是否大于设定阈值,若是,将所述极大值点作为所述多成分信号的Wigner-Ville分布中交叉项的位置。
在上述步骤S1123,该阈值可以根据信号情况视需要设定。在多成分信号的Wigner-Ville分布的频率边缘和该多成分信号的Wigner-Ville分布的模的频率边缘的差值的模大于设定阈值时,可以认为多成分信号的Wigner-Ville分布与其模的频率边缘差别较大,可以认为该极大值点为交叉项所在位置,具体地例如是交叉项的中心点位置(频率位置),从而找到交叉项所在位置。
在其他实施例中,可以通过观察多成分信号的Wigner-Ville分布的频率边缘谱和该多成分信号的Wigner-Ville分布的模的频率边缘谱的差异来确定该多成分信号的Wigner-Ville分布中的交叉项位置。实施例中,可以在多成分信号的Wigner-Ville分布的频率边缘谱中去除交叉项位置。
一些实施例中,在上述步骤S113中,基于所述多成分信号的Wigner-Ville分布中交叉项的位置,通过低通滤波将所述多成分信号分离为所述多个子带信号的方法,可包括:
步骤S1131:以所述多成分信号的Wigner-Ville分布中交叉项的位置为截止频率,对所述多成分信号进行低通滤波,得到第一子带信号,其中,所述多个子带信号包括所述第一子带信号。
图7是本发明一实施例中基于多成分信号的Wigner-Ville分布中交叉项的位置通过低通滤波将多成分信号分离为子带信号的方法流程示意图。如图7所示,在上述步骤S113中,基于所述多成分信号的Wigner-Ville分布中交叉项的位置,通过低通滤波将所述多成分信号分离为所述多个子带信号的方法,除了包含上述步骤S1131,还可包括:
步骤S1132:利用所述多成分信号减去所述第一子带信号,得到剩余信号;
步骤S1133:利用所述剩余信号的Wigner-Ville分布及其模的频率边缘确定所述剩余信号的Wigner-Ville分布中交叉项的位置;
步骤S1134:以所述剩余信号的Wigner-Ville分布中交叉项的位置为截止频率,对所述剩余信号进行低通滤波,得到第二子带信号,其中,所述多个子带信号还包括所述第二子带信号。
本实施例中,从多成分信号中减除已分离出的子带信号,再利用类似的方法重新确定剩余信号的Wigner-Ville分布中交叉项的位置,并以重新确定的交叉项的位置对剩余信号进行低通滤波,依次进行,可以依次分离出其余的子带信号。
图8是本发明另一实施例的多成分信号的时频谱获取方法的流程示意图。如图8所示,在上述步骤S110之前,即将多成分信号分离为多个子带信号之前,还可包括:
步骤S140:利用窗极值经验模态分解对原始信号进行分解得到内蕴模态函数IMF分量信号;
步骤S150:根据所述IMF分量信号得到所述多成分信号。
在上述步骤S140中,该窗极值经验模态分解可以依据中国发明专利申请(公开号:CN 106844293 A)中的窗极值经验模态分解(Window-Extreme Empirical ModeDecomposition,WE-EMD)的方式实施。在上述步骤S150中,可以直接将该IMF分量信号作为待处理的多成分信号,或者可以将光滑后的IMF分量信号作为待处理的多成分信号。
针对模态混叠问题,中国发明专利申请(公开号:CN 106844293 A)提出的经验模态分解中模态混叠问题的自适应解耦方法中,提出窗极值经验模态分解(WE-EMD)处理模态混叠问题。噪声辅助法和窗极值构成了WE-EMD两个核心技术,它们相互依赖,缺一不可。噪声创建的二进制背景基为窗长的自适应选择提供了依据和支撑;同时,引进窗极值代替标准的局部极值构造上下包络,保证信号相近的固有模态映射到同一个特征子空间中,能够有效避免相近的特征时间尺度分布在相邻的IMF中。利用WE-EMD对信号进行预处理分解,然后计算WVD可以有效避免模态混叠引起的交叉项。
本实施例中,利用窗极值经验模态分解对原始信号进行分解得到内蕴模态函数IMF分量信号,可以在消除IMF分量信号的WVD产生的交叉项干扰的同时,有效避免模态混叠引起的交叉项。
一些实施例中,上述步骤S150,即根据所述IMF分量信号得到所述多成分信号的方法,具体实施方式可以是:对所述IMF分量信号进行小波包滤波,得到所述多成分信号。本实施例中,利用小波包滤波对所述IMF分量信号进行光滑处理,可以尽可能去除窗极值经验模态分解时保留在IMF分量信号中的噪声。
图9是本发明一实施例中多成分信号的时频谱获取方法的流程示意图。如图9所示,本实施例的多成分信号的时频谱获取方法,依次采用窗极值经验模态分解、光滑化、子带分离、计算WVD的处理流程。结合窗极值经验模态分解(WE-EMD)、WVD和基于频率边缘的子带分离方法,采用窗极值Wigner-Ville分布(WE-WVD)表示信号的时频谱。首先对信号进行窗极值经验模态分解,得到IMF分量信号;然后逐个计算各IMF分量信号的WVD,并利用WVD和其模的频率边缘之差的性质对IMF进行子带分离;最后计算各子带分离信号的WVD并求和。
实施例中,信号x(t)的Wigner-Ville分布定义为:
其中,Wx(t,Ω)表示信号x(t)的Wigner-Ville分布;t表示时间变量;Ω表示频率变量;表示信号,表示的共轭;τ是积分变量。
由式(1)所示的Wigner-Ville分布的定义可知,Wigner-Ville分布为双线性函数。用两个单分量信号x1(t)和x2(t)表示信号x(t),可得到x(t)=x1(t)+x2(t),则两个信号x1(t)和x2(t)之和的Wigner-Ville分布可表示成:
Wx(t,Ω)=Wx1(t,Ω)+Wx2(t,Ω)+2Re(Wx1,x2(t,Ω)) (2)
式(2)中,Wx1(t,Ω)和Wx2(t,Ω)分别表示信号x1(t)和x2(t)的Wigner-Ville分布。两个信号之和的WVD并不等于它们各自WVD求和。式中2Re(Wx1,x2(t,Ω))是信号x1(t)和信号x2(t)的互WVD,称之为交叉项,它是由信号x1(t)和x2(t)相加所引进的干扰。当信号的成分增多时,交叉项干扰会造成分析结果出现频率混叠的现象,从而影响信号分析的精度。因此如何抑制交叉项干扰一直是WVD方法应用中需要解决的问题。
一个实施例中,结合窗极值经验模态分解、WVD和基于频率边缘的子带分离方法,提出一种新的信号时频谱获取方法,可称之为窗极值魏格纳-威尔分布(Window ExtremeWigner-Ville Distribution,WE-WVD)。对于单分量信号,式(2)不会出现交叉项干扰,WVD能得到聚集度高的时频分布。结合窗极值经验模态分解和基于频率边缘的子带分离方法,先将信号分解成单分量信号之和,然后计算各单分量信号WVD之和得到分辨率高的时频分布。窗极值经验模态分解主要用于解决模态混叠问题,避免IMF分量信号中包含多个不同时间特征尺度的成分信号。基于频率边缘的子带分离方法将第1个IMF分量信号和包含其它低能量特征时间尺度的IMF分解成单个成分的信号之和,进一步保证各分量信号的单成分特性,提高WVD的时频聚集度。
假设待分解的信号x(t)={xi,i=1,2,…N},其中N是信号的长度,xi表示第i个单分量信号。采用如图9所示的流程计算信号x(t)的窗极值魏格纳-威尔分布(WE-WVD),可包括步骤:
步骤1:利用窗极值经验模态分解(WE-EMD)对信号x(t)进行自适应分解得到IMF分量信号t为时间,n为正整数,i表示IMF分量信号的序号;
步骤2:利用小波包滤波得到光滑化后的IMF分量信号ci(t),i=1,2,…n;
步骤3:利用基于频率边缘的子带分离方法对所有的IMF分量信号进行再次分解,第1个IMF分量信号和包含低能量特征时间尺度的IMF分量信号将被分解成多个单成分信号之和,其它IMF信号将保持不变,记子带分离信号为di(t),i=1,2,…m,m为正整数,i表示分离出的子带信号的序号;
步骤4:计算子带分离信号di(t),i=1,2,…m的WVD,记为Wi(t,Ω),i=1,2,…m。
步骤5:计算各子带分离信号WVD之和得到信号x(t)的时频谱:
其中,m表示子带信号的个数,i表示分离出的子带信号的序号,t表示时间,Ω表示频率,Wi(t,Ω)表示第i个子带信号的WVD,Wx(t,Ω)表示信号x(t)的WVD。
实施例中,上述步骤S140中的窗极值经验模态分解的方法可以参照中国发明专利申请(公开号:CN 106844293 A)中的经验模态分解的方法实施,例如,利用窗极值经验模态分解(WE-EMD)对信号x(t)进行分解得到IMF分量信号,可包括步骤:
(1)在信号x(t)中添加噪声,记c为常数,σx是信号x(t)的方差,rand(t)是与x(t)长度相等的噪声,得到包含噪声的信号xn(t);
xn(t)=x(t)+c*σx*rand(t)
(2)确定信号xn(t)的所有局部极值点;
(3)从局部极值点中选出窗极值点;
(4)用三次样条曲线构造xn(t)的上(窗极大值点)包络线和下(窗极小值点)包络线
(5)计算上、下包络线的均值
(6)求信号与均值的差
(7)重复以上步骤,直到差满足停止条件,得到第1个含噪声的IMF分量信号它表示信号在局部时刻频率最高的成分;
(8)用信号xn(t)减去第1个含噪声的IMF分量信号继续重复以上过程,直到分离出全部IMF分量信号,得到:
其中,表示第i个含噪声的IMF分量信号,表示信号xn(t)中除IMF分量信号外的剩余信号。
取式(2)中交叉项Wx1,x2(t,Ω)模的平方可得(陈章位,路甬祥.Wigner分布中互谱项特征及其消除方法的探讨[J].数据采集与处理,1995(1):1-5):
其中,Wx1和Wx2分别表示信号x1(t)和x2(t)的自项的WVD,t表示时间,Ω表示频率,τ表示积分变量,ξ表示积分变量。
由上式可以看出,交叉项Wx1,x2(t,Ω)刚好分布在自项Wx1(t,Ω)和Wx2(t,Ω)的正中间。假定自项Wx1(t,Ω)和Wx2(t,Ω)分别存在于中心点和的邻域内,则交叉项存在于中心点的邻域内(孟小芬,杜文超,高学强,等.Wigner-Ville分布交叉项识别方法研究[J].海军航空工程学院学报,2006,21(1):187-191)。
如果x1(t)和x2(t)的频谱分别限定在[ω1,ω2]和[ω3,ω4]之间,而且ω2<ω3,由WVD的定义得到:
Wx1(t,Ω)可表示为
Wx2(t,Ω)可表示为
Wx1,x2(t,Ω)可表示为
对式(2)两边在时域内积分,由频率边缘性质可知:
其中,Wx(t,Ω)表示信号x(t)的WVD,X(Ω+θ/2)表示信号x(t)的功率谱,X*(Ω-θ/2)表示X(Ω+θ/2)的共轭,δ(θ)表示狄拉克函数,θ表示积分变量,X1(Ω)表示信号x1(t)的功率谱,X2(Ω)表示信号x2(t)的功率谱,Wx1(t,Ω)和Wx2(t,Ω)分别表示信号x1(t)和信号x2(t)的WVD,Re(Wx1,x2(t,Ω))表示信号x(t)的WVD中的交叉项,t表示时间。
由上式可推导得到:
∫2Re(Wx1,x2(t,Ω))dt=0 (8)
令ω1、ω2、ω3、ω4、ω5及ω6表示频率值,结合式(6)~式(8)可知:
由上式可知,在交叉项区域[ω5,ω6],WVD的频率边缘是0,它与WVD模的频率边缘的差刚好是交叉项模的频率边缘。类似可证明,在Wx1(t,Ω)和Wx2(t,Ω)中心点邻域内,WVD与其模的频率边缘相等。如此一来,通过比较WVD与其模的频率边缘就能够找到交叉项的中心点,由式(5)可知它刚好位于自项Wx1(t,Ω)和Wx2(t,Ω)的正中间。图1中泄漏电流信号的第1个IMF分量信号的WVD和其模的频率边缘分别如图10的第1个IMF分量信号的WVD的频率边缘曲线101和第1个IMF分量信号的WVD的模的频率边缘曲线102所示。可以看出,在自项的中心50Hz、150Hz、250Hz,两者的频率边缘值是相等的,但是在交叉项的中心100Hz、200Hz处,两者的频率边缘值相差很大。以此可以验证发明人的上述理论分析。
基于上述理论推导,一个实施例中,可以通过比较WVD和其模的频率边缘,找出交叉项的中心点,并以其作为上限频率,利用低通滤波将x(t)分解成两个单分量信号x1(t)和x2(t)之和。可以把本发明实施例的上述方法称之为基于频率边缘的子带分离方法,具体实施例中,可包括步骤:
(1)计算IMF分量信号ci(t)的WVD,记为Wi(t,Ω);
(2)计算Wi(t,Ω)和其模的频率边缘,分别记为Pi(Ω)和
(3)提取Wi(t,Ω)模的频率边缘的极大值点,记为ωi,max,i=1,2,…K,K表示极大值点的个数;
(4)给准备进行子带分离的初始信号赋初值hi(t)=ci(t),hi(t)表示需要进行子带分离的信号,初值为第i个IMF分量信号;
(5)对极大值点进行循环,计算极大值点处Pi(ωi,max)与差,若其模大于阈值,则以ωi,max为截断频率对h(t)低通滤波,得到子带信号di,j(t),更新hi(t)=hi(t)-di,j(t),di,j(t)表示从第i个IMF分量信号分离出的第i个子带信号。
子带分离完成后,得到IMF分量信号ci(t)可以表示为:
利用上述实施例的方法对图1中的第2个信号即第1个IMF分量信号进行分解,结果如图11所示。图11中第1个信号是待分解的第1个IMF分量信号,第2~第4个是分解后的子带信号。
对图11中的各子带分离信号分别计算WVD,然后合并得到第1个IMF分量信号的WVD,结果如图12(a)所示。从图11和图12(a)可以看出,基于频率边缘的子带分离方法能够完整地把第1个IMF分量信号分解成3个频率分别是50Hz、150Hz、250Hz的单成分分量信号,解决了第1个IMF分量信号频带过宽和低能量的高频信号叠加在一起引起的多成分分量问题,消除地直接计算WVD的交叉项干扰,大大提高了时频谱的聚集度。
利用本发明一实施例的窗极值魏格纳-威尔分布(WE-WVD)计算泄漏电流信号的时频谱,先计算各子带分离信号的WVD,然后求和得到它的时频谱,结果如图12(b)所示。利用其它方法,包括WVD、伪WVD、重排WVD、连续小波变换、基于EMD的WVD、基于EEMD的WVD,得到的时频谱如图13(a)至图13(f)所示。比较图12(b)和图13(a)至图13(f)可以看出,利用本发明实施例的WE-WVD方法得到的时频谱的分辨率是最高的,可以清晰地分辨出来泄漏电流信号包含频率为50Hz的基频信号和150Hz的3倍倍频信号,而且没有交叉项干扰,能有效诊断污闪病害。从WVD、伪WVD、重排WVD的时频谱中虽然也能看出基频和3倍倍频信号,但交叉项干扰严重,模糊了信号的物理意义。连续小波变换得到的时频谱没有交叉项的干扰,但时频聚集度明显比WE-WVD差,而且3倍倍频特征不明显,影响污闪病害的诊断。从图13(e)可以看出,直接利用EMD和WVD计算信号的时频谱,由于模态混叠、第1个IMF分量信号频带过宽的影响,得到的时频谱分辨率很低,而且物理意义模糊。结合EEMD和WVD虽然能部分解决模态混叠问题,但多成分的第1个IMF会导致时频谱中出现交叉项干扰。
一个实施例中,基于频率边缘的子带分离方法。根据IMF分量信号的WVD和其模的频率边缘的差别,来分辨交叉项的频率中心,并以其作为上限频率实行低通滤波,将多分量的IMF分量信号转换成多个单分量信号。一个实施例中,窗极值魏格纳-威尔分布(WindowExtreme Wigner-Ville Distribution,WE-WVD)。首先对信号进行窗极值经验模态分解,得到IMF分量信号;然后逐个计算各IMF分量信号的WVD,并利用WVD和其模的频率边缘之差的性质对IMF进行子带分离;最后计算各子带分离信号的WVD并求和得到信号的时频谱。
基于与图4所示的多成分信号的时频谱获取方法相同的发明构思,本申请实施例还提供了一种多成分信号的时频谱获取装置,如下面实施例所述。由于该多成分信号的时频谱获取装置解决问题的原理与多成分信号的时频谱获取方法相似,因此该多成分信号的时频谱获取装置的实施可以参见多成分信号的时频谱获取方法的实施,重复之处不再赘述。
图14是本发明一实施例的多成分信号的时频谱获取装置的结构示意图。如图14所示,本实施例的多成分信号的时频谱获取装置,可包括:子带分离单元210、WVD计算单元220及时频谱计算单元230,上述各单元顺序连接。
子带分离单元210,用于:将多成分信号分离为多个子带信号;
WVD计算单元220,用于:计算所述子带信号的Wigner-Ville分布;
时频谱计算单元230,用于:计算各所述子带信号的Wigner-Ville分布的和,得到所述多成分信号的时频谱。
本发明实施例还提供一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现上述各实施例所述方法的步骤。
本发明实施例还提供一种计算机设备,如图15所示,计算机设备300可包括存储器310、处理器320及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述程序时实现上述各实施例所述方法的步骤。
综上所述,本发明实施例的多成分信号的时频谱获取方法、装置、存储介质及计算机设备,通过先将包含多分量成分的多成分信号分离为仅包含但分量成分的子带信号,再对各子带信号的Wigner-Ville分布求和得到该多成分信号的时频谱,能够消除直接计算多成分信号的Wigner-Ville分布(时频分布)中产生的交叉项的干扰,从而能够解决谱分布难以分析的问题,并能够提高多成分信号时频分析的精度。
在本说明书的描述中,参考术语“一个实施例”、“一个具体实施例”、“一些实施例”、“例如”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中,对上述术语的示意性表述不一定指的是相同的实施例或示例。而且,描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。各实施例中涉及的步骤顺序用于示意性说明本发明的实施,其中的步骤顺序不作限定,可根据需要作适当调整。
本领域内的技术人员应明白,本发明的实施例可提供为方法、系统、或计算机程序产品。因此,本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。