循环平稳

循环平稳（共6篇）

循环平稳 篇1

0 引言

脉冲位置相位翻转键控 (Pulse Position Phase Reversal Keying, 3PRK) 是Walker先生提出的基于时刻相位翻转典型的超窄带 (Ultra Narrow Bandwidth, UNB) 调制技术, 其调制原理是, 将数字基带脉冲在固定时刻控制正弦载波翻转而形成的通信信号, 具有典型的循环平稳特性。

循环平稳信号是一种特殊的非平稳信号, 在分析和处理上既不同于平稳信号, 又有别于一般的非平稳信号。这主要是因为信号统计量变化的周期性是可以充分利用的重要信息[1]。文献[2]证明了3PRK已调信号是循环平稳信号并对其正交性进行了分析, 文献[3]对3PRK信号的功率谱进行了分析, 但对其循环谱密度等性质都没有论述。

本文通过对3PRK已调信号循环谱幅度的不同截面进行了系统的分析, 提出可以利用该类信号的f=±fc的a截面循环谱进行噪声中的信号检测和参数估计, 并对这一设想进行了仿真验证, 取得了较好效果。

1 循环谱密度函数

若x (t) 为循环平稳信号, 且具有循环遍历性, 以T为周期性的循环自相关函数Ras (τ) 定义为[2]:

$R{}_s^{\alpha }(\tau )=\underset{{\rm T}\to \infty }{{\displaystyle \lim }}{\displaystyle \int }\begin{array}{l}{\rm T}/2\\ -{\rm T}/2\end{array}\text{x}(\text{t}+\frac{\text{τ}}{2})\text{x}{}^{*}(\text{t}-\frac{\text{τ}}{2})e{}^{-j2\pi \text{α}\text{t}}d\text{t}$, (1)

式中, α为循环频率, R${}_{\text{s}}^{\text{α}}$ (τ) 为频率为α的循环自相关强度。一个循环平稳信号的循环频率α可能有多个, 其中α=0对应信号的平稳部分, 因此只有非零的循环频率才刻画信号的循环平稳性。

循环自相关函数R${}_{\text{s}}^{\text{α}}$ (τ) 的Fourier变换, 即循环谱密度 (cyclic spectrum density, CSD) 表示为:

$S{}_s^{\alpha }(f)={\displaystyle \int }\begin{array}{l}+\infty \\ -\infty \end{array}\text{R}{}_{\text{s}}^{\text{α}}(\text{τ})e{}^{-j2\pi \text{f}\text{τ}}d\text{τ}$。 (2)

工程上采用如下估算方法:

XT (t, f) =∫t+T/2-T/2x (u) e-j2πfudu, (3)

$\begin{array}{l}S{}_{X{}_{{\rm T}}}^{\alpha }(t,f)=\frac{1}{{\rm T}}X{}_{{\rm T}}(t,f+\frac{\alpha }{2})X{}_{{\rm T}}^{*}(t,f-\frac{\alpha }{2}),(4)\\ S{}_X^{\alpha }=\underset{{\rm T}\to \infty }{{\displaystyle \underset{\text{Δ}t\to \infty }{{\displaystyle \lim }}}}\frac{1}{\text{Δ}t}{\displaystyle \int }\begin{array}{l}\text{Δ}t/2\\ -\text{Δ}t/2\end{array}\text{R}{}_{\text{X}{}_{\text{Τ}}}^{\text{α}}(\text{t},\text{f})d\text{t}‚(5)\end{array}$

式中, Δt为接收数据的长度, T为短时傅里叶变换的窗长。

2 3PRK已调信号循环谱分析

设sn是输入比特流序列, 服从均匀分布。那么3PRK调制信号可以用如下式表示,

$x(t)={\displaystyle \sum _n(}s{}_n\cdot g{}_1(t-n{\rm T})+(1-s{}_n)\cdot g{}_0(t-n{\rm T}))$, (6)

式中, g0 (t) 、g1 (t) 分别是对应一个符号周期发送的“0”、“1”波形。

$\begin{array}{l}g{}_0(t)=A\sin 2\text{π}f{}_{c}t,0\le t\le {\rm T}‚\\ g{}_1(t)=\{\begin{array}{l}A\text{s}\text{i}\text{n}(2\text{π}f{}_{c}t+\text{π})‚\\ A\text{s}\text{i}\text{n}(2\text{π}f{}_{c}t)‚\end{array}\begin{array}{l}0\le t\le \tau \\ \tau \le t\le {\rm T}\end{array}?,\end{array}$

式中, fc为载波频率, A为载波的幅度。fc=N/T, τ=n0T/N。N表示每个符号周期中有N个载波周期数, 一般N≥10;n0表示相位翻转的载波个数, 一般取1～2。

3PRK已调信号的循环功率谱密度为:

$S{}_x^{\alpha }(f)=\frac{1}{{\rm T}}{\displaystyle \int }\begin{array}{l}+\infty \\ -\infty \end{array}{\displaystyle \int }\begin{array}{l}\text{Τ}/2\\ -\text{Τ}/2\end{array}\text{x}(\text{t}+\frac{\text{τ}}{2})\text{x}{}^{*}(\text{t}-\frac{\text{τ}}{2})e{}^{-j2\pi \text{α}\text{t}}e{}^{-j2\pi \text{f}\text{τ}}d\text{t}d\text{τ}$。 (7)

经推导, 上式可以化为,

$\begin{array}{l}S{}_x^{\alpha }(f)=\frac{1}{4{\rm T}}[G{}_1(f+\alpha /2)G{}_1^{*}(f-\alpha /2)-G{}_1(f+\alpha /2)\cdot \\ G{}_0^{*}(f-\alpha /2)]-G{}_0(f+\alpha /2)G{}_1^{*}(f-\alpha /2)+\\ G{}_0(f+\alpha /2)G{}_0^{*}(f-\alpha /2)]+\frac{1}{{\rm T}{}^2}{\displaystyle \sum _m\delta }(f+\\ \alpha /2-\frac{m}{{\rm T}})[G{}_1(f+\alpha /2)G{}_1^{*}(f-\alpha /2)+G{}_1(f+\alpha /2)\cdot \\ G{}_0^{*}(f-\alpha /2)+G{}_0(f+\alpha /2)G{}_1^{*}(f-\alpha /2)+\\ G{}_0(f+\alpha /2)G{}_0^{*}(f-\alpha /2)]‚(8)\end{array}$

其中, $G{}_0(f)=\frac{Af{}_c(\text{e}{}^{-\text{j}2\text{π}f\tau }-1)}{2\text{π}(f{}^2-f{}_c^2)}‚G{}_1(f)=\frac{Af{}_c(\text{e}{}^{-\text{j}2\text{π}f\tau }-2\text{e}{}^{-\text{j}2\text{π}f\tau }+1)}{2\text{π}(f{}^2-f{}_c^2)}$。

由式 (8) 可知, 3PRK信号具有S${}_x^{\alpha }$ (-f) =S-αx (f) 的性质。如图1 (a) 、图1 (b) , 当α=0时, $S{}_x^{\alpha }(f)=S{}_x(f)={\displaystyle \int }\begin{array}{l}+\infty \\ -\infty \end{array}\text{R}{}_{\text{x}}(\text{τ})e{}^{-j2\pi \text{f}\text{τ}}d\text{τ}$, 表示3PRK已调信号的平稳部分, 在f=mfb时取最大值, 其中fb为码元速率。当f=±fc时, α以fb为周期出现较大值。

3 3PRK信号循环平稳特性的噪声抑制性能

平稳噪声信号具有统计均值为常量, 自相关函数仅与时间间隔有关的特点, 因此其循环均值和循环自相关函数在非零循环频率上恒为0, 即平稳噪声不具有循环平稳性。因此可以利用循环平稳信号与平稳过程的这一不同, 将二者区分开来。考虑一个被平稳噪声污染的循环平稳信号, y (t) =x (t) +n (t) , 不失一般性, 假设信号x (t) 和噪声n (t) 循环独立, 则接收信号的循环自相函数为:

R${}_y^{\alpha }$ (τ) =R${}_s^{\alpha }$ (τ) +R${}_n^{\alpha }$ (τ) 。 (9)

