CN101709997A - 一种振动信号处理的谐波窗函数 - Google Patents

Abstract

本发明公布了一种振动信号处理的谐波窗函数，本发明利用谐波小波的基础函数，构造了一个谐波窗函数，它具有自适应窗功能。在不减少信息点数的情况下，用该方法更好地实现了振动信号的加窗分析，并且由于其紧支撑特性，其在时域、频域的泄漏都低于目前任何一种其他窗函数。该方法具有灵活变化的时频分析窗口和中心调节，可以选择任意感兴趣频段进行分析。本发明方案可精确地将对复杂振动信号进行检测分析。

Description

一种振动信号处理的谐波窗函数

技术领域

本发明涉及信号后处理技术，属于一种复杂振动信号处理新方法，适用于复杂振动信号的加窗分析。

背景技术

现代的复杂振动信号大多具有非平稳、强噪声的特征。振动信号的分析，就是指从噪声中提取有用信号、研究信号特征，进一步确定所测试结构目前所处状态。复杂振动信号(诸如存在于裂纹故障诊断等检测)的非平稳性增加了对结构故障进行有效诊断的困难，而强噪声环境又淹没了诸多有用的信号。因而，要确定结构的动态故障特征，需研究更加有效的信号分析方法。目前对于复杂振动信号的分析，在理论和应用上已取得了若干成果。常规的时域分析方法、频谱分析方法和故障诊断方法能够准确分析与诊断出若干导致故障的原因，在工业界已经取得若干满意的成果。而对于非平稳的故障信号的识别，采用STFT、时频分析或小波分析方法也能得到较好的解决。STFT是加了时间窗的傅立叶谱分析，但必须假定待定分析数据是分段平稳的，另一方面为了频率精度又需窗口足够长，就限制了STFT分析的应用。小波分析方法可以成功地进行非平稳信号的分析。为满足高精度动平衡测量的要求，以提高不平衡量及其相位的提取精度为研究目标，利用谐波小波的相位锁定和窄带分析能力，设计了谐波小波自适应滤波器，发展了确保高精度提取幅值和相位的自适应方法。仿真与实际应用表明，该方法满足高精度和实时性的要求，具有充裕的应用潜力(杨克己，基于谐波小波的自适应滤波在高精度动平衡检测系统中的应用。仪器仪表学报，2005，26(8)：133-135)。对旋转机械振动信号通过经验模态分解，然后经过希尔伯特变换获得信号频谱(杨世锡，旋转机械振动信号基于EMD的希尔伯特变换和小波变换时频分析比较。中国电机工程学报，2003，23(6)：102-107)，能够同时用于线性振动信号和非线性振动信号的分析。利用基于二进的小波分解变换技术，可以把不同频段的时域信号展现出来。但小波分析不具有自适应性的特点，一旦小波被选定，就必须用它来分析所有的待分析数据。同时随着分解层数的增加，由于二抽一的二进小波分解特性，密集特征信号的分析精细程度大大折扣(李舜酩等.行驶车辆振动信号的小波分析，汽车工程，1997.19(6)：370-375)。再一方面，二进小波分解在进行傅立叶变换后，相位发生变化，必须进行保相处理(李舜酩.二进小波振动检测中低频保相的二次变换处理.仪器仪表学报，2005，26(4)：370-373)。谐波小波分析具有频域相位锁定的能力，具有正交性、局部分解性及自适应性。因此，综合各种信号处理方法，本发明利用谐波小波的优良特性，构造一个更灵活的谐波窗函数。该窗具有谐波小波极好的紧支特性，同时也具有极其优秀抑制频率泄漏的能力，可以在任意感兴趣的分析频段展开分析，特别适合对复杂振动信号的加窗分析。

发明内容

本发明的目的就是在综合谐波小波的优点，构造一个谐波窗函数。该窗函数具有的自适应特性，具有灵活的变化窗口和移动中心，能优秀地抑制频率泄漏，对振动信号在任何感兴趣的频段进行加窗分析。

本发明为实现上述目的，采用如下技术方案：

本发明一种振动信号处理的谐波窗函数，其特征在于包括如下步骤：

