CN101256094A - 提取特定周期信号成分的信号消噪方法 - Google Patents

Abstract

本发明公开了一种提取信号中特定周期成分的信号消噪方法，该方法采用非递归滤波器技术构造梳状滤波器；用频域内插技术精确获得滤波器的中心频率；能有效消除信号中的各种噪声的干扰，提取出信号中感兴趣的周期成分。具备时域同步平均技术的优点，且不需要时标信号，需要的信号量小，计算速度快，适用范围广，实施成本低。可应用于机械动态信号分析和故障诊断领域，以及信号消噪、信号处理和特征提取等相关应用领域。

Description

提取特定周期信号成分的信号消噪方法

技术领域

本发明属于机械动态信号分析和故障诊断领域，进一步可扩展到信号消噪、信号处理和特征提取等相关应用领域，特别涉及一种提取周期信号成分的信号消噪方法。

背景技术

在自然界和工程应用领域，所涉及的信号往往受到各种确定的和随机的干扰源作用，使得信号中的有用成分被污染或淹没，严重影响了信号的分析精度，因此信号消噪在信号分析、检测和识别领域具有重要的应用价值。

在机械故障分析和诊断中，旋转和往复机械运行过程的各种动态信号是随着机器的运转而周期性重复的，其频率是机器旋转频率的整倍数，这种周期性信号特征反映了机器的运行状态。然而，在强噪声干扰的情况下，有用信息被削弱，不仅信号的时间历程显示不出规律性，而且由于常用的谱分析不能略去任何输入分量，在频谱图中这些周期成分很可能被淹没在噪声背景中。目前，具有独特优点的时域同步平均技术则是解决这一问题的主要手段。

时域同步平均(参见罗德阳.时域同步平均原理与应用.振动、测试与诊断，1999，19(3))是从混有噪声干扰的复杂周期信号中提取感兴趣的周期性成分的过程，也称相干检波(coherent detction)。方法的要点是对测得的振动信号以一定的周期为间隔截取若干段，然后将所截取的信号段中对应的离散点相加后取算术平均。这样可以消除原信号中的干扰噪声而保留指定的周期成分。

时域同步平均在信号处理和分析中具有如下的独特优点：

①、时域同步平均可以有效消除和抑制信号中的干扰噪声，而将感兴趣的(有用的)周期成分完整地保留下来。其实质是将信号中感兴趣的周期成分的基频及其倍频分量增强，且基频和各倍频之间的相位关系保持不变。这一优点能够将信号中感兴趣的周期成分完整地展现出来，从而揭示信号的周期性变化规律，籍此可以发现和分析信号中包含的重要信息。

②、时域同步平均对信号中周期性成分的保留和增强具有选择性，也就是说，它不仅可以消除随机干扰，还对不感兴趣的(无用的)其他周期成分有抑制和消除作用。这一优点一方面可以将分析的重点集中于感兴趣的周期成分，而排除其他周期成分的干扰。另一方面，也可以将不同的周期成分分别提取和展示出来。

正是由于时域同步平均具有上述两个独特的优点，使其在信号分析领域，尤其是机械信号分析和故障诊断领域，具有十分重要的应用价值，成为信号预处理不可或缺的重要手段，而且至今还没有有效的替代方法和技术出现。

然而，尽管时域同步平均有难以替代的独特优点，但也有其自身的局限性，常常由于环境和条件的制约而无法付诸实施。主要存在以下问题：

①、与通常的信号采集不同，时域同步平均不仅要拾取被分析信号，同时还要拾取旋转轴的时标脉冲(键相信号)，来锁定各信号段的起始点。这需要安装专门的传感器和相应的硬件装置，无疑会增加额外的费用，有时还会由于现场环境的制约而无法实现。这就给时域同步平均带来了具体实施上的难度，使其适用范围受到限制。

②、因为A/D变换器的采样频率一经设定是不变的，而且机器转速的微小波动使得时标脉冲频率及周期也随时在变化，常规采样不可能保证各数据段点数相等，这会对时域同步平均效果产生严重的影响(参见Blunt D M.Synchronous averaging of helicopter tail gearbox vibration：phase referenceconsiderations[R].Defense Science and Technology Organization of AustraliaDOD，DSTOTR20397，1998。)。解决这一问题的方法是采用频率跟踪技术，或等转角触发采样，使实际的采样频率实时跟踪回转频率。这无疑又增加了时域同步平均的花费和实施难度。