因为α≠0时, R${}_n^{\alpha }$ (τ) =0, 所以, R${}_y^{\alpha }$ (τ) =Rαx (τ) 。这样就可以利用循环自相关函数对信号进行检测。

由式 (8) 可以看出, 当f=±fc时, |Sαx (f) |α=2fc (mN-1) 取得较大值, 对于平稳噪声, 包括高斯噪声和非高斯噪声, 由于其不具有循环平稳特性, 故其循环谱在α=0存在非零值, 而在α≠0其值为零, 这个特点提供了在噪声中对3PRK已调信号检测的有效途径。

4 仿真结果

为验证以上分析, 使用Matlab进行仿真分析。选取的3PRK信号的参数为N=10, n0=1, 幅度为单位1, 码元数组为100, 码元速率为200 bit/s/Hz, 载波fc=2 000 Hz。仿真采样频率fs=5fc, 窗长T=500Ts。经分析得到结果如图1、图2、图3所示。

图1 (a) 、图2 (a) 、图3 (a) 分别为不同信噪条件下3PRK信号的循环谱, 图1 (b) 、图2 (b) 、图3 (b) 为f=fc时循环谱α的截面谱。

图1 (a) 显示的是没有噪声情况下, 3PRK信号循环谱, 图2 (a) 、图3 (a) 分别显示的是信噪比为1 dB、5 dB时, 3PRK信号的循环谱。图1 (b) 、图2 (b) 、

图3 (b) 为相应信噪条件下, f=fc时循环谱α的截面谱。

从图1 (b) 、图2 (b) 、图3 (b) 可以明显观察到, 截面谱以2/T的频率存在固定的周期成分。另外从图2 (b) 、图3 (b) 可以观察到, 当信噪比为1 dB、5 dB时噪声在α≠0的截面也存在非零值, 这是因为在实际操作中取窗等的原因造成的, 但观察发现, 即使信噪比为1 dB, 经过处理, 也能估计出噪声中的信号的参数。因此, 在检测3PRK信号时可以利用f=±fc的α截面循环谱进行噪声中的信号检测和参数估计。

5 结束语

本文在简要介绍二阶循环理论基础上对3PRK信号的循环谱相关函数进行了推导, 并对3PRK已调信号循环谱在无噪声环境、信噪比为1 dB、5 dB的高斯白噪声环境进行了仿真分析。结果表明, 由于取窗等原因, 噪声在a≠0的截面也存在非零值, 但是在f=±fc处的α截面能够较好地估计出3PRK信号的参数, 因此, 3PRK信号具有很好抑制噪声的能力。循环平稳特性在3PRK信号上的运用, 为进一步实现噪声中3PRK信号的优化处理提供了一种方法。

参考文献

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

[2]张士凯.EBPSK调制波的正交性及功率谱分析[J].应用科学学报, 2008, 26 (2) :127-131.

[3]戚晨皓.采用复合调制的AM广播功率谱分析[J].应用科学学报, 2007, 25 (6) :583-588.

[4]陈四根, 杨莘元.QPSK信号谱相关性质研究[J].哈尔滨工程大学学报, 2003.24 (2) :208-211.

循环平稳 篇2

脉搏信号的每个脉动周期不完全一样,不是严格的周期信号,属于非平稳信号,但却呈现出一定周期平稳性,即脉搏信号具有循环平稳特性。目前,有关脉搏信号质量评估与滤波方法并没有考虑信号的循环平稳特性[3,4]。近年来,随着循环平稳理论的完善和发展,在通信、雷达等领域得到广泛的应用[5]。相比之下,对人体生理信号的循环平稳特性研究刚刚起步。与脉搏信号类似,心音信号也属于循环平稳信号,有学者将循环平稳算法应用于心音信号的质量评估和包络估计,取得了一定的研究成果[6],可作为脉搏信号质量评估的参考。于是,本文提出基于脉搏信号循环平稳特性的质量评估与滤波方法,定义了可反映脉搏信号质量变化的质量系数,用于评估脉搏信号受噪声污染的程度;对于受噪声污染小的信号段利用循环平稳算法与最大信噪比准则设计循环相关匹配滤波器进行降噪,并将污染严重的脉搏信号剔除。采用建模产生的仿真脉搏信号与本课题组实际采集的脉搏信号验证所提出方法的准确性和实用性。

1 方法

1.1 脉搏信号的循环谱

脉搏信号为循环平稳信号(记为x(t)),其周期在一个常数上下波动(该常数为理想情况下的脉搏周期,记为T )。根据循环平稳理论[7],任意时刻t,x(t)的循环自相关函数定义为:

其中,α为循环频率,τ为时间延迟,上标*表示取共轭。如果循环频率α=0,那么循环自相关函数Rxa(τ)就变成普通的自相关函数Rx(τ)。对Rxa(τ)做傅里叶变换可得:

Sxa(f)被称为循环谱密度函数。f是谱频率,如果循环频率α=0,那么Sxa(f)就变成了普通的功率谱密度函数Sx(f)。由于噪声会影响脉搏信号的循环平稳特性,所以不同噪声环境下信号的循环谱不同,对应的循环频率也不同。

1.2 脉搏信号质量系数

对于受到噪声污染的脉搏信号,需要判断其质量是否符合后期分析的要求。本文根据脉搏信号的循环谱提出质量系数,反映脉搏信号质量变化。由于循环平稳算法基于循环频率α而非谱频率f,因此,Sxa(f)对循环谱积分后可消除谱频率f。对式(2)积分得到:

其中 ,γx( α ) 叫做循环 频率频谱 密度 ( Cycle Frequency Special Density,CFSD)。据此再定义质量系数为:

式中,μ是最大循环频率;λ是基本循环频率,即γx(α)波形的第一个峰值所对应的频率位置[8]。由于噪声会影响脉搏信号的循环平稳特性,因此,通过质量系数的大小可以有效地反映脉搏信号受噪声污染的程度。

1.2.1 滑窗的选取

将式(4)用于脉搏信号,可得到整段脉搏信号的质量系数,为一常数,不能反映脉搏信号局部质量变化。于是,引入滑窗的思想,给脉搏信号加入一个动态的窗口(窗宽记2ζ),采用式(4)计算该窗口([t-ζ, t+ζ])内信号的质量系数作为t时刻的质量系数。则式(4)变为:

其中,γx(λ, t)是表示在[t-ζ, t+ζ]窗口内脉搏信号的循环频率频谱密度。从时变质量系数d(λ, t)可以看出,滑动窗的长度太短,无法体现脉搏信号的循环平稳特性,也就无法得到循环频率α。所以,滑动窗的宽度不能小于2个完整的脉搏周期。

1.2.2 循环谱估计

对于实际脉搏信号,循环频率α是未知的,所以循环谱Sxa( f )也是未知的,质量系数d(λ , t)也无法得到。于是,需要对脉搏信号的循环谱Sxa( f )进行估计。本文采用常用的非参数化时域平滑方法——分段谱相关算法(Strip Spectral Correlation Algorithm,SSCA) [11]得到估计值

根据互谱分析理论[9]与时域平滑理论[10],得到分段谱相关算法表达式:

其中,Ts为采样周期,fs=1/Ts为采样频率;△t为时间长度;N为在△t内的数据采样点数;XT(n, f )是输入信号x(t)的复解调,即x(t)的N'点快速傅立叶变换,N'一般等于或者大于2的整数次幂;q是循环频率分辨率△α倍数;g(n-r)是P×L阶的矩形平滑窗;循环频率分辨率△α=1/△t=fs/N;频率分辨率△f=fs/N'是一个不随循环频率α的改变而改变的常数,△f=△α;△t·△f=N/N'为一个常量。

利用SSCA算法得到循环谱的估计值之后,可通过其峰值对应的横坐标得到循环频率α,再运用式(3)和式(4)就得到质量系数d(λ, t),对脉搏信号进行质量评估。

1.3 循环相关匹配滤波器

根据质量系数反映出脉搏信号受噪声污染的程度,剔除质量差的信号段,对剩余信号进行滤波处理。结合脉搏信号的循环平稳特性及滤波器输出信噪比最大的要求,在匹配滤波器的基础上引入循环平稳算法,设计循环相关匹配滤波器(Cyclic Correlation Matched Filter, CCMF)对脉搏信号进行降噪。