步骤101，构造谐波小波函数

谐波小波函数是由实函数ψ_e(x)和ψ_o(x)得到的傅立叶变换式：

角标e和o分别表示该实数是变量x的偶函数和奇函数；对式(1)和式(2)作逆傅立叶变换得到：

$\left\{\begin{array}{c}{\mathrm{\ψ}}_e\left(x\right)={\∫}_{-\∞}^{\∞}{\widehat{\mathrm{\ψ}}}_e\left(\mathrm{\ω}\right)\exp \left(\mathrm{i\ωx}\right)\mathrm{d\ω}\\ =[\sin \left(4\mathrm{\πx}\right)-\sin \left(2\mathrm{\πx}\right)]/\left(2\mathrm{\πx}\right)\\ {\mathrm{\ψ}}_o\left(x\right)={\∫}_{-\∞}^{\∞}{\widehat{\mathrm{\ψ}}}_o\left(\mathrm{\ω}\right)\exp \left(\mathrm{i\ωx}\right)\mathrm{d\ω}\\ =-[\cos \left(4\mathrm{\πx}\right)-\cos \left(2\mathrm{\πx}\right)]/\left(2\mathrm{\πx}\right)\end{array}\right.---\left(3\right)$

此定义的复函数：

ψ(x)＝[exp(i4πx)-exp(i2πx)]/(i2πx)(4)

为谐波小波，其中i为复数，ω为频率，exp为自然对数底数；

令m＝2^j，n＝2^j+1，j∈Z+，就得到小波变换在不同分解层上的频段分解结果：

${\mathrm{\ψ}}_{m,n}\left(x\right)=\frac{e^{\mathrm{in}2\mathrm{\πx}}-e^{\mathrm{im}2\mathrm{\πx}}}{i2\mathrm{\π}(n-m)x}---\left(5\right)$

式(5)即是谐波小波在时域上的一般表达式；

给定谐波小波位移步长k/(n-m)，k∈Z+且k≠0，Z为整数，则式(5)变为

${\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})=\frac{\{e^{\mathrm{in}2\mathrm{\π}[x-k/(n-m)]}-e^{\mathrm{im}2\mathrm{\π}[x-k/(n-m)]}\}}{i2\mathrm{\π}(n-m)[x-k/(n-m)]}---\left(6\right)$

这就是带宽为(n-m)2π、分析中心在x＝k/(n-m)的谐波小波的一般表达式；对时间离散信号f^d(r)，r＝0，1，2，...，N-1，d表示离散信号，其谐波小波变换可写为：

$W_f(m,n,k)=\frac{(n-m)}{N}\underset{r=0}{\overset{N-1}{\mathrm{\Σ}}}{f}^d\left(r\right){\stackrel{\‾}{\mathrm{\ψ}}}_{m,n}(r-\frac{k}{n-m})---\left(7\right)$

此即信号的离散谐波小波变换表达式，其中，ψ_m，n表示ψ_m，n的共轭；其频域表达式为：

$\hat{W}(m,n,\mathrm{\ω})=\hat{f}\left(\mathrm{\ω}\right){\widehat{\stackrel{\‾}{\mathrm{\ψ}}}}_{m,n}\left[(n-m)\mathrm{\ω}\right]---\left(8\right)$

步骤102，重新定义m、n的取值范围

m，n∈R+且m＜n，即m、n在正实数域内可以取非整数值，R为实数，这样在不进行任何分解的情况下，滑动窗口到所选择的频段上并伸缩窗口：

$\left.\begin{array}{c}m=\mathrm{qB}\\ n=(q+1)B\end{array}\right\}---\left(9\right)$

其中，q∈R+，B为分析频带宽度；

步骤103，构造谐波窗函数

重新定义m、n的取值后，谐波小波函数即为谐波窗函数，谐波窗函数的三角函数表达式为：

${\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})=\left\{\cos \right[2\mathrm{\πn}(x-\frac{k}{n-m})]-\cos [2\mathrm{\πm}(x-\frac{k}{n-m})\left]\right\}/\left[i2\mathrm{\π}(n-m)(x-\frac{k}{n-m})\right]$

