精确识别振动信号主频的方法及系统

1.本发明涉及数字信号处理领域，尤其涉及一种对振动信号主频精确识别的方法及系统。


背景技术：

2.信号频率识别作为数字信号处理研究领域中一项重要工作，具有显著的工程应用价值，它在非介入式内燃机转速测量、机械运行故障特征提取、电网谐波监测等领域都存在重要应用。信号频率识别根据是否对被测对象建立数学模型可分为参数化法与非参数化法，其中，非参数化法无需建立特定的数学模型，具有良好的适应性，但是此方法进行频率识别存在识别精度差、计算量大的问题；参数化法根据被测对象建立数学模型进行匹配，具有更高的识别精度和较小的计算量，但此方法存在模型建立难、参数适配、初值设定等问题。因此，从频率识别的普适化角度而言，利用非参数化法对信号主频进行快速精确的识别，具有更重要的研究意义。
3.目前，在非参数化法频率识别方面，国内外学者做了一定的研究工作。极大似然估计法和非线性最小二乘法
1.受限于计算量大且需要较高的采样频率，难以适应快速实时性的要求；fu h
2.基于时域进行信号分析，提出根据自相关函数的几何意义，简化频率估计过程，从而减少频率估计的计算量提高单个正弦波的频率估计的速度与精度；高志峰等人
3.提出根据迭代策略进行谱峰频率搜索的方法，能够直接估计信号频率，但此算法在极值点附近收敛速度较慢，并没有显著提高频率估计的精度。
4.【参考文献】
5.[1]stoica p,nehoral a.statistical analysis of two non
‑
linear least squares estimators of sine waves parameters in the colored noise[j].proceddings of the icassp,1998,4:2408
‑
2411.
[0006]
[2]fu h,kam p.sample autocorrelation function based frequency estimation of a single sinusoid in awgn[c].vehicular technology conference,ieee 75th,2012:1
‑
5.
[0007]
[3]高志峰,彭喜元,彭宇.基于迭代更新策略的快速高精度频率估计方法[j].振动与冲击,2015,34(14):16
‑
20.


技术实现要素：