③、时域同步平均是通过对截取的多段信号经平均运算实现的，因此，从原理上讲，截取的段数越多，干扰噪声的衰减越大，平均的效果越好。但是，从另一方面来看，截取的段数越多，需要存储的信号容量越大，计算时间越长。

发明内容

为了克服现有时域同步平均方法的上述缺陷或不足，本发明的目的在于，提供一种简单、实用、适应范围广的提取特定周期信号成分的信号消噪方法。

为了实现上述任务，本发明采取如下的技术解决方案：

一种提取特定周期信号成分的信号消噪方法，其特征在于，具体包括下列步骤：

步骤一，信号采集：

选取固定的采样时间间隔T，按常规采样将来自传感器的模拟信号转换为离散数字信号，即，假定来自传感器的模拟信号为x(t)，通过A/D转换器将x(t)转换为离散数字信号x_n，使得：

x_n＝x(nT)，n＝0，1，2，...，N-1

其中，N为样本长度，取N＝1024；

为保证提取的周期成分中包含足够的倍频数，采样时间间隔T按以下规则取值：

$T\≤\frac{1}{2.5K{f}_1}$

式中f₁为要提取的周期成分的基频；K为倍频的数目，通常为了保证提取的周期成分更为真实，K取值越大越好。根据经验，K的取值不小于5，即可满足工程要求；

步骤二，确定周期成分的基频：

对离散信号x_n，用快速富立叶变换获得离散频谱X_n，进而求得幅值谱|X_n|，n＝0，1，2，...，N-1。若幅值谱中周期成分基频附近的最大谱线为|X_最大|，其对应的频率可作为基频的近似值

采用频域内插技术对进行校正，以获得f₁的精确值，并按下列方法进行：

设|X_次大|为与|X_最大|两侧的左右两根谱线中最大的一根谱线，则可通过下式求得基频的频率精确值，即

$f_1={\hat{f}}_1+\mathrm{\δ}\·\mathrm{\Δf}$

式中，Δf为离散频谱的频率间隔，由下式表示：

$\mathrm{\Δf}=\frac{1}{N\·T}$

δ为频率校正系数，由下式计算：

上式中，如果|X_次大|位于|X_最大|的左侧，取负号；否则，如果|X_次大|位于|X_最大|的右侧，取正号；

步骤三，确定周期成分的倍频：

确定了周期成分的基频以后，通过下式确定周期成分的倍频：

f_k＝kf₁，k＝2，3，...，K

式中，K为倍频的数目；

步骤四，构造梳状滤波器：

为了保证提取的周期成分的波形特征，要求滤波器具备零相移特性，采用非递归数字滤波器技术构造梳状滤波器，为了消除吉布斯效应的影响，对滤波器的冲击响应进行加汉宁窗处理；

获得的冲击响应函数由下式给出：

$b_m=\left(\frac{\sin \left(2\mathrm{\π}(f+f_b)\mathrm{mT}\right)-\sin \left(2\mathrm{\π}(f-f_b)\mathrm{mT}\right)}{\mathrm{\πmT}}\right)(1+\cos (\mathrm{\πm}/M))$

式中，f为带通滤波器的中心频率；f_b为滤波器通带的半带宽，M为滤波器阶数，且M取值大于40即可满足要求；

用上述基本的带通滤波器以及步骤三获得的信号周期成分的基频和倍频f_k＝kf₁，k＝1，2，...，K，构造所需的梳状滤波器的单位冲击响应为：

$b_m=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\left(\frac{\sin \left(2\mathrm{\π}(f_k+f_b)\mathrm{mT}\right)-\sin \left(2\mathrm{\π}(f_k-f_b)\mathrm{mT}\right)}{\mathrm{\πmT}}\right)(1+\cos (\mathrm{\πm}/M))$

式中，K为信号中要提取的周期成分的倍频数目；

步骤五，周期成分的提取：

构造了梳状滤波器的冲击响应之后，按以下非递归过程完成对信号x_n的处理，即：

$y_n=\underset{m=-M}{\overset{M}{\mathrm{\Σ}}}{b}_m{x}_{n-m},$ n＝0，1，2，...，N-1

其中，y_n即为滤波后的离散信号，也就是经过消除噪声干扰后从信号x_n中提取出来的特定周期成分。

本发明的提取特定周期信号成分的信号消噪方法与传统的时域同步平均技术相比，不仅可以实现时域同步平均的功能，具备时域同步平均的优点