m，n∈R+，k∈Z+且k≠0(10)

其共轭表达式为

${\stackrel{\‾}{\mathrm{\ψ}}}_{m,n}(x-\frac{k}{n-m})=-{\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})---\left(11\right)$

式(10)、(11)就是谐波窗函数。

本发明与现有方案比较，本发明具有如下优势：

①采用谐波小波基函数构造，谐波窗变换简洁、直观、易实现；

②谐波窗具有时频窗分析中心自由移动和任意改变窗口尺度的优秀特性，能够更好的处理非平稳信号和频率密集信号；

③具有及其优秀的抗泄漏能力。与目前振动信号的加窗分析中最好的汉宁窗比较，其主瓣能量更集中，旁瓣能量衰减快特点，时域、频域的抗泄漏能力更强。。

附图说明

图1：汉宁窗函数在时频域的特性图，a为汉宁窗时域指示图，图b为汉宁窗频域指示图。

图2：谐波窗函数在时频域的特性图，图c为谐波窗时域指示图，图d谐波窗频率指示图。

图3：是本发明的方法流程图。

图4：复函数

的频谱特性图。

具体实施方式

下面结合附图对发明的技术方案进行详细说明：

如图1、图2所示，图a为汉宁窗时域指示图，图b为汉宁窗频域指示图，图c为谐波窗时域指示图，图d谐波窗频率指示图。

如图3所示，本发明包括如下步骤：

步骤101，谐波小波函数构造

通常意义上(也是实现上所需要)的小波变换和小波包变换是采用隔点采样，即二抽一采样。在数学意义上，这体现了“二进”的性质。之所以采用隔点采样，是为了进行数据编码的压缩，这虽然达到无冗余存贮且可以重建信号的目的，但随着分解层数的增加，各层、各频段序列的数据点数也减半、采样频率也减半，在数据点数太少时，信号细节也存在着失真问题。

另外，二进小波变换与二进小波包变换在每一层分解时都是“一分为二”的。即每进行一次分解都是进行低通滤波和高通滤波，把上一层信号分成低频段和高频段(在二进小波变换中叫做逼近信号和细节信号)。

针对二进小波变换存在的问题，1993年Newland在研究Daubechise小波的谱特征后，结合分析“窗口”的要求，提出了一种全新的小波，叫谐波小波。Newland给出的谐波小波函数是由实函数ψ_e(x)和ψ_o(x)(角标e和o分别表示该实数是变量x的偶函数和奇函数)得到的傅立叶变换式：

对式(10)和式(11)作逆傅立叶变换，得到

$\left\{\begin{array}{c}{\mathrm{\ψ}}_e\left(x\right)={\∫}_{-\∞}^{\∞}{\widehat{\mathrm{\ψ}}}_e\left(\mathrm{\ω}\right)\exp \left(\mathrm{i\ωx}\right)\mathrm{d\ω}\\ =[\sin \left(4\mathrm{\πx}\right)-\sin \left(2\mathrm{\πx}\right)]/\left(2\mathrm{\πx}\right)\\ {\mathrm{\ψ}}_o\left(x\right)={\∫}_{-\∞}^{\∞}{\widehat{\mathrm{\ψ}}}_o\left(\mathrm{\ω}\right)\exp \left(\mathrm{i\ωx}\right)\mathrm{d\ω}\\ =-[\cos \left(4\mathrm{\πx}\right)-\cos \left(2\mathrm{\πx}\right)]/\left(2\mathrm{\πx}\right)\end{array}\right.---\left(3\right)$

此定义的复函数

ψ(x)＝[exp(i4πx)-exp(i2πx)]/(i2πx)(4)

为谐波小波。其频域波形如图4所示。可以从此图看出，

具有极好的紧支特性特征。与二进小波变换比较，在数字计算中令m＝2^j，n＝2^j+1，就得到小波变换在不同分解层上的频段分解结果：

${\mathrm{\ψ}}_{m,n}\left(x\right)=\frac{e^{\mathrm{in}2\mathrm{\πx}}-e^{\mathrm{im}2\mathrm{\πx}}}{i2\mathrm{\π}(n-m)x}---\left(5\right)$

式(5)即是谐波小波在时域上的一般表达式。由式(5)可看出，谐波小波在分解过程中不产生移相，具有相位锁定功能。

给定谐波小波位移步长k/(n-m)，则式(5)变为

${\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})=\frac{\{e^{\mathrm{in}2\mathrm{\π}[x-k/(n-m)]}-e^{\mathrm{im}2\mathrm{\π}[x-k/(n-m)]}\}}{i2\mathrm{\π}(n-m)[x-k/(n-m)]}---\left(6\right)$

