CN113392813B

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

Info

Publication number: CN113392813B
Application number: CN202110872668.5A
Authority: CN
Inventors: 王选择; 王爱辉; 尹晋平; 张瑜灿; 王恒辉; 翟中生
Original assignee: Hubei University of Technology
Current assignee: Hubei University of Technology
Priority date: 2021-07-30
Filing date: 2021-07-30
Publication date: 2022-04-26
Anticipated expiration: 2041-07-30
Also published as: CN113392813A

Abstract

公开了一种精确识别振动信号主频的方法及方法，该方法包括：对振动信号以频率f_s进行采样，得到被测信号的AD转换数据序列；设定傅里叶时间窗长度k截取数据序列，采用线性插补压缩算法将序列长度压缩至固定长度N(N＝4ⁿ)，以进行傅里叶变换；根据傅里叶变换结果，以频谱图中幅值最大点的旁主比作为本次计算的结果ratio(k)，ratio(k)越小，说明对应主频点计算越准确；采用分支定界法快速寻找函数ratio(k)的最优解k₀即ratio(k)的最小值，此时k₀最接近主频周期采样点数的整倍数，应用k₀与主频位置信息计算出被测信号的主频和幅值。本发明可以有效的解决振动信号频率识别中频率泄漏和频率拖尾效应导致信号频率估计不准确的问题，精确识别出被测信号主频及信号频谱的幅值。

Description

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

技术领域

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

背景技术

信号频率识别作为数字信号处理研究领域中一项重要工作，具有显著的工程应用价值，它在非介入式内燃机转速测量、机械运行故障特征提取、电网谐波监测等领域都存在重要应用。信号频率识别根据是否对被测对象建立数学模型可分为参数化法与非参数化法，其中，非参数化法无需建立特定的数学模型，具有良好的适应性，但是此方法进行频率识别存在识别精度差、计算量大的问题；参数化法根据被测对象建立数学模型进行匹配，具有更高的识别精度和较小的计算量，但此方法存在模型建立难、参数适配、初值设定等问题。因此，从频率识别的普适化角度而言，利用非参数化法对信号主频进行快速精确的识别，具有更重要的研究意义。

目前，在非参数化法频率识别方面，国内外学者做了一定的研究工作。极大似然估计法和非线性最小二乘法^[1]受限于计算量大且需要较高的采样频率，难以适应快速实时性的要求；Fu H^[2]基于时域进行信号分析，提出根据自相关函数的几何意义，简化频率估计过程，从而减少频率估计的计算量提高单个正弦波的频率估计的速度与精度；高志峰等人^[3]提出根据迭代策略进行谱峰频率搜索的方法，能够直接估计信号频率，但此算法在极值点附近收敛速度较慢，并没有显著提高频率估计的精度。

【参考文献】

[1]Stoica P,Nehoral A.Statistical analysis of two non-linear leastsquares estimators of sine waves parameters in the colored noise[J].Proceddings of the ICASSP,1998,4:2408-2411.

[2]Fu H,Kam P.Sample autocorrelation function based frequencyestimation of a single sinusoid in AWGN[C].Vehicular Technology Conference,IEEE 75th,2012:1-5.

[3]高志峰,彭喜元,彭宇.基于迭代更新策略的快速高精度频率估计方法[J].振动与冲击,2015,34(14):16-20.

发明内容

本申请提供一种精确识别振动信号主频的方法，解决由于频率泄漏导致的信号主频识别精度较低，速度较慢，误差较大的问题。

根据本发明实施例的一方面，提供一种识别振动信号主频的方法，包括：

步骤1：对振动信号以频率f_s进行采样，得到被测信号的AD转换数据序列X(t)；

步骤2：使用随机数种子生成一个随机起始点k_start，在数据序列X(t)中截取t＝[1:k_start]，长度为k的序列X₁(t)；使用线性插补压缩算法将长度为k的序列X₁(t)压缩为长度为N的序列G(t)；

步骤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)的函数值；

步骤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)所对应的自变量k₀为被测信号能够被正周期截取的点数，k₀处的N点离散傅里叶变换频谱图中主频信号的数字角频率为w₀，幅值为L(w₀)，被测信号的频率为

根据本发明实施例的第二方面，提供一种识别振动信号主频的系统，包括：处理器；用于存储处理器可执行指令的存储器；其中，所述处理器被配置为执行所述方法的全部或部分步骤。

根据本发明实施例的第三方面，提供一种非临时性计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现项所述方法的全部或部分步骤。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍。

图1是根据本发明一实施例提供的一种精确识别振动信号主频的方法流程图。

图2是采样得到原始信号X(t)后截取k＝4500点后的信号X₁(t)的函数图像。

图3是分别截取长度为k＝4100，k＝4300，k＝4500，k＝4700，k＝4900时的X₁(t)序列后进行压缩得到的固定长度信号G(t)的函数图像。

图4是对k＝4500时信号G(t)的4096点快速傅里叶变换频谱图。