(参见背景技术中的优点①和优点②)而且扩展了时域同步平均的功能和适应范围，具有以下独特的技术特点：

①.时域同步平均需要时标信号(键相信号)，以便锁定各信号段的起始点。而本发明不需时标信号，因此省去了安装专门的传感器和相应的硬件装置的花费，具有实施成本低的特点；同时，也不受安装时标传感器的现场环境和条件的制约，拓宽了适用范围。

②.时域同步平均需要截取多段信号，然后对各段信号进行平均运算，因此，需要的信号量和存储空间大，计算时间长。而本发明只需一段信号，存储容量小，计算速度快。

③.此外，本发明可以根据需要灵活构造梳状滤波器，可以选取基频和倍频中的任意组合来构造梳状滤波器。例如，选取奇数频率(即，选取基频和3、5、7...倍频为中心频率)或偶数频率(即，选取2、4、6...倍频为中心频率)等等。甚至可以构造不等距梳状滤波器，例如，选取基频、2、4、5等等。这一优点，在机械信号分析和故障诊断中常常是特别有用的。

本发明的方法采用非递归滤波器技术构造梳状滤波器；用频域内插技术精确获得滤波器的中心频率；能有效消除信号中的各种噪声的干扰，提取出信号中感兴趣的周期成分。具备时域同步平均技术的优点，且不需要时标信号，需要的信号量小，计算速度快，适用范围广，实施成本低。可应用于机械动态信号分析和故障诊断领域，以及信号消噪、信号处理和特征提取等相关应用领域。

附图说明

图1为实施例中齿轮振动信号的原始波形图；

图2为图1信号的频谱；

图3为通过本发明获得的梳状滤波器的频率响应函数；

图4为用本发明提取的特定周期成分的波形图。

为了更清楚地理解本发明，以下结合附图和发明人依本发明的技术方案给出的实施例对本发明作进一步地详细描述。

具体实施方式

一、本发明的技术方案构思基础及要点

时域同步平均之所以能够保留感兴趣的完整周期信号成分而将干扰噪声和其它周期成分去除，原因在于同步平均等价于一种梳状数字滤波器。该滤波器由一系列等距分布的带通滤波器及旁瓣组成，带通滤波器的中心频率是特定周期成分的频率的整数倍，并能保证通过的信号各频率之间无相移。此外，平均次数越大，通带宽度越窄，消除干扰噪声的效果越好。时域同步平均的这种梳状滤波器特性，正是本发明技术方案的技术构思基础。

本发明的核心在于利用数字滤波器技术构造梳状滤波器，对信号进行滤波处理，从而达到提取感兴趣的周期成分、滤除随机噪声和其它不感兴趣的周期成分的目的。要点在于：

①、梳状滤波器的构造：为保证特定周期成分的基频及其倍频之间无相移，首先采用有限冲击响应滤波器技术设计出基本的零相移带通滤波器。该滤波器的带宽应尽可能地窄，以确保对提取的周期成分有最好的效果。然后根据信号中要提取的周期成分的基频及其倍频，构造出梳状滤波器。图3给出了构造的梳状滤波器频率特性的一个实例。

②、梳状滤波器各中心频率的确定：首先需要事先根据要提取的周期成分的基频及其倍频来确定梳状滤波器各中心频率的值，此外，要点①中要求带通滤波器的带宽尽可能地窄，因此精确确定信号中周期成分的基频及其倍频十分重要。本发明通过用快速富立叶变换(FFT)求解信号的频谱来获得信号周期成分的基频及其倍频值。但是，由于频率分辨率的制约，由FFT获得的频率值，常常存在很大误差，因此本发明采用频域内插方法提高基频和倍频的估值精度。

二.具体实现步骤

<1>信号采集

选取固定的采样时间间隔T，按常规采样将来自传感器的模拟信号转换为离散数字信号。即，假定来自传感器的模拟信号为x(t)，通过A/D转换器将x(t)转换为离散数字信号x_n

x_n＝x(nT)，n＝0，1，2，...，N-1

其中，N为样本长度，通常取为1024。

为保证提取的周期成分中包含足够的倍频数，采样时间间隔T按以下规则取值：

$T\&le;\frac{1}{2.5K{f}_1}$

式中f₁为要提取的周期成分的基频(1/每秒)；K为倍频的数目，通常为了保证提取的周期成分更为真实，K取值越大越好。根据经验，K的取值不小于5，即可满足工程要求。

<2>确定周期成分的基频