这就是带宽为(n-m)2π、分析中心在x＝k/(n-m)的谐波小波的一般表达式。对时间离散信号f^d(r)，r＝0，1，2，...，N-1，其谐波小波变换可写为：

$W_f(m,n,k)=\frac{(n-m)}{N}\underset{r=0}{\overset{N-1}{\mathrm{\Σ}}}{f}^d\left(r\right){\stackrel{\‾}{\mathrm{\ψ}}}_{m,n}(r-\frac{k}{n-m})---\left(7\right)$

此即信号的离散谐波小波变换表达式。其中，ψ_m，n表示ψ_m，n的共轭。其频域表达式为

$\hat{W}(m,n,\mathrm{\ω})=\hat{f}\left(\mathrm{\ω}\right){\widehat{\stackrel{\‾}{\mathrm{\ψ}}}}_{m,n}\left[(n-m)\mathrm{\ω}\right]---\left(8\right)$

步骤102，重新定义m、n的取值范围

谐波小波的时频表达式(7)或式(8)中的m、n决定了小波变换中的层次，起着与二进小波变换中2^j中的j相同的作用。具体地说，若f_h为最高分析频率，则

n-m＝2^-jf_h  n，m，j∈Z⁺  (9)

很显然n≥m。在谐波小波变换中，可设如下条件

$\left\{\begin{array}{cc}n=2m& m\≠0\\ n=1& m=0\end{array}\right.---\left(10\right)$

由此在式(10)中令m＝2-jf_h，则n＝2m＝2^1-jf_h。

同样，若令m＝0，则由式(9)得

n＝2^-jf_h，m＝0           (11)

或者若令n＝f_h，则由式(10)得

m＝f_h(1-2^-j)             (12)

放松m、n的取值关于2^j的约束，重新定义m、n的取值：m，n∈R+且m＜n，即m、n在正实数域内可以取非整数值。这样就可以在不进行任何分解的情况下，滑动窗口到所选择的频段上并伸缩窗口：

$\left.\begin{array}{c}m=\mathrm{qB}\\ n=(q+1)B\end{array}\right\}q\∈R+---\left(13\right)$

其中，B为分析频带宽度。

步骤103，构造谐波窗函数

重新定义m、n的取值后，谐波小波函数就变成了一个新型的窗函数，这里定义为谐波窗函数。如果我们对哪个频段的信号感兴趣，可将信号分解到由式(13)确定的信号频段上。

谐波窗函数的三角函数表达式为

${\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})=\left\{\cos \right[2\mathrm{\πn}(x-\frac{k}{n-m})]-\cos [2\mathrm{\πm}(x-\frac{k}{n-m})\left]\right\}/\left[i2\mathrm{\π}(n-m)(x-\frac{k}{n-m})\right]$

m，n∈R+，k∈Z+且k≠0     (14)

其共轭表达式为

${\stackrel{\‾}{\mathrm{\ψ}}}_{m,n}(x-\frac{k}{n-m})=-{\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})---\left(15\right)$

式(14)、(15)就是谐波窗函数。可以看出，其变换简洁、直观、易实现，且由于在分解中不需要进行二抽一运算，避免了基于二进的小波分解过程频率泄漏的问题。

Claims (1)

1.一种振动信号处理的谐波窗函数，其特征在于包括如下步骤：

步骤101，构造谐波小波函数

谐波小波函数是由实函数ψ_e(x)和ψ_o(x)得到的傅立叶变换式：

角标e和o分别表示该实数是变量x的偶函数和奇函数；对式(1)和式(2)作逆傅立叶变换得到：