设输入信号x(t)=s(t)+n(t),其中,s(t)是已知循环谱Sxa(f)和循环频率α的循环平稳信号;n(t)是噪声信号,与s(t)统计独立。根据循环相关匹配滤波器设计算法[11]得到滤波器传递函数:

其中Ssa( f )是信号s(t)的循环谱密度。

对应的滤波器时域传递函数为:

所以,循环相关匹配滤波器输出为:

其中,c是常数;Rs (τ0 - τ)是脉搏信号s(t)滞后时间常数为τ0的循环自相关函数。

2 数据

由于无法获得实际脉搏信号的真值,所以,无法对算法的准确性进行评估。于是,本文给干净的脉搏信号(真值信号)叠加噪声信号产生实验数据。干净脉搏波由三个高斯函数来合成[12,13],分别对应于脉搏波的主波、重搏波及重搏前波。每个高斯函数由3个参数确定,即幅度V、时间T和宽度U。合成的脉搏波p(t)由下式表示:

其中,V1= 0.8 , V2= 0.5 , V3= 0.4 , T1= 0.2 5 ,T2=0.45,T3=0.7,U1=0.012,U2=0.01,U3=0.03。

对式(11)的单个脉搏波进行延拓,产生一组脉搏信号,将其作为干净脉搏信号,1 s内信号含有250个数据点,信号总长度为10 s。为了消除信号幅值不同的影响,将其归一化(均值为0,标准差为1),如图1所示。噪声信号取自PysioNet BIH Noise Stress Test数据库[14]。该数据库包含了健康受试者在脉搏采集的3.5 h,采样频率为250 Hz的噪声信号,主要有基线漂移(bwm)和肌电干扰(mam)。截取与脉搏信号相同长度的噪声信号,加入干净的脉搏信号,如图2所示。

3 结果

3.1 质量系数结果

利用分段谱相关算法计算出被噪声污染脉搏信号循环谱估计值),分别通过式(3)与(5)得到循环频率频谱密度γx(α)与质量系数d(λ , t),如图3(a)与图3(b)所示。从图3(a)中看出,被噪声污染的脉搏信号循环频率范围是0~2.4 Hz,峰值点对应基本循环频率η=1.2 Hz。从图3(b)中看出质量系数在0~4 s段明显下降,在4~10 s段上升。其中质量系数d(λ, t)最小值为-1.7。而对于图2(a)中干净的脉搏信号,其基本循环频率η=1.25 Hz,质量系数d(λ, t)=-1.2。实验结果表明,质量系数随噪声的增大而减小,即质量系数可以准确地反映脉搏信号质量的变化。

3.2 滤波结果

3.2.1 循环相关匹配滤波器

将被噪声淹没的脉搏信号(如图2(b)中0~4 s段)剔除后,将剩余部分(如图2(b)中4~10s段)通过循环相关匹配滤波器滤波,结果如图4(a)所示。由图4可见,脉搏信号中的噪声被有效去除。

3.2.2 三种滤波方法的比较

采用整系数滤波法 ( 实时性高 , 但滤波效果一般)、EMD分解滤波法(滤波效果好,但实时性差)和本文方法进行对比。将图2(b)中4~10 s段脉搏信号分别通过设计的整系数滤波器与EMD滤波器[15,16]后得到的输出信号如图4(b)和4(c)所示。3种滤波器输出信号与原始干净的脉搏信号(图3(a))比较后可以看到,整系数滤波存在明显的群时延并且滤波效果不好;EMD滤波没有群时延但噪声未完全滤除;相比之下CCMF滤波效果明显好于整系数滤波和EMD分解滤波。

分别采用均方根误差(MSE)、信噪比(SNR)、运算时间(TIME)对三种滤波器的性能进行评价。计算三种滤波方法的MSE,SNR及运算时间(软件:Matlab R2010a,计算机配置:Intel E4600双核处理器,主频2.4 GHz,内存1.5 GB),对比三种方法的滤波性能。结果如表1所示。

注:MSE:均方根误差;SNR:信噪比(d B);TIME:运算时间(s)

通过表1可以看到CCMF与其它两种滤波方法相比,MSE最小(0.292 2),SNR最大(7.968 1),但是由于CCMF的输出需要计算循环频率,所以运算时间比整系数滤波器时间长。

综上所述,整系数的运算量小但是准确性不高,同时存在群时延;EMD准确性较好但是算法复杂度也高,不利于脉搏信号的实时分析处理;相比之下,CCMF不仅在准确性和时效性上较好,而且信噪比高,可以很好地保留脉搏信号的特征信息并有效地去除噪声。

3.3 应用

选取本课题组实际采集的一组脉搏信号,如图5(a)所示,前后部分存在受噪声严重污染的信号段,特征已完全丢失,中间部分(4~8 s)受噪声污染小。利用质量系数对该段脉搏信号的质量进行评估,结果如图5(b)所示,由于中间段受信号噪声污染小,所以其质量系数比两端大。通过质量系数可以反映脉搏信号质量的变化,将信号质量不好的,即信号被噪声淹没的部分(0~4s与8~10 s)剔除;对剩余部分(4~8 s,图6(a)所示)采用CCMF滤波,结果如图6(b)所示,可以看到CCMF能够有效地去除脉搏信号中的噪声。

4 总结

循环平稳 篇3

EMD[1]去噪分解方法应用广泛, 在地质勘探[2]、电力检测[3]、音视频处理技术[4]等领域都取得了较好效果。如一些大型发电厂、化工企业的大型齿轮箱有上百个齿轮, 发生故障时, 技术人员很难直接判断出故障部位, 如果逐个部位拆卸检查, 不仅影响工作进度, 而且会对其它好的齿轮产生影响。在EMD分解与滤波去噪的基础上, 通过分析大型齿轮箱振动信号并解调处理[5], 可以准确判断出故障部位。大型故障齿轮箱振动信号表现为故障调幅信号、调频信号、调频调幅信号等[5]。EMD分解算法虽然在一定程度上去除了噪声并分离有用信号, 但过度分解会对有用信号产生破坏。

1 单个调频、调幅、调频调幅信号EMD分解局限性

1.1 单个调幅信号EMD分解局限性

理论上, EMD去噪分解过程中, 单个调幅信号可以作为一个IMF分量[1]被分解出, 但SD值[1]越小, 信号被迭代次数越多, 信号的包络越趋于对称化, 因此在某一SD下能作为IMF分出的调幅信号在另一更小的SD下将被分解。从式 (1) 看出, 调幅信号可看作3个单频成分的叠加。因此, SD值很小, 调幅信号迭代次数过多, 原调幅信号被EMD分解完全破坏掉, 变成3个单频分量析出。如图1所示, 当SD=0.2时, 信号被完整析出, 当SD=0.002时, 信号被破坏分解成两部分 (7Hz、8Hz作为一部分, 6Hz作为另一部分) 。

1.2 单个调频信号EMD分解局限性

单个调频信号与单频信号的包络特性与信号局部均值特性完全一致, 因此SD值对信号分解效果影响不大, 即使SD值设置很小, 迭代次数增加, 但其信号仍然能被完整解析出。

对比图2与图3, 仿真信号x (t) =cos[2π10t+5sin (2π1t) ]采样频率60Hz, 采样点数1024点。可以看出, 取SD=0.2与SD=0.002对信号的分解无太大影响。

1.3 单个调频调幅信号EMD分解局限性

调频调幅信号的包络特性类似于调幅信号, 又类似于调频信号, 如式 (2) 所示。一个同时调频调幅的信号可看作3个频率成分交叉的纯调频信号叠加。

从式 (2) 可以看出, 上述3个调频信号的频率成分相互干扰, EMD分解先将高频成分析出, 再析出低频成分。从式 (2) 可以看出, 上述3个调频信号的频率成分相互干扰, EMD分解先将高频成分析出, 再析出低频成分。因此, 无论SD值小, 增加迭代次数, 都不能将调频调幅信号分解, 但过度分解会使信号包络均匀化, 破坏原信号, 出现虚假频率成分。仿真信号:

