CN103811018B - 一种基于Yule-Walker法的冲床噪声功率谱估计改进方法 - Google Patents

Abstract

一种基于Yule-Walker法的冲床噪声功率谱估计改进方法，步骤如下：初始化功率谱检测设备；将冲床冲裁控制信号作为开始采样触发信号；将冲床收回液压锤控制信号作为结束采样触发信号；对采集得到的噪声样本加窗处理；用迭代方法计算平均自相关函数近似值存储平均自相关函数近似值结果；重复以上步骤直到迭代计数器i=S；利用自相关函数估值求解Yule-Walker方程得到干扰白噪声的方差和p阶AR模型参数；计算冲裁噪声的功率谱估计值P_xx(ω)。

Description

一种基于Yule-Walker法的冲床噪声功率谱估计改进方法

技术领域

本发明涉及一种用于噪声控制领域的功率谱估计方法，具体是一种基于Yule-Walker法的冲床噪声功率谱估计改进方法。

背景技术

在工业生产过程中，具体的在冲床的工作过程中，冲床撞击材料会产生一种冲裁噪声。这种噪声具有：重复性、短时性、高强度性的特征。这种重复撞击噪声会导致机器设备的声疲劳，长期作用将会缩短其使用寿命，甚至发生生产事故。强烈的噪声极易形成差拍型次声波，作用与人的躯体。人体各个部位都存在固有频率，身体为7-13HZ，内脏为4-6HZ，头部为8-12HZ，这些固有频率刚好在次声波频带内，所以冲压工人在强烈噪声环境中工作，常有头昏脑胀、恶心和心悸之感。降低冲床噪声已成为噪声控制工程中的当务之急。

不论是用传统的无源消噪技术还是新型的有源消噪技术都需要对噪声进行检测，为噪声控制提供噪声的先验信息。其中最主要的信息是噪声的功率谱信息。功率谱信息能反映出噪声所含有的主要频率成分，以及各个频率成分的大小。传统的无源消噪技术对功率谱信息的依赖不是很强，部分新型的有源消噪技术需要更多的噪声信息。所以功率谱估计方法的精度直接影响这类依赖噪声先验信息的新型有源消噪技术的性能。

近几十年，已有许多学者提出了各种经典的功率谱估计方法和现代功率谱估计方法，并对其进行了深入的研究，取得了一些重要的成果。其中有一种基于参数模型估计的现代功率谱估计方法称为Yule-Walker法。这种方法首先计算信号的自相关函数，根据信号的自相关函数计算信号AR模型参数，最后由AR模型参数求得信号的功率谱估计值。但是这种方法应用到具有重复、短时、高强度噪声背景下存在一些局限：

1.由于Yule-Walker功率谱估计方法先求噪声的自相关函数，而实际中测得的噪声是离散的有限长的，计算得到的自相关函数值是估计值。在对短时噪声进行自相关函数值的计算中，自相关数据太少误差会很大。自相关函数的估计值误差会直接影响功率谱估计的误差。

2.Yule-Walker功率谱估计方法应用在短时高强度（幅值变化剧烈）噪声中严重依赖噪声检测设备的采样频率。只有噪声检测设备达到足够高的采样频率才能有效测得这类噪声的功率谱。而提高设备采样频率的成本高昂。

3.Yule-Walker功率谱估计方法不能利用噪声的重复性这一重要的先验信息。

4.Yule-Walker功率谱估计方法对高阶信号进行估计需要大量自相关函数数据计算高阶AR模型参数。

冲床作业中产生的冲裁噪声就是一类重复、短时、高强度噪声，用Yule-Walker功率谱估计方法无法有效测得冲裁噪声的功率谱信息。如何将噪声的重复性信息利用起来，在保持噪声检测设备采样率不变的前提下，提高功率谱估计性能成为冲床噪声控制工程需要解决的一个问题。

发明内容

本发明要克服现有Yule-Walker功率谱估计方法在处理重复、短时、高强度噪声时的不足，提出一种基于Yule-Walker法的冲床噪声功率谱估计改进方法。

