python处理扭振信号

一、概述

扭转振动是旋转机械轴系一种特殊的振动形式，它本质上是由于轴系存在弹性，当曲轴在以平均速度进行的旋转过程中，各弹性部件间会因各种原因而产生不同大小、不同相位的瞬时速度的起伏，形成沿旋转方向的来回扭动。

那么如何获取扭振信号？又如何进行扭振分析呢？本文章将由模拟方波信号出发，为您详细讲解扭振处理过程，并配备详细的python程序。

二、扭振分析流程

扭振信号无法直接测得，通常是由瞬时转速转换得来。因此，最重要的步骤就是得到瞬时转速。在工程中，常采用齿轮盘配合电涡流传感器或码带配合激光传感器获取方波信号，关于传感器部分不是本篇文章重点，这里不详述。

在获取方波信号后通过瞬时频率求解得到瞬时频率，再有瞬时频率计算扭转角，最后对扭转角进行FFT变换，就可以得到等角度采样下的扭振谐次分析结果，如下图所示。

三、关键算法解释

这里对其中的其中的关键算法进行解释。

3.1 模拟方波信号

本文章用到的方波信号有三角函数叠加生成，具体公式如下：
$$signal=\frac{\pi }{4}\times \sum _{k=0}^{\infty }\frac{sin[(2k-1)\times (2\pi f)\times t]}{(2k-1)\times \pi }\text{\hspace{1em}(1)}$$
其中f为转频

a.简单方波信号

下列代码首先采用了简单方波信号验证瞬时转速计算的正确性，损失转频大小为20Hz，晶振频率（采样频率）为100Hz。

b.简单复杂信号

为了更贴近实际应用，在简单方波信号验证瞬时转速计算正确后，采用复杂方波信号进行扭振信号处理。假设用到的齿轮盘齿数为100个，轴转一圈得到前50个方波信号频率为5000Hz，后50个方波信号频率为5063Hz，晶振频率（采样频率）为400kHz。

3.2 瞬时转速计算

采用上升沿穿越方式确定穿越点位置。具体的穿越时间，通过穿越前后时间点的一次拟合得到，具体见代码。

3.3 扭转角计算

将瞬时转速减去平均转速（平均转速的求法有多种，这里直接求整个采样周期内的平均转速）得到转速波动，在对转速波动求积分得到扭转角。

3.4 谐次分析

关于谐次分析，内容涉及等角度采样内容，才疏学浅，篇幅受限，实在是讲不透，请各位看官自行学习。

四、代码实现

代码中有详细备注，如有疑问或错误，请留言。

# 导入模块
import numpy as np
import matplotlib.pyplot as plt
# 1 定义穿越幅值
pass_value = 0

# 2 构造简单虚拟方波信号
# 2.1 设置采样率
dt = 1 / 100
# 2.2 定义时间序列
t = np.arange(0, 1, dt)
# 2.3 定义拟合k值范围
k = np.arange(1, 100)
# 2.4 定义转频
f = 20
# 2.5 构造方波信号
sig = np.zeros_like(t)
for i, __ in enumerate(t):
    sig[i] = (4.0 / np.pi) * np.sum(np.sin((2 * k - 1) * (2 * np.pi) * f * t[i]) / ((2 * k - 1) * np.pi))

def cal_fre(t, sig, pass_value):
    # 3 计算转频
    # 3.1 计算穿越点位置
    judge = (sig[:-1] <= pass_value) & (sig[1:] >= pass_value)
    _index = np.where(judge == 1)[0]
    # 3.2 计算穿越点时间
    t_pass = []
    for i in _index:
        #如果穿越点上数值为pass_value，则穿越位置确认
        if sig[i] == pass_value:
            t_pass.append(i * dt)
        elif sig[i + 1] == pass_value:
            t_pass.append((i + 1) * dt)
        #如果穿越点上数值均不为pass_value，则穿越位置通过两点直线拟合
        else:
            t_pass.append((i + (-1) * (pass_value - sig[i])/(sig[i] - sig[i + 1])) * dt)

    # 4 计算转频
    # 4.1 计算时间间隔
    t_pass = np.array(t_pass)
    deta_t = t_pass[1:] - t_pass[:-1]
    # 4.2 计算转频
    fre_res = 1 / deta_t
    
    return t_pass, fre_res

t_pass, fre_res = cal_fre(t, sig, pass_value)
# 5 查看效果
plt.figure(1)
plt.plot(t, sig)
plt.xlabel('t(s)')
plt.ylabel('sig')
plt.figure(2)
plt.plot(t_pass[1:], fre_res)
plt.xlabel('t(s)')
plt.ylabel('frequence')
plt.ylim(0,30)
plt.show()
# 6 构造复杂方波信号
# 6.1 设置采样率
dt = 1 / 400000
# 6.2 设置方波重复次数
num = 4
# 6.2 定义拟合k值范围
k = np.arange(1, 100)
# 6.3 定义转频
f1 = 5000
f2 = 5063
# 6.4 构造方波信号,前50个为5000Hz方波，后50个为5063Hz方波
# a. 计算转频f1下方波信号
t1 = np.arange(0, 50 / f1, dt)
sig1 = np.zeros_like(t1)
for i, __ in enumerate(t1):
    sig1[i] = (4.0 / np.pi) * np.sum(np.sin((2 * k - 1) * (2 * np.pi) * f1 * t1[i]) / ((2 * k - 1) * np.pi))
# b. 计算转频f2下方波信号
t2 = np.arange(0, 50 / f2, dt)
sig2 = np.zeros_like(t2)
for i, __ in enumerate(t2):
    sig2[i] = (4.0 / np.pi) * np.sum(np.sin((2 * k - 1) * (2 * np.pi) * f2 * t2[i]) / ((2 * k - 1) * np.pi))
# c. 构造f1,f2混合方波,并重复num次
sig = np.hstack((sig1, sig2))
sig = np.tile(sig, num)
# d. 构造时间序列
t = np.linspace(0, dt * sig.shape[0], sig.shape[0])
t_pass, fre_res = cal_fre(t, sig, 0.1)
# 7 查看效果
plt.figure(1)
plt.plot(t, sig)
plt.xlabel('t(s)')
plt.ylabel('sig')
plt.figure(2)
plt.plot(t_pass[1:], fre_res)
plt.xlabel('t(s)')
plt.ylabel('frequence')
plt.ylim(4500,6000)
plt.show()
# 8 计算扭振信号
# 8.1 计算时间差
deta_t = t_pass[1:] - t_pass[:-1]
# 8.2 计算瞬时扭转角
w = 2 * np.pi *(fre_res / 100 - fre_res.mean() / 100)
# 8.3 计算随时扭转角度
theta = deta_t * w
for i in range(theta.shape[0] - 1):
    theta[i + 1] = theta[i + 1] + theta[i]
# 8.4 查看计算效果
plt.figure(5)
plt.plot(t_pass[1:], theta)
plt.xlabel('t(s)')
plt.ylabel('degree')
plt.show()
# 9 阶次分析
# 9.1 计算FFT
fft_res = np.fft.fft(theta) / (theta.shape[0] / 2)
fft_res = np.fft.fftshift(fft_res)[int(fft_res.shape[0] / 2):]
# 9.2 获取阶次数
ranks = np.arange(fft_res.shape[0]) * (100 / theta.shape[0])
# 9.3 查看计算效果
plt.plot(ranks[1:], fft_res[1:])
plt.xlabel('ranks')
plt.ylabel('amplitude')
plt.show()