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

Abstract

基于Burg法的冲床噪声功率谱估计改进方法，具体步骤如下：初始化功率谱检测设备；将冲床冲裁控制信号作为开始采样触发信号；将冲床收回液压锤控制信号作为结束采样触发信号；对采集得到的噪声样本加窗处理；用迭代方法计算有效噪声的平均值序列；存储有效噪声平均值序列；重复以上步骤直到迭代计数器i=S；计算反射系数K_p、干扰噪声方差和AR模型参数{a_p,1,a_p,2,...,a_p,p}；计算冲裁噪声的功率谱估计值P_xx(ω)，公式如下：<maths num="0001"></maths>。

Description

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

技术领域

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

背景技术

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

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

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

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

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

3.Burg功率谱估计方法受Levinson迭代关系式的约束，功率谱估计结果存在谱线分裂和频率偏移现象。

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

发明内容

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

改进方法首先利用加窗法截取多段有效噪声，然后求取有效噪声的平均值序列，再从均值序列求得反射系数，利用Levinson递推算法与反射系数来求得AR参数，最后根据AR参数求得噪声功率谱。该方法间接优化了Burg算法中的反射系数，提高了功率谱估计的精度和分辨率。该方法主要针对具有重复、短时、高强度噪声的功率谱估计，除了能有效应用于冲床的冲裁噪声功率谱估计，还适用于其它具有重复、短时、高强度噪声的功率谱估计，如打桩机、锻造机、射击靶场的噪声的功率谱估计。用该方法改进传统功率谱估计设备不需要改变硬件设备，只需要更新软件中的计算方法，成本低。

本发明是通过以下技术方案实现的，本发明在Burg功率谱估计方法的基础上，根据冲床的冲裁噪声的特性优化Burg算法中的反射系数，间接提高检测设备对冲裁噪声的功率谱估计精度。冲床的冲裁噪声具有重复性和短时性，所以本发明用以下数学公式描述冲裁噪声信号：

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

T₁=ξT,(ξ≤1)

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

本发明所述的基于Burg法的冲床噪声功率谱估计改进方法，具体步骤如下：

(1)初始化功率谱检测设备；设定传感器采样频率、窗函数类型、窗函数长度、迭代次数计数器i初始值、总迭代次数S、AR模型阶数p；

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

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

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

(5)用迭代方法计算有效噪声的平均值序列；迭代更新公式如下：