对离散信号x_n，用快速富立叶变换(FFT)获得离散频谱X_n，进而求得幅值谱|X_n|，n＝0，1，2，...，N-1。若幅值谱中周期成分基频附近的最大谱线为|X_最大|，其对应的频率可作为基频的近似值

通常，对于离散频谱，由于频率分辨率的制约，

与真实值f₁并不相等，存在很大偏差。为此，本发明采用频域内插技术对

进行校正，以获得f₁的精确值。

设|X_次大|为与|X_最大|两侧的左右两根谱线中最大的一根谱线，则可通过下式求得基频的频率精确值，即

$f_1={\hat{f}}_1+\mathrm{\&delta;}\&CenterDot;\mathrm{\&Delta;f}$

式中，Δf为离散频谱的频率间隔，由下式表示：

$\mathrm{\&Delta;f}=\frac{1}{N\&CenterDot;T}$

δ为频率校正系数，由下式计算：

上式中，如果|X_次大|位于|X_最大|的左侧，取负号；否则，如果|X_次大|位于|X_最大|的右侧，取正号。

<3>确定周期成分的倍频

确定了周期成分的基频以后，可以通过下式确定其倍频：

f_k＝kf₁，k＝2，3，...，K

式中，K为倍频的数目。

<4>构造梳状滤波器

为了保证提取的周期成分的波形特征，要求滤波器必须具备零相移特性，因此采用非递归数字滤波器技术构造梳状滤波器。首先设计一个基本的具有理想通带特性的带通滤波器。为了消除吉布斯效应的影响，对滤波器的冲击响应进行加汉宁(Hanning)窗处理。获得的冲击响应函数由下式给出：

$b_m=\left(\frac{\sin \left(2\mathrm{\&pi;}(f+f_b)\mathrm{mT}\right)-\sin \left(2\mathrm{\&pi;}(f-f_b)\mathrm{mT}\right)}{\mathrm{\&pi;mT}}\right)(1+\cos (\mathrm{\&pi;m}/M))$

式中，f为带通滤波器的中心频率；f_b为滤波器通带的半带宽，M为滤波器阶数，通常滤波器阶数M取得越大，精度越高，但计算量越大。根据经验，M取为大于40，即可满足工程应用要求。

用上述基本的带通滤波器以及步骤<3>获得的信号周期成分的基频和倍频f_k＝kf₁，k＝1，2，...，K，可以方便地构造所需的梳状滤波器的单位冲击响应，即

$b_m=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\left(\frac{\sin \left(2\mathrm{\&pi;}(f_k+f_b)\mathrm{mT}\right)-\sin \left(2\mathrm{\&pi;}(f_k-f_b)\mathrm{mT}\right)}{\mathrm{\&pi;mT}}\right)(1+\cos (\mathrm{\&pi;m}/M))$

式中，K为信号中要提取的周期成分的倍频数目。

值得指出的是，上述由基本的带通滤波器构造梳状滤波器的过程是非常简单、灵活的，这意味着，可以根据应用需要，构造不同用途的梳状滤波器。例如，申请人可以只取信号周期成分中的奇数倍频来构造梳状滤波器，获得仅有奇数倍频的信号波形变化规律，这一点在机械故障诊断中是非常有价值的，这也是本发明的特点之一。

<5>周期成分的提取

构造了梳状滤波器的冲击响应之后，借助于以下非递归过程完成对信号x_n的处理，即

$y_n=\underset{m=-M}{\overset{M}{\mathrm{\&Sigma;}}}{b}_m{x}_{n-m},$ n＝0，1，2，...，N-1

其中，y_n即为滤波后的离散信号，也就是经过消除噪声干扰后从信号x_n中提取出来的特定周期成分。

实施例：

此处描述的实施例是用某低速齿轮齿面磨损损伤模拟试验的振动加速度数据进行的。该齿轮的转速为n₀＝120转/分，齿轮的齿数z＝25，齿轮各齿的齿面由于长期运行而处于均匀磨损状态。由齿轮的转速和齿数可以计算出齿轮的啮合频率为

f_z＝z·n₀/60＝25×120/60＝50(Hz)

根据步骤<1>对信号采样获得离散信号x_n，采样时间间隔T＝0.001秒，样本长度N＝1024。图1是采集的原始信号时域波形图。

按步骤<2>采用FFT计算信号x_n的频谱|X_n|，参见图2。在齿轮的啮合频率f_z＝50Hz附近找出最大谱线，并通过频域内插技术求得基频f₀(此例算得的精确值为51.2Hz)。