[0008]
本技术提供一种精确识别振动信号主频的方法，解决由于频率泄漏导致的信号主频识别精度较低，速度较慢，误差较大的问题。
[0009]
根据本发明实施例的一方面，提供一种识别振动信号主频的方法，包括：
[0010]
步骤1：对振动信号以频率f
s
进行采样，得到被测信号的ad转换数据序列x(t)；
[0011]
步骤2：使用随机数种子生成一个随机起始点k
start
，在数据序列x(t)中截取t＝[1:k
start
]，长度为k的序列x1(t)；使用线性插补压缩算法将长度为k的序列x1(t)压缩为长度为
n的序列g(t)；
[0012]
步骤3：对压缩后的序列g(t)进行傅里叶变换得到g(t)的频谱函数l(w)，均零化处理并计算幅值最大处的频率值w和幅值l(w)，再分别计算幅值最大处的频率值w左侧频率点的幅值l(w
‑
1)与最大幅值l(w)的比值l(w
‑
1)/l(w)，幅值最大处的频率值w右侧频率点的幅值l(w+1)与最大幅值l(w)的比值l(w+1)/l(w)，取两个比值中较小的值作为ratio(k)的函数值；
[0013]
步骤4：以起始点k
start
为中心分别沿数轴向左取点k
left
，沿数轴向右取点k
right
，以k
start
作为父节点，k
left
和k
right
作为两个子节点，重复步骤2所述序列截断和线性插补压缩算法以及步骤3，得出ratio(k
left
)，ratio(k
right
)，比较父节点和两个子节点的ratio(k)，若最小ratio(k)存在于子节点则以此子节点为下一分支的父节点重复步骤4，直到父节点的ratio(k)小于两个子节点的ratio(k)时，历遍两子节点所限定的k值范围，此时最小的ratio(k)所对应的自变量k0为被测信号能够被正周期截取的点数，k0处的n点离散傅里叶变换频谱图中主频信号的数字角频率为w0，幅值为l(w0)，被测信号的频率为
[0014][0015]
根据本发明实施例的第二方面，提供一种识别振动信号主频的系统，包括：处理器；用于存储处理器可执行指令的存储器；其中，所述处理器被配置为执行所述方法的全部或部分步骤。
[0016]
根据本发明实施例的第三方面，提供一种非临时性计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现项所述方法的全部或部分步骤。
附图说明
[0017]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍。
[0018]
图1是根据本发明一实施例提供的一种精确识别振动信号主频的方法流程图。
[0019]
图2是采样得到原始信号x(t)后截取k＝4500点后的信号x1(t)的函数图像。
[0020]
图3是分别截取长度为k＝4100，k＝4300，k＝4500，k＝4700，k＝4900时的x1(t)序列后进行压缩得到的固定长度信号g(t)的函数图像。
[0021]
图4是对k＝4500时信号g(t)的4096点快速傅里叶变换频谱图。
[0022]
图5是本发明实施例提供的一种精确识别振动信号主频的方法中分支定界法剪枝搜索的过程示意图。
[0023]
图6是使用分支定界法快速搜索ratio(k)的函数图像。
[0024]
图7是对原始信号分别截取点数为k＝4455、k＝4500、k＝4620、k＝4680、k＝4720、k＝4765的信号再经过压缩后进行fft的频谱图。
具体实施方式
[0025]
在基于频域的信号主频估计算法中，由于被测信号的采样频率限制以及时域离散化和频域离散化进行傅里叶变换，特别是当信号采样长度为主频信号周期的非整数倍时，信号频谱在主频附近产生频率泄漏和频率拖尾效应，导致信号频率估计不准确。因此对主
频的估计，一方面可以通过截取合适的信号长度进行傅里叶变换，最大程度抑制频率泄露的影响(如图7所示，不同截取长度，频率泄露程度不同)，提高信号频率识别精度；另一方面运用插补运算等手段提高计算速度，为快速准确的主频估计提供一种有效的方法。
[0026]
图1是根据本发明一实施例提供的一种精确识别振动信号主频的方法流程图。如图1所示，该方法在测量时，为减少迭代计算的次数，提高计算速度，采用分支定界法对被测信号的主频所在区间进行分支预测，迭代推算最优解所在区间；对被测信号截取不同长度的傅里叶时间窗，为实现n点基4快速傅里叶变换，在时域将截断后的信号采用线性插补压缩算法压缩至n点(n＝4^
n
)；以最小的旁瓣主瓣比(旁主比)为评价该点的频率识别准确度的参考依据，旁主比越小则频率泄漏越小，识别越精确，进而可以准确计算出被测信号的主频和频谱幅值；该方法能够以较少的迭代运算次数，精确的寻找到最佳的离散傅里叶时间窗截取长度，有效的抑制信号的频率泄漏和“栅栏效应”，实现信号主频的快速高精度识别。下文通过步骤1～步骤4对本发明进行详细说明。
[0027]
步骤1：在满足采样定律的条件下，对振动信号以频率fs进行采样，得到被测信号的ad转换数据序列x(t)。
[0028]
本实施例中，假设被测量的信号为x(t)＝100sin(2πf0t)+30ξ，仿真中，将信号的频率设为f0＝330hz，噪声ξ为(0～1)之间的随机数。以采样频率fs＝5000hz对待测信号进行等间隔采样，获得5000个原始采样数据组成的时间序列x(t)。
[0029]
步骤2：在k可取的范围内使用随机数种子生成一个随机起始点k
start
，在数据序列x(t)中截取t＝[1:k
start
]，长度为k的序列x1(t)；使用线性插补压缩算法将长度为k的序列x1(t)压缩为长度为n的序列g(t)。
[0030]
本实施例中，在k可取的范围内使用随机数种子生成一个随机起始点k
start
＝4500，取初始傅里叶时间窗截取长度k＝4500截取x(t)得到序列x1(t)，如图2所示；将长度为k＝4500的序列x1(t)压缩至4096点的序列g(t)，其中g(t)与x1(t)的压缩点的关系满足：
[0031][0032]
式中，k是信号x(t)被截取的长度，n是信号压缩长度，t是时间序列的自变量，t∈[1:n]，n是原始信号采样的点数，比如采样频率fs是5000hz，则n也是5000，[]表示取整，{}表示取小数；经过此方法压缩后的序列g(t)与x1(t)相比具有相同的周期数与幅值，如图3所示。
[0033]
步骤3：对步骤2压缩后的序列g(t)进行基4fft得到g(t)的频谱函数l(w)，均零化处理并计算幅值最大处的频率值w和幅值l(w)，再分别计算幅值最大处的频率值w左侧频率点的幅值l(w
‑
1)与最大幅值l(w)的比值l(w
‑
1)/l(w)，幅值最大处的频率值w右侧频率点的幅值l(w+1)与最大幅值l(w)的比值l(w+1)/l(w)，取两个比值中较小的值作为ratio(k)的函数值。
[0034]
本实施例中，对g(t)进行4096点基4fft得到g(t)的频谱函数l(w)，均零化处理并计算幅值最大处的频率值w＝15和幅值l(15)＝161400，如图4所示，再分别计算幅值最大处的频率值w左侧频率点的幅值l(14)＝95160与最大幅值l(15)＝161400的比值l(14)/l(15)
＝0.590，幅值最大处的频率值w＝15右侧频率点的幅值l(16)＝42280与最大幅值l(15)＝161400的比值l(16)/l(15)＝0.262，取两个比值中较小的值0.262，ratio(k＝4500)＝0.262。
[0035]
步骤4：以起始点k
start
为中心分别沿数轴向左取点k
left
，沿数轴向右取点k
right
，以k
start
作为父节点，k
left
和k
right
作为两个子节点，重复步骤2所述序列截断和线性插补压缩算法以及步骤3，得出ratio(k
left
)，ratio(k
right
)，比较父节点和两个子节点的ratio(k)，若最小ratio(k)存在于子节点则以此子节点为下一分支的父节点重复步骤4，直到父节点的ratio(k)小于两个子节点的ratio(k)时，历遍两子节点所限定的k值范围，此时最小的ratio(k)所对应的自变量k0为被测信号能够被正周期截取的点数，k0处的n点离散傅里叶变换频谱图中主频信号的数字角频率为w0，幅值为l(w0)，被测信号的频率为
[0036]
本实施例中，采用分支定界法以起始点k
start
＝4500为中心分别沿数轴向左取点k
left
＝4490，沿数轴向右取点k
right
＝4510，以k
start
做为父节点，k
left
和k
right
做为两个子节点，重复步骤2所述序列截断和线性插补压缩算法以及步骤3，得出ratio(k
left
)＝0.2786，ratio(k
right
)＝0.2443，ratio(k
start
)＝0.262，比较父节点和两个子节点的ratio(k)。
[0037]
此时，最小下界位于ratio(k
right
)中，可以预测，最优解k位于k≥k
right
＝4510的范围内，以k
right
做为新的分支的父节点，重复步骤4，如图5所示，设定左节点和右节点的k值计算ratio(k)，直到父节点的ratio(k)小于两个子节点的ratio(k)时，此时，ratio(k
left
＝4610)＝0.0315，ratio(k
right
＝4630)＝0.0282，ratio(k
start
＝4620)＝0.000264，可以预测，最优解位于区间k∈[4610:4630]，历遍此最优解所在区间，如图6所示，此时最小的ratio(k)＝0.000264所对应的自变量k0＝4620为被测信号能够被正周期截取的点数，k0处的n＝4096点离散傅里叶变换频谱图中主频信号的数字角频率为w0＝14，幅值l(w0)＝264900，被测信号的频率f＝(14x5000)/4620＝15.15hz。
[0038]
在示例性实施例中，还提供一种系统，该系统包括处理器其中，所述处理器被配置为执行所述方法的全部或部分步骤。
[0039]
在示例性实施例中，还提供了一种非临时性计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现所述方法的全部或部分步骤。例如，所述非临时性计算机可读存储介质可以是rom、ram、cd
‑
rom、磁带、软盘和光数据存储设备等。