$\stackrel{\&OverBar;}{\hat{x}\left(n\right)}=\left\{\begin{array}{cc}{\hat{x}}_{\left(1\right)}\left(n\right)& i=1\\ \frac{(i-1)\stackrel{\&OverBar;}{\hat{x}\left(n\right)}+{\hat{x}}_{\left(i\right)}\left(n\right)}{i}& i>1\end{array}\right.$

其中，i表示当前迭代计数器次数；

(6)存储有效噪声平均值序列，迭代计数器i加1，i=i+1；

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

(8)计算反射系数K_p、干扰噪声方差和AR模型参数{a_p,1,a_p,2,...,a_p,p}；

(i)初始化前向误差e₀(n)，后向误差b₀(n)，干扰方差估计值迭代次数计数器k=1，具体公式如下：

$e_0\left(n\right)=\stackrel{\&OverBar;}{\hat{x}\left(n\right)}$

$b_0\left(n\right)=\stackrel{\&OverBar;}{\hat{x}\left(n\right)}$

${\mathrm{\&sigma;}}_0^2=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{x}^2\left(n\right)$

(ii)计算K_k，计算公式如下：

$K_k=-\frac{2\underset{n=k}{\overset{N-1}{\mathrm{\&Sigma;}}}\left[e_{k-1}\left(n\right)b_{k-1}(n-1)\right]}{\underset{n=k}{\overset{N-1}{\mathrm{\&Sigma;}}}[e_{k-1}^2\left(n\right)+b_{k-1}^2(n-1)]}$

(iii)计算k阶AR模型参数a_k,i(i=1,2,...,k-1)，公式如下：

a_k,k=K_k

a_k,i=a_k-1,i+K_ka_k-1,k-i,(i=1,2,...,k-1)

(iv)更新前向误差e_k(n)和后向误差b_k(n)，干扰方差估计值迭代更新公式如下：

e_k(n)=e_k-1(n)+K_kb_k-1(n-1)

b_k(n)=b_k-1(n-1)+K_ke_k-1(n)

${\mathrm{\&sigma;}}_k^2=(1-K_k^2){\mathrm{\&sigma;}}_{k-1}^2$

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

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

$P_{\mathrm{xx}}\left(\mathrm{\&omega;}\right)=\frac{{\mathrm{\&sigma;}}_p^2}{{|1+\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_k{e}^{-\mathrm{j\&omega;k}}|}^2}$

其中，ω表示频率；

在冲床噪声检测控制工程中采用本发明提出的方法能够获得足够的功率谱估计精度和分辨率，能抑制白噪声的干扰。本发明最大的特点就是：用加窗法和平均法计算冲裁噪声平均值，间接优化Burg算法中的反射系数，解决了传统方法的功率谱估计存在谱线分裂和频率偏移现象的缺陷，且方法简单，易于实现。

附图说明

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

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

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

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

具体实施方式

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

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

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

(1)初始化功率谱检测设备。设定传感器采样频率为40KHz，窗函数选择长度为4000的矩形窗，迭代次数计数器i为1，总迭代次数为10，AR模型阶数为400。

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

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

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

(5)用迭代方法计算有效噪声的平均值序列。迭代更新公式如下：

$\stackrel{\&OverBar;}{\hat{x}\left(n\right)}=\left\{\begin{array}{cc}{\hat{x}}_{\left(1\right)}\left(n\right)& i=1\\ \frac{(i-1)\stackrel{\&OverBar;}{\hat{x}\left(n\right)}+{\hat{x}}_{\left(i\right)}\left(n\right)}{i}& i>1\end{array}\right.$

其中，i表示当前迭代计数器次数。

(6)存储有效噪声平均值序列，迭代计数器i加1，i=i+1。

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

(8)计算反射系数K₄₀₀、干扰噪声方差和AR模型参数{a_400,1,a_400,2,...,a_400,400}。

(i)初始化前向误差e₀(n)，后向误差b₀(n)，干扰方差估计值迭代次数计数器k=1，具体公式如下：

$e_0\left(n\right)=\stackrel{\&OverBar;}{\hat{x}\left(n\right)}$

$b_0\left(n\right)=\stackrel{\&OverBar;}{\hat{x}\left(n\right)}$

${\mathrm{\&sigma;}}_0^2=\frac{1}{4000}\underset{n=0}{\overset{3999}{\mathrm{\&Sigma;}}}{x}^2\left(n\right)$

(ii)计算K_k，计算公式如下：

$K_k=-\frac{2\underset{n=k}{\overset{3999}{\mathrm{\&Sigma;}}}\left[e_{k-1}\left(n\right)b_{k-1}(n-1)\right]}{\underset{n=k}{\overset{3999}{\mathrm{\&Sigma;}}}[e_{k-1}^2\left(n\right)+b_{k-1}^2(n-1)]}$

(iii)计算k阶AR模型参数a_k,i(i=1,2,...,k-1)，公式如下：

a_k,k=K_k

a_k,i=a_k-1,i+K_ka_k-1,k-i,(i=1,2,...,k-1)

(iv)更新前向误差e_k(n)和后向误差b_k(n)，干扰方差估计值迭代更新公式如下：

e_k(n)=e_k-1(n)+K_kb_k-1(n-1)

b_k(n)=b_k-1(n-1)+K_ke_k-1(n)

${\mathrm{\&sigma;}}_k^2=(1-K_k^2){\mathrm{\&sigma;}}_{k-1}^2$

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

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

$P_{\mathrm{xx}}\left(\mathrm{\&omega;}\right)=\frac{{\mathrm{\&sigma;}}_{400}^2}{{|1+\underset{k=1}{\overset{400}{\mathrm{\&Sigma;}}}{a}_k{e}^{-\mathrm{j\&omega;k}}|}^2}$

其中，ω表示频率。

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

Claims (1)

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

(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{\&OverBar;}{\hat{x}\left(n\right)}=\left\{\begin{array}{cc}{\hat{x}}_{\left(1\right)}\left(n\right)& i=1\\ \frac{(i-1)\stackrel{\&OverBar;}{\hat{x}\left(n\right)}+{\hat{x}}_{\left(i\right)}\left(n\right)}{i}& i>1\end{array}\right.$

其中，i表示当前迭代计数器次数；

(6)存储有效噪声平均值序列，迭代计数器i加1，i＝i+1；

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

(8)计算反射系数K_p、干扰噪声方差和AR模型参数{a_p,1,a_p,2,...,a_p,p}；

(i)初始化前向误差e₀(n)，后向误差b₀(n)，干扰方差估计值迭代次数计数器k＝1，具体公式如下：

$e_0\left(n\right)=\stackrel{\&OverBar;}{\hat{x}\left(n\right)}$

$b_0\left(n\right)=\stackrel{\&OverBar;}{\hat{x}\left(n\right)}$

${\mathrm{\&sigma;}}_0^2=\frac{1}{N}\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}^2\left(n\right)$

(ii)计算K_k，计算公式如下：

$K_k=-\frac{2\underset{n=k}{\overset{N-1}{\&Sigma;}}\&lsqb;e_{k-1}\left(n\right)b_{k-1}(n-1)\&rsqb;}{\underset{n=k}{\overset{N-1}{\&Sigma;}}\&lsqb;e_{k-1}^2\left(n\right)+b_{k-1}^2(n-1)\&rsqb;}$

(iii)计算k阶AR模型参数a_k,i(i＝1,2,...,k-1)，公式如下：

a_k,k＝K_k

a_k,i＝a_k-1,i+K_ka_k-1,k-i,(i＝1,2,...,k-1)

(iv)更新前向误差e_k(n)和后向误差b_k(n)，干扰方差估计值迭代更新公式如下：

e_k(n)＝e_k-1(n)+K_kb_k-1(n-1)

b_k(n)＝b_k-1(n-1)+K_ke_k-1(n)

${\mathrm{\&sigma;}}_k^2=(1-K_k^2){\mathrm{\&sigma;}}_{k-1}^2$

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

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

$P_{xx}\left(\mathrm{\&omega;}\right)=\frac{{\mathrm{\&sigma;}}_p^2}{|1+\underset{k=1}{\overset{p}{\&Sigma;}}{a}_k{e}^{-j\mathrm{\&omega;}k}{|}^2}$

其中，ω表示频率。