$\left\{\begin{array}{c}{\mathrm{\ψ}}_e\left(x\right)={\∫}_{-\∞}^{\∞}{\widehat{\mathrm{\ψ}}}_e\left(\mathrm{\ω}\right)\exp \left(\mathrm{i\ωx}\right)\mathrm{d\ω}\\ =[\sin \left(4\mathrm{\πx}\right)-\sin \left(2\mathrm{\πx}\right)]/\left(2\mathrm{\πx}\right)\\ {\mathrm{\ψ}}_o\left(x\right)={\∫}_{-\∞}^{\∞}{\widehat{\mathrm{\ψ}}}_o\left(\mathrm{\ω}\right)\exp \left(\mathrm{i\ωx}\right)\mathrm{d\ω}\\ =-[\cos \left(4\mathrm{\πx}\right)-\cos \left(2\mathrm{\πx}\right)]/\left(2\mathrm{\πx}\right)\end{array}\right.----\left(3\right)$

此定义的复函数：

ψ(x)＝[exp(i4πx)-exp(i2πx)]/(i2πx)    (4)

为谐波小波，其中i为复数，ω为频率，exp为自然对数底数；

令m＝2^j，n＝2^j+1，j∈Z+，就得到小波变换在不同分解层上的频段分解结果：

${\mathrm{\ψ}}_{m,n}\left(x\right)=\frac{e^{\mathrm{in}2\mathrm{\πx}-e^{\mathrm{im}2\mathrm{\πx}}}}{i2\mathrm{\π}(n-m)x}---\left(5\right)$

式(5)即是谐波小波在时域上的一般表达式；

给定谐波小波位移步长k/(n-m)，k∈Z+且k≠0，Z为整数，则式(5)变为

${\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})=\frac{\{e^{\mathrm{in}2\mathrm{\π}[x-k/(n-m)]}-e^{\mathrm{im}2\mathrm{\π}[x-k/(n-m)]}\}}{i2\mathrm{\π}(n-m)[x-k/(n-m)]}---\left(6\right)$

这就是带宽为(n-m)2π、分析中心在x＝k/(n-m)的谐波小波的一般表达式；对时间离散信号f^d(r)，r＝0，1，2，...，N-1，d表示离散信号，N为自然数，其谐波小波变换可写为：

$W_f(m,n,k)=\frac{(n-m)}{N}\underset{r=0}{\overset{N-1}{\mathrm{\Σ}}}{f}^d\left(r\right){\stackrel{\‾}{\mathrm{\ψ}}}_{m,n}(r-\frac{k}{n-m})---\left(7\right)$

此即信号的离散谐波小波变换表达式，其中，ψ_m，n表示ψ_m，n的共轭；其频域表达式为：

$\hat{W}(m,n,\mathrm{\ω})=\hat{f}\left(\mathrm{\ω}\right){\widehat{\stackrel{\‾}{\mathrm{\ψ}}}}_{m,n}\left[(n-m)\mathrm{\ω}\right]---\left(8\right)$

步骤102，重新定义m、n的取值范围

m，n∈R+且m＜n，即m、n在正实数域内可以取非整数值，R为实数，这样在不进行任何分解的情况下，滑动窗口到所选择的频段上并伸缩窗口：

$\left.\begin{array}{c}m=\mathrm{qB}\\ n=(q+1)B\end{array}\right\}---\left(9\right)$

其中，q∈R+，B为分析频带宽度；

步骤103，构造谐波窗函数

重新定义m、n的取值后，谐波小波函数即为谐波窗函数，谐波窗函数的三角函数表达式为：

${\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})=\left\{\cos \right[2\mathrm{\πn}(x-\frac{k}{n-m})]-\cos [2\mathrm{\πm}(x-\frac{k}{n-m})\left]\right\}/\left[i2\mathrm{\π}(n-m)(x-\frac{k}{n-m})\right]$

其共轭表达式为

${\stackrel{\‾}{\mathrm{\ψ}}}_{m,n}(x-\frac{k}{n-m})=-{\mathrm{\ψ}}_{m,n}(x-\frac{k}{n-m})---\left(11\right)$

式(10)、(11)就是谐波窗函数。

Images