x (t) =[1+cos (2π1t) ]cos[2π10t+2sin (2π1t) ]对比图4中 (a) 、 (b) 、 (c) 小图, SD值越小, 过度迭代越多次, 信号变形破坏变形越严重, 信号产生虚假频率成分越多且复杂。

2 两个调频调幅叠加信号EMD分解局限性

调频调幅信号集合了调频信号与调幅信号的特征, 本文主要分析两个调频调幅信号叠加的EMD分解局限性, 如式 (3) :

通过仿真得出以下结论:

(1) 两调频调幅信号的中心频率离得越远 (即f2远大于f1) , 两个信号越容易分解开。因此, 在工程实践中, 当两个或多个信号中心频率接近时, 通过EMD分解法很难将它们分开, EMD分解对这种中心频率很接近的实际叠加信号的分解往往无能为力。

(2) 只有当高频信号能量比低频信号能量高很多时 (即A2远大于A1, B2远大于B1) , 才能将信号有效分开, 因此对实际工程信号要求很高, EMD方法对低频分量能量高于高频分量能量的实际工程信号分解效果不明显。

仿真信号1:

仿真信号2:

如式 (4) 、图5所示, 取低频成分能量为高频成分的2倍, 并不能将两个叠加信号分开, 还因过度迭代产生了大量的虚假频率成分。如式 (5) 、图6所示, 取高频成分能量为低频成分的2倍, 能将两个调频调幅信号分开。

采样点数1 024点, 采样频率60Hz SD值0.02, 低频部分极值点数374点, 高频部分极值点数511点, 叠加信号迭代后为375点。

(3) 调制系数值越小 (即C2、C1都很小) , 则信号调制边频带越窄, 两个调频调幅成分越不容易发生干涉, 分解效果越好。EMD分解对各频率分量边频带较宽的实际工程信号的分解效果较差。

综上所述, 当高频调频调幅分量和低频调频调幅分量由于中心频率接近或边频带较宽产生频率干扰时, 或者叠加信号中低频成分能量高于高频成分能量时, 叠加信号并不能将其分开。此时, 过度分解将出现假频。当两调频调幅叠加信号不能分解时, SD值越小, 迭代次数越多, 则过度迭代产生的假频将越多越复杂。

3 结语

(1) 对于单个信号, 单个调幅信号可以看作是3个单频信号的合成, 如被EMD过度分解, 将解体为3个单频信号;单个调频信号最稳定, 不易被EMD过度迭代影响;单个调频调幅信号可以看作是3个频带重叠的单独调频成分的叠加, 过度EMD分解不能使其分解, 但会产生复杂的虚假频率成分。

(3) 对于叠加信号, 无论哪种信号, 两个叠加信号被分开的条件是两个信号频率成分隔开, 即高、低频率分量中心频率离得越远越好, 高、低频率分量各自边频带越窄越好, 边频带不发生干扰, 而且要求高频信号分量的能量大于低频信号分量。但实际工程信号很难满足上述条件, 因此往往造成EMD去噪分解效果不佳。

(3) 当叠加信号不能被分解时, 信号被过度迭代, 严重变形, 产生假频成分, SD设置值越小, 假频现象越严重。

(4) EMD分解方法适合于频率成分少, 各IMF分量间频率成分差异较大, 不产生频率干涉, 各IMF分量局部均值差异较大, 能量从高频分量到低频分量逐步减小, 各分量包络不易干扰的非平稳信号。故障齿轮箱振动调制信号是循环平稳信号, 频率成分非常复杂, 各分量局部均值都可看作是零, 各分量间频率干涉与包络干涉非常严重, 其EMD分解去噪算法对这类相互重叠干涉的复合循环平稳信号效果较差。

摘要：利用matlab软件进行仿真分析, 探讨各类循环平稳信号EMD去噪分解算法的局限性。总结了单独调幅信号、单独调频信号、单独调频调幅信号EMD去噪分解的局限性。

关键词：EMD分解,调幅信号,调频信号,调频调幅信号

参考文献

[1]HUANG NE, SHEN Z, LONG SR, et al.The pirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[J].The Rayal Society, 1998 (454) :903-995.

[2]陈凯.基于经验模式分解的去噪方法[J].石油地球物理勘探, 2009 (5) :603-608.

[3]孙金宝, 朱永利, 刘丽轻, 黄建才.基于EMD的绝缘子泄漏电流去除噪声研究[J].华北电力大学学报:自然科学版, 2010 (6) :1-5.

[4]李烨.基于EMD的图像拼接和图像识别研究与实现[J].重庆:重庆交通大学, 2010.

循环平稳 篇4

正交频分复用(OFDM)宽带通信技术作为一种高频带利用率的多载波调制技术,具有很强的抗多径、抗衰落能力,将成为下一代无线通信系统物理层关键技术之一,但由于多载波调制信号的特征不易提取,因此越来越多的人关注于OFDM信号的盲识别算法研究。传统的OFDM信号盲识别算法大致可以分为两类:从时域上考虑,OFDM信号表现出渐近高斯性,而单载波信号是非高斯的,而高斯随机变量的高阶(>2)累积量恒等于零,所以利用高阶累积量对OFDM信号进行识别,而这些算法有一个严格的限制条件,就是当子载波数目很大时,接收信号才具有渐近高斯性;从频域上考虑,单载波信号在频域上,一般只是很窄的一个脉冲,而多载波信号则会呈现出比较多的窄脉冲,通过统计脉冲的个数,就可以判断出信号所使用的调制是否为多载波,由于基于频域幅度分析的方法受观察者的感官判断影响,不易得到准确的结果。

针对以上识别算法的不足,本文将利用信号循环平稳性知识对该问题进行分析,因为数字通信信号一个突出特点表现为循环平稳性[1],将信号建模为循环平稳信号更为合适,本文从理论上分析了可以区分单载波调制信号和OFDM信号的方法,最后通过仿真对分析结论进行验证。

2 循环累积量知识

循环平稳信号x(t)的n阶时变矩函数(TMF)定义为n阶滞后积的期望值:

其中τ是一个滞后向量:τ=[τ1,τ2,…,τn],τj表示延迟;q(q=0,1,…,n/2)表示n阶滞后向量积中取共轭因子总数,(*)j是根据共轭因子总数q选择x(t)的第j个因子x(t+τj)是否取共轭;mxα(τ)n,q表示信号x(t)的n阶循环矩,其表示为:

与n阶时变矩函数(TMF)相对应的n阶时变累积量函数(TCF)可以由M-C (Moment to Cumulants)转换得到

集合I={1,2,…,n},Ip表示集合I的无交连的非空分割,分割为Q。x(t)在循环频率β处的n阶循环累积量可表示为:

其傅里叶变换表示为:

称为n阶循环多谱。

循环累积量最重要的性质就是对感兴趣信号的选择性,在多种调制类型信号环境下、不同符号周期的信号模型可表示为:

不失一般性,假定感兴趣的信号为r0(t,T0),并且T0≠nTj,j≠0,n是正整数。由于k个信号只有r0(t,T0)在循环频率1/T0处有非零的循环累积量值,可以得到利用循环累积量的这个性质,可以将感兴趣的信号类型从多种调制类型信号中检测识别出来。

3 单载波调制信号模型及其循环累积量

经过AWGN信道的单载波接收信号表示为

式中,sk表示第k个符号周期内的发送符号,它是零均值独立同分布的复符号序列;a,T,Δfc,θ,t0分别表示符号幅度、周期、载波频率偏差、相位旋转和时延:w(t)是宽平稳复加性高斯白噪声,p(t)表示脉冲成型函数。

单载波信号rSCLD(t)在循环频率集γ处的n阶循环累积量为:

式中,Cs,n,q表示序列符号Sk的n阶累积量(q个共轭);q=0,1,...,n/2表示n阶滞后积中取共轭的因子总数;(*)m是根据共轭因子总数q选择p(t)的第m个因子p(t+τm)是否取共轭;(一)m是为负值和(*)m是否取共轭一致;τ=[τ1,…,τn-1,0]表示延迟向量;循环频率集γ=β+(n-2q)Δfc;循环频率β=k/T,k为整数。

4 OFDM信号模型及其循环累积量

4.1 OFDM信号模型

经过AWGN信道的OFDM接收信号表示为