按步骤<3>计算f₀的倍频值。由于采样频率为1/0.001＝1000Hz，有效分析频率范围为0～500Hz，因此可取倍频数目K＝9。

用计算获得的基频和倍频值，按步骤<4>构造梳状滤波器。图3是构造的梳状滤波器的频率响应函数。

从该实施例可以看出，原始信号(图1)受到强大噪声的污染，信号杂乱无章，从中无法得出有价值的信息。而经过本发明的方法处理以后(图4)，信号中周期成分明显地展显现出来，干扰噪声受到抑制和消除。可见本发明的方法有很强的消噪能力，对周期成分的提取有很好的效果。

Claims (1)

1.一种提取特定周期信号成分的信号消噪方法，其特征在于包括下列步骤：

步骤一，信号采集：

选取固定的采样时间间隔T，按常规采样将来自传感器的模拟信号转换为离散数字信号，即，假定来自传感器的模拟信号为x(t)，通过A/D转换器将x(t)转换为离散数字信号x_n，使得：

x_n＝x(nT)，n＝0，1，2，...，N-1

其中，N为样本长度，取N＝1024；

为保证提取的周期成分中包含足够的倍频数，采样时间间隔T按以下规则取值：

$T\&le;\frac{1}{2.5K{f}_1}$

式中f₁为要提取的周期成分的基频；K为倍频的数目，且K的取值不小于5；

步骤二，确定周期成分的基频：

对离散信号x_n，用快速富立叶变换获得离散频谱X_n，进而求得幅值谱|X_n|，n＝0，1，2，...，N-1。若幅值谱中周期成分基频附近的最大谱线为|X_最大|，其对应的频率可作为基频的近似值

采用频域内插技术对

进行校正，以获得f₁的精确值，并按下列方法进行：

设|X_次大|为与|X_最大|两侧的左右两根谱线中最大的一根谱线，则可通过下式求得基频的频率精确值，即

$f_1={\hat{f}}_1+\mathrm{\&delta;}\&CenterDot;\mathrm{\&Delta;f}$

式中，Δf为离散频谱的频率间隔，由下式表示：

$\mathrm{\&Delta;f}=\frac{1}{N\&CenterDot;T}$

δ为频率校正系数，由下式计算：

上式中，如果|X_次大|位于|X_最大|的左侧，取负号；否则，如果|X_次大|位于|X_最大|的右侧，取正号；

步骤三，确定周期成分的倍频：

确定了周期成分的基频以后，通过下式确定周期成分的倍频：

f_k＝kf₁，k＝2，3，...，K

式中，K为倍频的数目；

步骤四，构造梳状滤波器：

为了保证提取的周期成分的波形特征，要求滤波器具备零相移特性，采用非递归数字滤波器技术构造梳状滤波器，为了消除吉布斯效应的影响，对滤波器的冲击响应进行加汉宁窗处理；

获得的冲击响应函数由下式给出：

$b_m=\left(\frac{\sin \left(2\mathrm{\&pi;}(f+f_b)\mathrm{mT}\right)-\sin \left(2\mathrm{\&pi;}(f-f_b)\mathrm{mT}\right)}{\mathrm{\&pi;mT}}\right)(1+\cos (\mathrm{\&pi;m}/M))$

式中，f为带通滤波器的中心频率；f_b为滤波器通带的半带宽，M为滤波器阶数，且M取值大于40即可满足要求；

用上述基本的带通滤波器以及步骤三获得的信号周期成分的基频和倍频f_k＝kf₁，k＝1，2，...，K，构造所需的梳状滤波器的单位冲击响应为：

$b_m=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\left(\frac{\sin \left(2\mathrm{\&pi;}(f_k+f_b)\mathrm{mT}\right)-\sin \left(2\mathrm{\&pi;}(f_k-f_b)\mathrm{mT}\right)}{\mathrm{\&pi;mT}}\right)(1+\cos (\mathrm{\&pi;m}/M))$

式中，K为信号中要提取的周期成分的倍频数目；

步骤五，周期成分的提取：

构造了梳状滤波器的冲击响应之后，按以下非递归过程完成对信号x_n的处理，即：

$y_n=\underset{m=-M}{\overset{M}{\mathrm{\&Sigma;}}}{b}_m{x}_{n-m},$ n＝0，1，2，...，N-1

其中，y_n即为滤波后的离散信号，也就是经过消除噪声干扰后从信号x_n中提取出来的特定周期成分。

Images