图5是本发明实施例提供的一种精确识别振动信号主频的方法中分支定界法剪枝搜索的过程示意图。

图6是使用分支定界法快速搜索ratio(k)的函数图像。

图7是对原始信号分别截取点数为k＝4455、k＝4500、k＝4620、k＝4680、k＝4720、k＝4765的信号再经过压缩后进行FFT的频谱图。

具体实施方式

在基于频域的信号主频估计算法中，由于被测信号的采样频率限制以及时域离散化和频域离散化进行傅里叶变换，特别是当信号采样长度为主频信号周期的非整数倍时，信号频谱在主频附近产生频率泄漏和频率拖尾效应，导致信号频率估计不准确。因此对主频的估计，一方面可以通过截取合适的信号长度进行傅里叶变换，最大程度抑制频率泄露的影响(如图7所示，不同截取长度，频率泄露程度不同)，提高信号频率识别精度；另一方面运用插补运算等手段提高计算速度，为快速准确的主频估计提供一种有效的方法。

图1是根据本发明一实施例提供的一种精确识别振动信号主频的方法流程图。如图1所示，该方法在测量时，为减少迭代计算的次数，提高计算速度，采用分支定界法对被测信号的主频所在区间进行分支预测，迭代推算最优解所在区间；对被测信号截取不同长度的傅里叶时间窗，为实现N点基4快速傅里叶变换，在时域将截断后的信号采用线性插补压缩算法压缩至N点(N＝4^ⁿ)；以最小的旁瓣主瓣比(旁主比)为评价该点的频率识别准确度的参考依据，旁主比越小则频率泄漏越小，识别越精确，进而可以准确计算出被测信号的主频和频谱幅值；该方法能够以较少的迭代运算次数，精确的寻找到最佳的离散傅里叶时间窗截取长度，有效的抑制信号的频率泄漏和“栅栏效应”，实现信号主频的快速高精度识别。下文通过步骤1～步骤4对本发明进行详细说明。

步骤1：在满足采样定律的条件下，对振动信号以频率fs进行采样，得到被测信号的AD转换数据序列X(t)。

本实施例中，假设被测量的信号为x(t)＝100sin(2πf₀t)+30ξ，仿真中，将信号的频率设为f₀＝330Hz，噪声ξ为(0～1)之间的随机数。以采样频率fs＝5000Hz对待测信号进行等间隔采样，获得5000个原始采样数据组成的时间序列X(t)。

步骤2：在k可取的范围内使用随机数种子生成一个随机起始点k_start，在数据序列X(t)中截取t＝[1:k_start]，长度为k的序列X₁(t)；使用线性插补压缩算法将长度为k的序列X₁(t)压缩为长度为N的序列G(t)。

本实施例中，在k可取的范围内使用随机数种子生成一个随机起始点k_start＝4500，取初始傅里叶时间窗截取长度k＝4500截取X(t)得到序列X₁(t)，如图2所示；将长度为k＝4500的序列X₁(t)压缩至4096点的序列G(t)，其中G(t)与X₁(t)的压缩点的关系满足：

式中，k是信号X(t)被截取的长度，N是信号压缩长度，t是时间序列的自变量，t∈[1:n]，n是原始信号采样的点数，比如采样频率fs是5000Hz，则n也是5000，[]表示取整，{}表示取小数；经过此方法压缩后的序列G(t)与X₁(t)相比具有相同的周期数与幅值，如图3所示。

步骤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)的函数值。

本实施例中，对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。

本实施例中，采用分支定界法以起始点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)。

此时，最小下界位于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所对应的自变量k₀＝4620为被测信号能够被正周期截取的点数，k₀处的N＝4096点离散傅里叶变换频谱图中主频信号的数字角频率为w₀＝14，幅值L(w₀)＝264900，被测信号的频率f＝(14x5000)/4620＝15.15Hz。

在示例性实施例中，还提供一种系统，该系统包括处理器其中，所述处理器被配置为执行所述方法的全部或部分步骤。

在示例性实施例中，还提供了一种非临时性计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现所述方法的全部或部分步骤。例如，所述非临时性计算机可读存储介质可以是ROM、RAM、CD-ROM、磁带、软盘和光数据存储设备等。

Claims

1.一种识别振动信号主频的方法，其特征在于，包括：

步骤2：使用随机数种子生成一个随机起始点k_start，在数据序列X(t)中截取t＝[1:k_start]，长度为k的序列X₁(t)；使用线性插补压缩算法将长度为k的序列X₁(t)压缩为长度为N的序列G(t)，对于G(t)中的任意t∈[1:n]，满足：

式中，t是时间序列的自变量，t∈[1:n]，n是原始信号采样的点数，[]表示取整，{}表示取小数；

2.一种识别振动信号主频的系统，其特征在于，包括：

处理器；

用于存储处理器可执行指令的存储器；

其中，所述处理器被配置为执行权利要求1所述的方法的步骤。

3.一种非临时性计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1所述方法的步骤。