改进方法首先利用加窗法截取多段有效噪声，然后求取每段噪声的自相关函数，再计算冲裁噪声的平均自相关函数，结合Yule-Walker法求得噪声功率谱。该方法提高了自相关函数估计精度，间接提高了功率谱估计的精度。该方法主要针对具有重复、短时、高强度噪声的功率谱估计，除了能有效应用于冲床的冲裁噪声功率谱估计，还适用于其它具有重复、短时、高强度噪声的功率谱估计，如打桩机、锻造机、射击靶场的噪声的功率谱估计。用该方法改进传统功率谱估计设备不需要改变硬件设备，只需要更新软件中的计算方法，成本低。

本发明是通过以下技术方案实现的，本发明在Yule-Walker功率谱估计方法的基础上，根据冲床的冲裁噪声的特性改进噪声自相关函数估计方法，提高检测设备对冲裁噪声的功率谱估计精度。冲床的冲裁噪声具有重复性和短时性，所以本发明用以下数学公式描述冲裁噪声信号：

x(t)=s(t)+u(t),t∈[0,∞)

T₁=ξT,(ξ≤1)

其中x(t)表示含有高斯白噪声干扰的冲裁噪声信号，s(t)表示冲裁噪声，u(t)表示高斯白噪声干扰，t表示时间，T₁表示一次冲裁噪声有效时长，T表示冲裁周期，ξ表示冲裁噪声占空比。

冲床冲裁噪声功率谱估计的具体步骤如下：

(1)初始化功率谱检测设备。设定传感器采样频率、自相关函数变量最大值、离散傅里叶变换参数、窗函数形状、窗函数长度、迭代次数计数器i初始值、总迭代次数S、AR模型阶数p。

(2)将冲床冲裁控制信号作为开始采样触发信号。等待触发信号，触发传感器开始采集冲裁噪声样本序列。

(3)将冲床收回液压锤控制信号作为结束采样触发信号。等待触发信号，触发传感器结束采集冲裁噪声x(t)的样本。

(4)对采集得到的噪声样本加窗处理。默认选择长度为N的矩形窗，也可选择改变窗的长度和形状。比较有效的窗还有汉明窗和布莱克曼窗。加窗的具体做法是对步骤(3)中得到的噪声样本进行截取或补零，样本长度大于N则截取长度为N的样本序列，样本长度小于N则在样本序列末尾补零；然后将样本序列与窗函数序列做点乘，得到一次冲裁噪声样本 $\{{\hat{x}}_{\left(i\right)}\left(n\right),n=\mathrm{0,1},...,N-1\}.$

(5)用迭代方法计算平均自相关函数近似值更新公式如下：