式中,sk,l表示第l个符号周期内第k个载波的发送符号,它是零均值独立同分布的复符号序列;a,θ,t0分别表示符号幅度、相位旋转和时延,w(t)是宽平稳复加性高斯白噪声;T=Tu+Tcp表示符号周期;Tu表示有用符号时长;Tcp表示循环前缀;ΔfK=I/Tu表示子载波间隔;p(t)表示脉冲成形函数。

4.2 OFDM信号循环累积量推导[4]

根据累积量的共线性特征[5],可以得到接收到的OFDM信号的n阶(q个共轭)时变累积量:

式中⊗表示卷积运算;对于高斯白噪声而言,当n≥3时,cw(t,τ)n,q≡0,在下面的推导过程中,将不考虑高斯白噪声部分。对式(10)进行傅里叶变换得:

利用等式:

代入式(11)得:

从式中可以看出,只有当时,,对式(13)进行傅里叶反变换可表示为:

根据循环累积量理论[6],时变累积量可以表示为傅里叶级数的形式:

式中,kn,q={γ:cγ(τ)n,q≠0}为循环频率集,c(γ,τ)n,q表示为n阶(q个共轭)循环累积量。

可以得到OFDM信号n阶循环累积里为:

5 循环累积量特征分析及仿真验证

现在选取接收信号二阶(n=2,q=1)循环累积量来研究信号特征量之间的区别,这样选取的好处是:当q=n/2时可消除载波频率偏差对循环频率的影响。则接收的单载波信号和OFDM信号二阶循环累积量分别表示为:

仿真参数设置:

单载波信号:调制方式为MQAM或MPSK,符号长度T=1μs,升余弦成型函数。

OFDM信号:子载波调制方式为MQAM,子载波数N=16,符号长度T=1μs,循环前缀长度Tcp=0.2μs,矩形脉冲成型函数。

比较式(17)、(18)可以看出,OFDM信号二阶循环累积量表达式多一个乘积因子:

由上式可知,当τ=l/ΔfK,l为整数时,OFDM信号二阶循环累积量幅值会有显著峰值,而单载波信号则不存在这一峰值,如图1、图2所示:分别表示单载波信号与OFDM信号二阶循环累积量幅度图,其中,时延τ相对T进行了归一化,循环频率α相对T-1进行了归一化。

进一步分析OFDM信号的二阶循环累积量

,p(t)为矩形脉冲,定义为:,代入式(18)可得:

由于OFDM信号传输中循环前缀的引入,结合上式分析得,当τ/△fK时,T≠|τ|,因此OFDM信号峰值离散,当β=n/T时,有最大值,其他条件下为零(如图3所示),我们正是利用这些离散的峰值来识别单载波信号和OFDM信号。

6 结语

数字通信信号为一类特殊的非平稳信号,本身更适宜建模成循环平稳信号。本文从理论上分析了单载波信号和OFDM信号二阶循环累积量区别:由于OFDM信号传输中循环前缀的引入,当τ=l/ΔfK,T≠|τ|,β=n/T时,OFDM信号循环累积量幅值有离散峰值,在并通过计算机仿真,仿真结果与理论结果一致,为下一步工作,单载波信号与OFDM信号检测和识别做了理论基础工作。

摘要：文章介绍了信号循环累积量知识,推导了OFDM信号循环累积量的表达式,并以二阶循环累积量为例,从理论上分析了单载波信号和OFDM信号的区别,最后通过仿真验证了分析结论。

关键词：循环累积量,单载波信号,OFDM信号

参考文献

[1]黄奇珊,彭启琮,路友荣等.OFDM循环谱结构分析[J].电子与信息学报,2008,30(1):134-138

[2] Octavia A.Dobre and Yeheskel Bar-Ness.Higher-order Cyclic cumulants for high order modulation classification,In Proc.IEEE MILCOM.2003:112-117

[3] M.Oner and F.Jondral.Air interface recognition for a software radio system exploiting cyclostationarity.ln Proc. IEEE PIMRC.2004

循环平稳 篇5

由于人为调制等原因,大多数雷达和通信信号都体现出一种特殊的周期性,即循环平稳特性,不同调制样式的信号具有不同的循环平稳特性[1],因此,基于信号循环平稳特性的参数估计方法在信号的选择性和噪声抑制方面具有显著优势[2]。

阵列信号处理是信号循环平稳特性成功应用的一个典型例子[2]。GARDNER W A首先把信号的循环平稳特性引入到阵列测向方法中,提出了相应的Cyclic MU-SIC(简记为C-MUSIC)方法[3]。文献[4,5]分析了C-MUSIC方法对高斯和非高斯窄带循环平稳信号测向的性能,但该分析都是基于理想阵列结构进行的。随后,CHARGE P等通过综合考虑入射信号的循环相关函数和共轭循环相关函数,对已有的C-MUSIC方法进行了扩展,提出了Extended Cyclic MUSIC(简记为EC-MUSIC)方法[6]。与C-MUSIC方法相比,EC-MUSIC方法更有效地利用了入射信号的信息,因此具有更优的测向性能。当入射信号带宽较为显著时,YAN H对C-MUSIC方法进行了修正,提出了Improved Cyclic MUSIC方法[7],极大地减小了实际信号的非零带宽在C-MUSIC方法中所带来的测向偏差。该修正思想同样可用于EC-MUSIC方法中,借以减小对非零带宽近似窄带信号的测向偏差[7]。

在实际系统中,接收阵列可能存在各种模型误差,如阵列互耦、通道不一致性等。当此类误差存在时,C-MUSIC方法和EC-MUSIC方法可能产生测向偏差。本文以常规阵列中很难消除的阵列互耦效应为对象,通过建立互耦条件下Cyclic MUSIC方法和EC-MUSIC方法的伪数据模型,考查这两种方法的测向性能受互耦效应的影响情况,给出了均匀线阵和一般线阵中测向偏差的解析表达式,并通过仿真实验验证了信号入射方向和阵元间距变化时该理论偏差的准确性。

1 互耦条件下窄带循环平稳信号测向方法的伪数据模型

假设K个循环频率为α的窄带信号从不同方向Θ=[θ1,θ2,…,θK]同时入射到M元均匀线阵上,则阵列接收数据模型为:

其中x(t)=[x1(t),x2(t),…,xM(t)]T为阵列在t时刻的接收数据构成的向量,s(t)=[s1(t),s2(t),…,sK(t)]T为各入射信号在t时刻的波形,(△k为第k个信号在相邻阵元间的传播时延)为K个信号的导向矢量,为表述方便,以下的讨论中简记A(f,Θ)为A(f),并特别标记A0=A(f0,Θ),n(t)=[n1(t),n2(t),…,nM(t)]T为阵列加性噪声。另外,记信号sk(t)的幅度包络为gk(t),即。

互耦条件下阵列观测数据为[8]:

其中分别为互耦条件下和理想情况下t时刻阵列的接收数据,C为阵列互耦矩阵。在以下的叙述中,如无特别说明,均表示变量Φ在互耦条件下的取值。

第p个阵元t时刻的接收数据为:

其中cpu为C的第(p,u)个元素。

1.1 互耦条件下C-MUSIC方法的伪数据模型

互耦条件下第p个阵元和第q个阵元接收数据的循环互相关函数为:

其中为共轭运算符,为共轭转置运算符,且:

进一步地,。

结合文献[3]给出的无互耦伪数据矩阵,得到互耦条件下考虑窄带信号非零带宽时的阵列伪数据矩阵为:

1.2 互耦条件下EC-MUSIC方法的伪数据模型

经过类似的推导,得到互耦条件下第p个阵元和第q个阵元接收数据的共轭循环互相关函数为:

因此,为转置运算符。

进而得到互耦条件下EC-MUSIC方法的阵列伪数据矩阵为:

其中,。

结合文献[6]给出的理想情况下EC-MUSIC方法的伪数据矩阵,得到互耦条件下考虑窄带信号非零带宽时的阵列伪数据矩阵为:

2 互耦条件下C-MUSIC方法的测向误差分析

对YC(α,τ)进行特征值分解,得到互耦条件下C-MUSIC方法的信号子空间和噪声子空间:

