振动信号转速提取-基于STFT

前言

之前曾介绍过阶次分析的步骤，为了能进行角域重采样，需要获得转速信息，如果安装了转速传感器，我们可以通过转速脉冲信号计算出转速，但是在很多实际工业环境中，振动信号采集和转速采集可能不是一个厂家提供的，因此很可能存在对时的问题，这种情况下就需要我们从振动信号中直接提取转速信息了。
目前比较常用的方法是先将时域振动信号转换为时频信号，然后在转频所在的大致区域内进行峰值提取，下面我们介绍基于短时傅里叶变换（STFT）的转速提取方法。

短时傅里叶变换原理

我们知道傅里叶变换（FFT）可以将时域振动信号转换到频域，但它是对整个信号进行全局频率分析，无法捕捉频率随时间变化的细节。因此为了能够在变转速情况下提取到转速，需要使用时频分析方法，STFT就是一种经典的时频分析工具。
它的原理也很简单，对于一整段振动信号，我们只需要在时间轴上滑动固定大小的窗口，然后对每个窗口的数据再单独进行傅里叶变换即可得到不同时间段的频率信息了。这里面一个很重要的参数就是窗口大小的选择，窗口太小，时间分辨率高，能够捕捉信号的快速变化，但是频率分辨率低，无法区分相近的频率成分，窗口太大，频率分辨率高，能够区分相近的频率成分，但是时间分辨率低，无法捕捉信号的快速变化。这里再顺便提一下频率分辨率的概念，它是指频谱上两根谱线的距离，由信号时长决定。
不同的窗函数选择也会对分析结果有影响，常见的窗口函数有：

• 矩形窗：时间分辨率高，但频率分辨率低，频谱泄漏严重。
• 汉宁窗：频率分辨率较高，频谱泄漏较少，适用于大多数应用。
• 汉明窗：类似于汉宁窗，但频谱泄漏更少。
• 高斯窗：提供较好的时间分辨率和频率分辨率平衡。

在Python中，scipy的singal模块提供了短时傅里叶变换的接口：

scipy.signal.stft(x, fs, window, nperseg, noverlap, nfft, detrend, return_oneside, boundary, padded, axis)

x就是输入的信号，fs是采样频率，window是窗口函数，nperseg是窗口大小，noverlap是窗口的重叠大小（类似于滑动窗口的步长，默认为窗口大小的一半）。

转速提取步骤

介绍了短时傅里叶变换的基本原理，就可以正式进入转速提取的环节了。我们仍然用加拿大渥太华的数据集来演示一下，此次使用H-A-1.mat数据。

时频转换

读取数据的代码如下所示：

def load_vibration_data(file_path):
    """加载振动数据"""
    data = loadmat(file_path)
    vibration_data = data['Channel_1'].flatten()
    return vibration_data

STFT时频转换的代码如下所示：

def stft(vibration_data, fs):
    """
    短时傅里叶变换
    :param vibration_data: 振动信号数据
    :param fs: 采样率
    :return:
    """
    # STFT参数
    window_time = 0.5  # 500ms窗口
    overlap_ratio = 0.95  # 95%重叠以获得更高的时间分辨率
    nperseg = int(window_time * fs)
    nperseg = 2 ** int(np.log2(nperseg))  # 确保是2的幂
    noverlap = int(nperseg * overlap_ratio)

    # 执行STFT
    f, t, Zxx = signal.stft(vibration_data, fs=fs, nperseg=nperseg,
                            noverlap=noverlap, window='hann')
    return f, t, Zxx

为了方便演示，这里对于频率做了截断，只选取了0-30区间内的频率，可以很明显地看到在10-30区间内有一条很明显的能量带，实际上就是转速频率所在。

峰值提取

得到时频图之后，就可以提取转速信息了，最简单的就是使用singnal模块的find_peaks方法：

scipy.signal.find_peaks(x, height=None, threshold=None, distance=None, prominence=None, width=None, wlen=None, rel_height=0.5, plateau_size=None)

各主要参数的意义如下：

• x：输入的一维信号数组。
• height：峰值的最小高度。
• threshold：峰值与其相邻样本的最小差值。
• distance：相邻峰值之间的最小水平距离。
• prominence：峰值的最小突出度，突出度是峰值与严格较低（不包括距离小于distance的点）的局部最小值之间的垂直距离，具体可以参照Peak prominences
• width：峰值的最小宽度，峰值周围至少要有这么多（不包括峰值本身）的连续点。

为了进一步地去除噪声的影响，还需要做一定的平滑处理，可以使用signal模块的savgol_filter方法：

scipy.signal.savgol_filter(x, window_length, polyorder)

各主要参数的意义如下：

• x：输入的一维信号数组。
• window_length：滤波窗口长度，取值为奇数，且不能超过原始信号长度。
• polyorder：多项式拟合的参数，其值越大，最后的曲线越平滑。

最终的代码如下所示：

def estimate_speed(f, t, Zxx):
    """
    转速估计
    :param f: 频率
    :param t: 时间
    :param Zxx: 时频幅值矩阵
    :return:
    """
    # 选择转频所在的区间
    freq_mask = (f > 0) & (f <= 30)
    freq_band = f[freq_mask]

    speeds = []
    max_freqs = []

    # 对每个时间窗口进行分析
    for i in range(len(t)):
        spectrum = np.abs(Zxx[freq_mask, i])

        # 峰值检测
        peaks, properties = signal.find_peaks(spectrum,prominence=np.max(spectrum) * 0.1)

        if len(peaks) > 0:
            peak_idx = peaks[np.argmax(properties['prominences'])]
            freq = freq_band[peak_idx]
        else:
            peak_idx = np.argmax(spectrum)
            freq = freq_band[peak_idx]

        max_freqs.append(freq)
        # 转换为RPM
        speed = freq * 60  
        speeds.append(speed)

    speeds = np.array(speeds)

    window_length = min(501, len(speeds) - 2)
    # 保证window_length为奇数
    if window_length % 2 == 0:
        window_length -= 1
    if window_length >= 5:
        speeds = signal.savgol_filter(speeds, window_length, 3)

    return speeds, t

我们可以将提取的转速与真实转速做一个对比：

可以看出差距还是比较小的，但是要想应用到实际工作中，还是需要做一些参数优化及其他工作。

总结

本篇文章介绍了基于短时傅里叶变换的转速提取工作，但整个过程还是较为粗糙，譬如短时傅里叶变换如何平衡时间分辨率和频率分辨率，转速提取也可以使用现有的脊线提取方法替代。