$\stackrel{\‾}{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(m\right)}=\left\{\begin{array}{cc}{\widehat{\mathrm{\φ}}}_1\left(m\right),& i=1\\ \frac{(i-1)\stackrel{\‾}{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(m\right)}+{\widehat{\mathrm{\φ}}}_i\left(m\right)}{i},& i>1\end{array}\right.$

${\widehat{\mathrm{\φ}}}_{\left(i\right)}\left(m\right)=\frac{1}{N}\underset{n=0}{\overset{N-\left|m\right|-1}{\mathrm{\Σ}}}{\hat{x}}_{\left(i\right)}\left(n\right){\hat{x}}_{\left(i\right)}(n+m)$

其中，m=-N+1,-N+2,...,N-1，i表示当前迭代计数器次数，表示第i次冲裁噪声的自相关函数近似值，表示第i次冲裁噪声的采样值。

(6)存储平均自相关函数近似值结果，迭代计数器加1，i=i+1。

(7)重复步骤(2)到(6)直到迭代计数器i=S。

(8)利用自相关函数估值求解Yule-Walker方程得到干扰白噪声的方差和p阶AR模型参数。Yule-Walker方程如下：

求解Yule-Walker方程的具体步骤如下：

(i)计算一阶模型AR(1)的参数a₁₁和对应方差解如下方程：

$\left[\begin{array}{cc}{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)& {\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(1\right)\\ {\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(1\right)& {\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_{11}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\σ}}^2\\ 0\end{array}\right]$

解为：

$a_{11}=\frac{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(1\right)}{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)},{\mathrm{\σ}}^2=(1-{\left|a_{11}\right|}^2){\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)$

(ii)设定迭代计数器k=2，用迭代法计算p阶模型AR(p)的参数{a_p,1,a_p,2,...,a_p,p}和对应的方差迭代更新公式如下：

$a_{k,k}=-\frac{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(k\right)+\underset{l=1}{\overset{k-1}{\mathrm{\Σ}}}{a}_{k-1,l}{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}(k-1)}{{\mathrm{\σ}}_{k-1}^2}$

a_k,j=a_k-1,j+a_k,ka_k-1,k-1,j=1,2,...,k-1

${\mathrm{\σ}}_k^2=(1-{{|a}_{\mathrm{kk}}|}^2){\mathrm{\σ}}_{k-1}^2$

(iii)迭代计数器k加1，k=k+1，重复步骤(ii)，直到k=p。

(9)计算冲裁噪声的功率谱估计值P_xx(ω)，公式如下：

$P_{\mathrm{xx}}\left(\mathrm{\ω}\right)=\frac{{\mathrm{\σ}}_p^2}{{|1+\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{a}_k{e}^{-\mathrm{j\ωk}}|}^2}$

其中，ω表示频率；

在冲床噪声检测控制工程中采用本发明提出的方法能够获得足够的功率谱估计精度和分辨率，能抑制白噪声的干扰。

本发明的优点是：用加窗法和平均法计算冲裁噪声自相关函数近似值，解决了传统方法因自相关函数近似值方差大导致的功率谱估计精度低的缺陷，且方法简单，易于实现。

附图说明

图1为采用本发明方法的程序流程图。

图2为本发明实施例中冲裁噪声十个周期内的检测图。

图3为本发明实施例中第一次冲裁噪声的检测图。

图4为本发明实施例中未改进方法和改进方法得到的功率谱比较图。

具体实施方式

以下结合附图和实施例对本发明的技术方案作进一步描述。

如图1所示，功率谱检测设备首先初始化设备参数，然后等待冲床冲裁控制信号，触发传感器采集噪声样本序列，当冲床收回液压锤控制信号发出时，触发传感器停止采集噪声样本，这样就完成了冲裁噪声截取和加矩形窗的步骤。加窗程序不做任何处理就是选择矩形窗，也可根据噪声的特点选择非矩形窗，进一步提高算法性能。这里功率谱检测设备使用了迭代的方法计算平均自相关函数近似值，不需要存储多次冲裁噪声，节省内存空间。然后求解Yule-Walker方程得到干扰白噪声的方差和p阶AR模型参数，以此求得冲裁噪声的功率谱估计值。

如图2所示，冲裁噪声周期为1秒，一次冲裁噪声持续时间为0.1秒，信噪比为10dB。以此为实施例，本发明的冲裁噪声功率谱估计流程如下：

(1)初始化功率谱检测设备。设定传感器采样频率为40KHz，自相关函数变量最大值为3999，离散傅里叶变换数据点数为4096，窗函数选择长度为4000的矩形窗，迭代次数计数器i为1，总迭代次数为10，AR模型阶数为400。

(2)将冲床冲裁控制信号作为开始采样触发信号。等待触发信号，触发传感器开始采集冲裁噪声样本序列。

(3)将冲床收回液压锤控制信号作为结束采样触发信号。等待触发信号，触发传感器结束采集冲裁噪声样本。

(4)对采集得到的噪声样本加窗处理。默认选择长度为4000的矩形窗，也可选择改变窗的长度和形状。比较有效的窗还有汉明窗和布莱克曼窗。加窗的具体做法是对步骤(3)中得到的噪声样本进行截取或补零，样本长度大于4000则截取长度为4000的样本序列，样本长度小于4000则在样本序列末尾补零。然后将样本序列与窗函数序列做点乘，得到一次冲裁噪声样本如图3所示为第一次冲裁噪声的检测图。

(5)用迭代方法计算平均自相关函数近似值更新公式如下：

$\stackrel{\‾}{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(m\right)}=\left\{\begin{array}{cc}{\widehat{\mathrm{\φ}}}_1\left(m\right),& i=1\\ \frac{(i-1)\stackrel{\‾}{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(m\right)}+{\widehat{\mathrm{\φ}}}_i\left(m\right)}{i},& i>1\end{array}\right.$

${\widehat{\mathrm{\φ}}}_{\left(i\right)}\left(m\right)=\frac{1}{4000}\underset{n=0}{\overset{4000-\left|m\right|-1}{\mathrm{\Σ}}}{\hat{x}}_{\left(i\right)}\left(n\right){\hat{x}}_{\left(i\right)}(n+m)$

其中，m=-3999,-3998,...,3999，i表示当前迭代计数器次数，表示第i次冲裁噪声的自相关函数近似值，表示第i次冲裁噪声的采样值。

(6)存储平均自相关函数近似值结果，迭代计数器加1，i=i+1。

(7)重复步骤(2)到(6)直到迭代计数器i=10。

(8)利用自相关函数估值求解Yule-Walker方程得到干扰白噪声的方差和400阶AR模型参数。

(i)计算一阶模型AR(1)的参数a₁₁和对应方差解如下方程：

$\left[\begin{array}{cc}{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)& {\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(1\right)\\ {\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(1\right)& {\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_{11}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\σ}}^2\\ 0\end{array}\right]$

解为：

$a_{11}=\frac{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(1\right)}{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)},{\mathrm{\σ}}^2=(1-{\left|a_{11}\right|}^2){\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(0\right)$

(ii)设定迭代计数器k=2，用迭代法计算400阶模型AR(400)的参数{a_400,1,a_400,2,...,a_400,400}和对应的方差迭代更新公式如下：

$a_{k,k}=-\frac{{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}\left(k\right)+\underset{l=1}{\overset{k-1}{\mathrm{\Σ}}}{a}_{k-1,l}{\widehat{\mathrm{\φ}}}_{\mathrm{xx}}(k-1)}{{\mathrm{\σ}}_{k-1}^2}$

a_k,j=a_k-1,j+a_k,ka_k-1,k-1,i=1,2,...,k-1

${\mathrm{\σ}}_k^2=(1-{{|a}_{\mathrm{kk}}|}^2){\mathrm{\σ}}_{k-1}^2$

(iii)迭代计数器k加1，k=k+1，重复步骤(ii)，直到k=400。

(9)计算冲裁噪声的功率谱估计值P_xx(ω)，公式如下：

$P_{\mathrm{xx}}\left(\mathrm{\ω}\right)=\frac{{\mathrm{\σ}}_{400}^2}{{|1+\underset{k=1}{\overset{400}{\mathrm{\Σ}}}{a}_k{e}^{-\mathrm{j\ωk}}|}^2}$

其中，ω表示频率；

结果显示在图4中，其中实线是未改进算法对冲裁噪声的功率谱估计结果，虚线是本发明的改进算法对冲裁噪声的功率谱估计结果。虚线的波峰明显，通过虚线可以测出噪声中含有9个主要频率分量：130Hz,290Hz,400Hz,500Hz,611Hz,772Hz,810Hz,881Hz,1000Hz。

Claims (1)

1.一种基于Yule-Walker法的冲床噪声功率谱估计改进方法，具体步骤如下：

(1)初始化功率谱检测设备；设定传感器采样频率、自相关函数变量最大值、离散傅里叶变换参数、窗函数形状、窗函数长度、迭代次数计数器i初始值、总迭代次数S、AR模型阶数p；

(2)将冲床冲裁控制信号作为开始采样触发信号；等待触发信号，触发传感器开始采集冲裁噪声样本序列；

(3)将冲床收回液压锤控制信号作为结束采样触发信号；等待触发信号，触发传感器结束采集冲裁噪声x(t)的样本；

(4)对采集得到的噪声样本加窗处理；选择长度为N的矩形窗、汉明窗或布莱克曼窗；加窗的具体做法是对步骤(3)中得到的噪声样本进行截取或补零，样本长度大于N则截取长度为N的样本序列，样本长度小于N则在样本序列末尾补零；然后将样本序列与窗函数序列做点乘，得到一次冲裁噪声样本 $\{{\hat{x}}_{\left(i\right)}\left(n\right),n=0,1,...,N-1\};$

(5)用迭代方法计算平均自相关函数近似值更新公式如下：

$\stackrel{\‾}{{\widehat{\mathrm{\φ}}}_{xx}\left(m\right)}=\left\{\begin{array}{cc}{\widehat{\mathrm{\φ}}}_1\left(m\right),& i=1\\ \frac{\left(i-1\right)\stackrel{\‾}{{\widehat{\mathrm{\φ}}}_{xx}\left(m\right)}+{\widehat{\mathrm{\φ}}}_i\left(m\right)}{i},& i>1\end{array}\right.$

${\widehat{\mathrm{\φ}}}_{\left(i\right)}\left(m\right)=\frac{1}{N}\underset{n=0}{\overset{N-\left|m\right|-1}{\Σ}}{\hat{x}}_{\left(i\right)}\left(n\right){\hat{x}}_{\left(i\right)}(n+m)$

其中，m＝-N+1,-N+2,...,N-1，i表示当前迭代计数器次数，表示第i次冲裁噪声的自相关函数近似值，表示第i次冲裁噪声的采样值；

(6)存储平均自相关函数近似值结果，迭代计数器加1，i＝i+1；

(7)重复步骤(2)到(6)直到迭代计数器i＝S；

(8)利用自相关函数估值求解Yule-Walker方程得到干扰白噪声的方差和p阶AR模型参数，Yule-Walker方程如下：

求解Yule-Walker方程的具体步骤如下：

(i)计算一阶模型AR(1)的参数a₁₁和对应方差解如下方程：

$\left[\begin{array}{cc}{\widehat{\mathrm{\φ}}}_{xx}\left(0\right)& {\widehat{\mathrm{\φ}}}_{xx}\left(1\right)\\ {\widehat{\mathrm{\φ}}}_{xx}\left(1\right)& {\widehat{\mathrm{\φ}}}_{xx}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_{11}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\σ}}^2\\ 0\end{array}\right]$

解为：

$a_{11}=\frac{{\widehat{\mathrm{\φ}}}_{xx}\left(1\right)}{{\widehat{\mathrm{\φ}}}_{xx}\left(0\right)},{\mathrm{\σ}}^2=(1-|a_{11}{|}^2){\widehat{\mathrm{\φ}}}_{xx}\left(0\right)$

(ii)设定迭代计数器k＝2，用迭代法计算p阶模型AR(p)的参数{a_p,1,a_p,2,...,a_p,p}和对应的方差迭代更新公式如下：

$a_{k,k}=-\frac{{\widehat{\mathrm{\φ}}}_{xx}\left(k\right)+\underset{l=1}{\overset{k-1}{\Σ}}{a}_{k-1,l}{\widehat{\mathrm{\φ}}}_{xx}(k-1)}{{\mathrm{\σ}}_{k-1}^2}$

a_k,j＝a_k-1,j+a_k,ka_k-1,k-1,j＝1,2,...,k-1

${\mathrm{\σ}}_k^2=(1-|a_{kk}{|}^2){\mathrm{\σ}}_{k-1}^2$

(iii)迭代计数器k加1，k＝k+1，重复步骤(ii)，直到k＝p；

(9)计算冲裁噪声的功率谱估计值P_xx(ω)，公式如下：

$P_{xx}\left(\mathrm{\ω}\right)=\frac{{\mathrm{\σ}}_p^2}{|1+\underset{k=1}{\overset{p}{\Σ}}{a}_k{e}^{-j\mathrm{\ω}k}{|}^2}$

其中，ω表示频率。