结合式(6)、式(10)不难看出,信号子空间与真实信号方向矩阵A(f0+α/2)之间存在如下关系:

其中分别表示由矩阵的各列所张成的子空间。

如果互耦条件下,空间多个目标仍然可以分辨,则由空间谱函数:

估计得到的信号波达方向所构成的阵列导向矢量与真实导向矢量间满足如下关系:

其中ρ为幅度调整系数。

记,为相应的互耦条件下对应的变量。由互耦条件下阵列接收数据的循环互相关协方差矩阵的特征值分解得到的估计值之后,可用如下罚函数确定λk:

参考文献[9]中的数学推导,可以由以上罚函数得到互耦条件下对第k个信号角度θk的估计误差为:

其中,c0为电磁波传播速度,1为全1列向量,其维数由上下文确定,分别为复数的实部和虚部,⊙表示矩阵或向量点乘运算符。

3 互耦条件下EC-MUSIC方法的测向误差分析

对YEC(α,τ)进行特征值分解,得到互耦条件下EC-MUSIC方法中伪数据协方差矩阵的信号子空间和噪声子空间:

互耦条件下波达方向的估计结果由如下空间谱估计函数得到[6]:

为了消除加权向量h的影响,将上述空间谱函数转化为:

由于为二阶单位矩阵乘以系数a(θ)22,因此当h等于矩阵的最小特征值对应的特征向量时,取最小值,且该最小值等于矩阵的最小特征值。定义该矩阵为E,则:

由于,因此:

其中。

经过类似参考文献[6]附录中的推导过程可以得到结论:,进而有:

因此,该矩阵的最小特征值为:

最终的空间谱函数为:P(θ)=γ-1min。

通过搜索P(θ)的峰值可以得到信号波达方向的估计值,由于:

因此γmin≥0。

同时,观察YEC(α,τ)的特征分解式(16)不难看出:

所以,当时,有:

进而P(θ)→+∞,即空间谱函数P(θ)的峰值出现在满足式的θk,k=1,2,…,K处,即:

其中ρ为幅度调整系数。

由阵列接收数据的循环互相关协方差矩阵的特征值分解得到的估计值之后,可用如下罚函数确定λk:

经过一系列的数学推导[9],可以由以上罚函数得到互耦条件下对第k个窄带循环平稳信号的入射角度θk的估计误差为:

4 仿真

为了验证以上理论分析结果的正确性,在以下的实验中均假设一个远场BPSK信号入射到均匀线性阵列上,并假设窄带接收阵列由多个细线偶极子天线沿x轴排列而成。信号带宽与载波频率之比为0.02,可近似看作窄带信号,但为了减小非零带宽给窄带循环平稳信号测向模型所带来的误差,测向方法选择Improved Cyclic MUSIC及其扩展形式。阵列的阻抗由开环计算方法得到[10],细线偶极子天线之间的阻抗参考文献[11]进行计算。

(1)实验1:假设一个相对带宽为2%的BPSK信号入射到8元均匀线阵上,相邻阵元间距等于半波长,仿真得到C-MUSIC方法、EC-MUSIC方法的角度估计偏差和由式(15)、式(29)给出的均匀线阵的理论偏差的对比情况如图1所示。

(2)实验2:在实验1的基础上保持信号入射方向为20°,相邻阵元间距从0.1倍波长到半波长之间变化,仿真得到C-MUSIC方法、EC-MUSIC方法的角度估计偏差和由式(15)、式(29)给出的均匀线阵的理论偏差的对比情况如图2所示。

上述仿真结果表明,在阵列互耦效应存在的情况下,本文所得到的窄带循环平稳信号阵列测向方法及其共轭扩展方法的测向偏差的理论结果都具有较高的准确度,很好地反映了两种方法的测向误差随信号入射方向和阵列结构的变化情况。

5 结束语

本文通过分析互耦条件下窄带循环平稳信号测向方法C-MUSIC以及EC-MUSIC的伪数据矩阵模型,借助数学分析得到了这两种方法的测向偏差受互耦效应影响情况的解析表达式,并借助仿真实验验证了所得理论偏差的准确性。该解析结果对阵列误差条件下窄带循环平稳信号阵列测向方法的性能分析以及实际系统的误差校正都具有较强的指导意义。

摘要：常规阵列各阵元感应电流间的耦合作用是很难消除的,互耦效应的存在会导致窄带循环平稳信号阵列测向方法产生角度估计偏差。为了分析这一偏差,首先基于互耦条件下的阵列接收数据建立了用于窄带循环平稳信号测向的阵列伪数据矩阵模型,并以此为基础得到了均匀线阵和一般线阵中测向误差的解析表达式。仿真结果表明该理论偏差在各种信号环境和阵列结构中都具有很高的准确性。

关键词：阵列信号处理,波达方向估计,性能分析,阵列互耦,循环平稳信号

参考文献

[1]GARDNER W A.Signal interception:A unifying theoreti cal framework for feature detection[J].IEEE Transactions on Communications,1988,36(8):897-906.

[2]GARDNER W A,NAPOLITANO A,PAURA L.Cyclostationarity:Half a century of research[J].Signal Processing,2006,86(4):639-697.

[3]GARDNER W A.Simplification of MUSIC and ESPRIT by exploitation of cyclostationarity[J].Proceedings of the IEEE,1988,76(7):845-847.

[4]SCHELL S V,GARDNER W A.The Cramer-Rao lower bound for directions of arrival of Gaussian cyclostationary signals[J].IEEE Transactions on Information Theory,1992,38(4):1418-1422.

[5]SCHELL S V.Performance analysis of the cyclic MUSICmethod of direction estimation for cyclostationary signals[J]IEEE Transactions on Signal Processing,1994,42(11):3042-3050.

[6]CHARGE P,Wang Yide,SAILLARD J.An extended cyclic MUSIC algorithm[J].IEEE Transactions on Signal Processing,2003,51(7):1695-1701.

[7]YAN H,FAN H.On improvements of cyclic MU SIC[J].Eurasip Journal on Applied Signal Processing,2005,50(1)61-68.

[8]YEH C,LEOU M,UCCI D.Bearing estimations with mutua coupling present[J].IEEE Trans.on Antennas and Propagation,1989,37(10):1332-1335.

[9]LIU Z,HUANG Z,ZHOU Y.Bias analysis of MUSIC in the presence of mutual coupling[J].IET Signal Process ing,2009,3(1):74-84.

[10]GUPTA I,KSIENSKI A.Effect of mutual coupling on the performance of adaptive arrays[J].IEEE Trans.An tennas Propagat.,1983,9(31):785-791.

循环平稳 篇6

活塞销磨损故障是柴油发动机常见的机械故障之一, 主要表现为异响或振动异常, 严重时还可能引起拉缸乃至发动机停机, 因此必须作为发动机状态监测的重点[1]。活塞销振动是指活塞销冲击活塞座孔或连杆小头座孔而产生的振动, 当活塞销和座孔之间的配合间隙由于磨损而逐渐增大时, 往往会产生异响或异常振动, 配合间隙越大, 活塞销产生的冲击力越大, 异响或异常振动会越明显。活塞销的冲击不能被直接测量, 而是通过活塞、活塞环传向缸体。利用加速度传感器拾取的发动机缸体振动信号包含了与活塞销技术状况相对应的信息, 但是缸体振动信号是由内部构件旋转性运动和往复性运动传递到缸体的合成, 有用的信息常常淹没在强大的背景噪声中, 不利于特征的提取。

工程上存在一类统计特征函数随时间呈现周期性的非平稳信号, 这类信号被称为循环平稳信号。高阶循环平稳理论构筑在循环平稳信号的高阶循环统计量的基础上, 吸取了高阶统计量和循环平稳两者的优点, 理论上可以完全抑制任何高斯和非高斯噪声以及非平稳的高斯噪声的干扰[2], 因此是一种理想的信号特征表示方法。目前研究较多的是循环平稳信号的一阶 (均值) 、二阶 (相关函数) 循环统计特性, 而高阶循环统计量的研究还很少, 仅在通信[3]、旋转机械[4,5]等领域有所应用, 在内燃机振动方面的应用则鲜有报道。

柴油机稳速时振动信号常常表现为统计特性的周期性, 可以证明是典型的循环平稳信号[6], 而在加速过程中, 振幅逐渐增大, 各缸做功间隙时间变短, 但是发动机内部各配合副的固有振动频率不会改变, 较短的时间间隔内仍可认为是循环平稳信号。目前活塞销故障诊断的研究方法主要有虚拟样机仿真[7]、轴心轨迹计算[8]、时频分析[9]、双谱分析[10]等方法, 尚未见到利用高阶循环平稳 (higher order cyclostationary, HOCS) 理论分析活塞销振动信号的研究报道。本文尝试利用高阶循环累积量谱理论处理活塞销加速振动信号, 分析了三种不同的活塞销磨损间隙下的振动信号的循环双谱, 讨论了发动机转速对振动信号循环双谱的影响, 有效地提取出了能够反映活塞销故障状况的特征值。

1 信号的高阶循环平稳及循环双谱

1.1 高阶循环平稳

高阶循环平稳理论是以循环平稳信号的高阶统计量为理论基础的, 包括高阶循环矩、高阶循环累积量以及对应的高阶循环矩谱、高阶循环累积量谱。

k阶循环矩定义为循环平稳信号x (t) 对时间t的k阶延迟积的正弦抽取运算[8]:

Mαk x (τ1, τ2, …, τk-1) =〈x (t) x (t+τ1) …

x (t+τk-1) exp (-i2π α t) 〉t (1)

式中, α为循环频率;〈·〉t表示时间平均;τ1, τ2, …, τk-1为时间延迟。

定义使Mαk x (τ1, τ2, …, τk-1) ≠0的频率α为信号x (t) 的循环频率, 显然α=0对应于信号的平稳部分。信号的循环平稳性可以通过循环矩来判断, 如果存在不为零的循环频率, 使Mαk x (τ1, τ2, …, τk-1) ≠0, 则可以判定该信号是k阶循环平稳的。

循环平稳信号x (t) 的循环累积量可以通过信号的低阶循环矩求得:

Cαk x (τ1, τ2, …, τk-1) =

${\displaystyle \sum _{{\displaystyle \underset{p=1}{\overset{q}{\cup }}}{\rm I}{}_p={\rm I}}}$$[(-1){}^q(q-1)!{\displaystyle \sum _{\alpha {}_1+\cdots +\alpha {}_q=\beta }{\displaystyle \prod _{p=1}^q{\rm M}}}{}_{n{}_p,x}^{\alpha {}_p}(\tau {}_{{\rm I}{}_p})](2)$

式中, ${\displaystyle \sum _{{\displaystyle \underset{p=1}{\overset{q}{\cup }}}{\rm I}{}_p={\rm I}}(}\cdot )$表示在指数符集I={0, 1, …, k-1}的所有无交连的非空集合分割 (1≤q≤k) 内的求和;β为在滞后集${\displaystyle \prod _{j\in {\rm I}{}_p}x}(t+\tau {}_j)$内的纯k阶正弦波的频率;np为分割Ip中的元素个数, np=|Ip|。

循环平稳信号x (t) 的k阶循环矩谱Pαkx (f1, f2, …, fk-1) 和k阶循环累积量谱Sαkx (f1, f2, …, fk-1) 分别定义为

$\begin{array}{l}{\rm P}{}_{kx}^{\alpha }(f{}_1,f{}_2,\cdots ,f{}_{k-1})=\\ {\displaystyle \sum _{\tau {}_1=-\infty }^{\infty }\cdots }{\displaystyle \sum _{\tau {}_{k-1}=-\infty }^{\infty }{\rm M}}{}_{kx}^{\alpha }(\tau {}_1,\tau {}_2,\cdots ,\tau {}_{k-1})\exp (-\text{i}{\displaystyle \sum _{j=1}^{k-1}2}\text{π}f{}_j\tau {}_j)\end{array}$

$\begin{array}{l}S{}_{kx}^{\alpha }(f{}_1,f{}_2,\cdots ,f{}_{k-1})=\\ {\displaystyle \sum _{\tau {}_1=-\infty }^{\infty }\cdots }{\displaystyle \sum _{\tau {}_{k-1}=-\infty }^{\infty }C}{}_{kx}^{\alpha }(\tau {}_1,\tau {}_2,\cdots ,\tau {}_{k-1})\exp (-\text{i}{\displaystyle \sum _{j=1}^{k-1}2}\text{π}f{}_j\tau {}_j)\end{array}$

1.2 三阶循环累积量及循环双谱

三阶循环平稳是高阶循环平稳的最低阶形式, 计算量最小, 但继承了高阶循环平稳的优良特性。任何平稳的高斯和非高斯噪声以及非平稳的高斯噪声, 其高阶循环累积量恒等于零, 而高阶循环矩却没有这样的性质, 因此从实用的角度考虑, 高阶循环统计量最常用的是三阶循环累积量、三阶循环累积量谱即循环双谱。

由式 (2) , 信号的三阶循环累积量可以通过信号的一阶、二阶、三阶循环矩求得:

Cα3x (τ1, τ2) =Mα3x (τ1, τ2) +2 (Mα1x) 3-

Mα1x[Mα2x (τ1) +Mα2x (τ2) +Mα2x (τ2-τ1) ]

循环双谱是三阶循环平稳的主要应用, 定义为

$\begin{array}{l}S{}_{3x}^{\alpha }(f{}_1,f{}_2)=\\ {\displaystyle \sum _{\tau {}_1=-\infty }^{\infty }{\displaystyle \sum _{\tau {}_2=-\infty }^{\infty }C}}{}_{kx}^{\alpha }(\tau {}_1,\tau {}_2)\exp (-\text{i}(f{}_1\tau {}_1+f{}_2\tau {}_2))\end{array}$

由定义可知, 循环双谱是含循环频率α、频率f1、f2三个参变量的多维函数, 只能表示一个循环频率下的信号特征。当信号的循环频率为已知时, 该表示方法是有效的, 但信号的循环频率未知或者循环频率个数较多时, 就不能有效地检测出信号的特征[5], 柴油机振动信号就存在这样的问题。针对上述问题, 本文提出了在特定的循环频率段内, 以定步长搜索, 计算每一循环频率的循环双谱, 并将这些循环双谱累加, 以此来表征信号在该频段内多个循环频率中的特征。具体实现步骤如下:

(1) 选择特定的循环频率集{α}, 以定步长Δα搜索, 分别计算每一循环频率αj∈{α}的三阶循环累积量Cαj3x (τ1, τ2) ;

(2) 将循环频率集{α}内的每一循环频率αj的三阶循环累积量Cαj3x (τ1, τ2) 进行累加, 得到$C{}_{3x}^{\{\alpha \}}(\tau {}_1,\tau {}_2)={\displaystyle \sum _{j}C}{}_{3x}^{\alpha {}_j}(\tau {}_1,\tau {}_2)$;

(3) 计算C{α}3x (τ1, τ2) 的二维Fourier变换, 得到信号在循环频率区间 (α1, α2) 的循环双谱:

$\begin{array}{l}S{}_{3x}^{\{\alpha \}}(f{}_1,f{}_2)=\\ {\displaystyle \sum _{\tau {}_1}{\displaystyle \sum _{\tau {}_2}C}}{}_{3x}^{\{\alpha \}}(\tau {}_1,\tau {}_2)\exp (-\text{i}2\text{π}(f{}_1\tau {}_1+f{}_2\tau {}_2))\end{array}$

2 数值仿真实验

发动机振动信号含有多个调制源, 而相应部件的故障信号往往隐含在调制信号中。但是当调制信号较弱而被其他信号淹没时, 传统的解调方法难以奏效。而循环双谱不仅对噪声免疫, 且能够识别微弱的调制信号。设调制信号为

y (t) = (1+cos (2πfat) ) cos (2πfbt) +sin (2πfct)

y′ (t) =y (t) +r (t)

其中, 调制频率fa=30Hz;载波频率fb=100Hz;频率fc=20Hz;r (t) 为均值为0、方差为1的白噪声;y′ (t) 为调制信号y (t) 受噪声污染后的信号。采样频率fs=600Hz, 采样点数N=1024。仿真信号y (t) 及y′ (t) 的循环双谱图见图1。

由图1a可以清晰地观察到, 在 (100Hz, 100Hz) 、 (20Hz, 20Hz) 处存在谱峰, 还可以在双频率平面中得到 (100±30) Hz的频率成分, 且这些频率成分之间存在着频率的相互作用。对比图1b和图1a, 可以获得同样的频率成分, 这说明循环双谱不仅能够清晰地反映信号的特征, 而且具有较强的抑制噪声的能力。

与双谱不同, 在循环双谱中不存在频率耦合现象, 只存在频率的相互作用。如上述算例的特征频率为20Hz和100Hz, 则在循环双谱图中 (20Hz, 20Hz) 、 (20Hz, 100Hz) 、 (100Hz, 20Hz) 以及 (100Hz, 100Hz) 处有明显的峰值。正因为循环双谱中特征频率存在直接的相互作用, 因而使得循环双谱图能够比较直观地表示信号的特征频率。

3 基于循环双谱的发动机活塞销加速振动信号分析

3.1 活塞销加速振动信号的采集

试验对象为东风EQ6BT型六缸四冲程柴油发动机。试验时将第三缸活塞销与座孔的配合间隙g设置为0.01mm、0.10mm、0.20mm, 分别模拟活塞销正常、中度磨损以及严重磨损三种工况。将振动加速度传感器放置在第三缸缸体右侧与油底壳结合处, 根据先验知识, 该位置测取的振动信号对活塞销磨损故障较为敏感, 可以认为是最佳测试位置[1]。

发动机处于加速状态时, 内部运动件会产生更强烈的激励, 使故障暴露得更明显, 这是汽车维修专家所形成的共识。因此本文采用定转速触发采样方式采集加速振动信号, 即发动机在加速过程中, 如果达到预先设好的触发转速时, 采集器开始工作。试验时采集器触发转速设为1800r/min, 对应发动机中速状态。采集器采样频率设为25 600Hz, 采样点数为16 384。试验过程中, 忽略负载的影响, 发动机处于空载状态, 实际故障检测时将实车上的离合器置于分离状态。

3.2 基于循环双谱的发动机活塞销加速振动信号分析

计算采集到的发动机振动信号的循环双谱。当发动机转速为1800r/min时, 发动机曲轴旋转频率f0=30Hz, 在较短的时间内发动机的速度波动可以忽略不计, 故以f0为频率间隔, 重点考察循环频率集{α|α=kf0, k=1, 2, 3, …}内振动信号的循环双谱的特征, 即计算所有循环频率集内以曲轴旋转频率f0及其倍频为循环频率的循环双谱值并累加。三种状况下时域波形如图2所示, 其对应的的循环双谱图见图3。

由图2可知, 随着活塞销与座孔的配合间隙g的增大, 发动机振幅A也在增大, 但不易于提取相关信息。在图3中, 当活塞销正常时, 发动机振动信号的循环双谱图频率成分比较少, 随着磨损状况的加重, 配合间隙增大, 循环双谱图中频率成分能量随配合间隙的增大而增加, 但是峰值较多, 因此需要在整个双频率区间进行搜索, 找到与活塞销磨损状况相对应的故障特征值。

从循环双谱图可以看出, 大量的频率成分遍布在整个双频率区间, 因此采用常用的对角线切片搜索方法可能会丢失对角线以外的重要信息, 且循环双谱图中频率峰值太多, 故典型峰值搜索法也不便使用。本文提出一种全局搜索策略:在循环双谱图中沿f2方向的频带宽度固定取为100Hz, 即该方向上频带区间为[f′2-50, f′2+50] (Hz) , 沿f1方向进行搜索, 当搜索到频带区间为[f1a, f1b]时, 可以得到一个矩形平面, 其左下角的坐标为 (f1a, f′2-50) (Hz) , 右上角坐标为 (f1b, f′2+50) (Hz) 。将该平面内循环双谱值进行累加, 如果该平面内的平均累加值能够反映活塞销磨损故障的变化规律, 则认为该平面是循环双谱特征频率平面。每次采样时, 由于发动机转速的波动性, 特征频率峰值会发生偏移, 故提取特征频率平面比提取特征频率峰值更具有一般性和稳定性。

按照上述全局搜索策略, 搜索图3所示的活塞销三种磨损间隙下的振动信号的循环双谱图的特征频率平面。为了防止遗漏, 将f2方向上相邻的频带设置为50%的重叠, 搜索得到的前5个最明显的特征频率平面如图4所示, 其对应的平面内循环双谱幅值的平均累加值变化规律如表1所示。可以看出, 搜索得到的前5个特征频率平面内的循环双谱幅值的平均累加值随活塞销磨损间隙的增大而显著增大, 能够反映活塞销配合间隙的变化规律, 因此可以作为诊断活塞销磨损间隙的特征量。

4 结论

(1) 当信号的循环频率未知或循环频率太多时, 可在特定的循环频率频段内以定步长进行搜索, 计算所有循环频率的循环双谱并累加, 以此估计信号的循环双谱。

(2) 柴油发动机活塞销振动信号具有高阶平稳性, 利用高阶循环平稳理论分析振动信号, 能够消除噪声的影响, 可以提取出隐含在振动信号中且能够反映柴油发动机活塞销磨损状况的信息。

(3) 活塞销振动信号的循环双谱图中频率成分分布在整个频率区间内, 利用全局搜索策略, 获得5个特征频率平面, 该特征频率平面内的循环双谱幅值的平均累加值随活塞销磨损增大而明显增大, 可以将其作为对活塞销磨损故障进行诊断的特征值。

(4) 与双谱等高阶累积量谱相比, 高阶循环累积量谱计算量巨大, 不适合于故障诊断的在线应

用, 因此必须寻找高阶循环累积量谱的简化算法或快速算法, 以弥补目前算法的不足, 使之适合工程在线应用。

摘要：分析了柴油发动机活塞销振动信号的高阶循环平稳性, 提出了一种基于高阶循环平稳理论提取柴油发动机活塞销磨损故障特征的新方法。利用循环双谱分析活塞销振动信号, 考察了振动信号循环双谱图中频率成分的分布情况。定义循环双谱图中幅值的平均累加值与发动机活塞销磨损间隙变化相对应的区域为特征频率平面, 提出一种全局搜索策略以搜索特征频率平面。仿真结果及试验结果表明了循环双谱的有效性。利用所述方法在实际活塞销振动信号的循环双谱图中提取的5个特征频率平面内幅值的平均累加值可以作为活塞销磨损故障诊断的特征值。

关键词：高阶循环平稳,循环双谱,柴油发动机,活塞销,故障诊断

参考文献

[1]肖云魁.汽车故障诊断学[M].北京:北京理工大学出版社, 2006.

[2]Spooner C M.Theory and Application of Higher-order Cyclostationarity[D].Berkeley:University ofCalifornia, 1992.

[3]冯祥, 李建东.基于高阶循环累积量的SQAM信号调制识别算法[J].电子与信息学报, 2007, 29 (1) :125-128.

[4]Zhu Z K, Feng Z H, Kong F R.CyclostationarityAnalysis for Gearbox Condition Monitoring:Ap-proaches and Effectiveness[J].Mechanical Systemand Signal Processing, 2005, 19 (2) :467-482.

[5]朱忠奎, 孔凡让, 王建平, 等.循环双谱及其在齿轮箱故障识别中的应用研究[J].振动工程学报, 2004, 17 (2) :224-227.

[6]Antoni J, Daniere J, Guillet F.Effective VibrationAnalysis of IC Engines Using Cyclostationarity.PartI:a Methodology for Condition Monitoring[J].Jour-nal of Sound and Vibration, 2002, 257 (5) :815-837.

[7]白欣, 黄博, 王立臣.基于虚拟样机的发动机活塞销故障诊断研究[J].汽车技术, 2005 (4) :20-23.

[8]丁恒鉴, 肖云魁, 阳建成, 等.汽油机活塞销轴心轨迹计算及在故障诊断中的应用[J].内燃机学报, 2001, 19 (2) :187-191.

[9]郑海波, 李志远, 陈心昭, 等.基于时频分布的发动机异响特征分析及故障诊断研究[J].内燃机学报, 2002, 20 (3) :267-272.