python 振动信号分析(1):时频域指标

时域统计指标计算公式
振动信号原始统计特征分为两类：时域统计特征、频域统计特征。

时域指标

信号的时域特征是通过统计分析信号的各种时域参数、指标的估计或计算得到的，如表所示，分为有量纲参数和无量纲参数两种，其中1-9为有量纲参数和10-15无量纲参数。

python程序

def get_time_domain_features(data):
    '''data为一维振动信号'''
    x_rms = 0
    absXbar = 0
    x_r = 0
    S = 0
    K = 0
    k = 0
    x_rms = 0
    fea = []
    len_ = len(data.iloc[0, :])
    mean_ = data.mean(axis=1)  # 1.均值
    var_ = data.var(axis=1)  # 2.方差
    std_ = data.std(axis=1)  # 3.标准差
    max_ = data.max(axis=1)  # 4.最大值
    min_ = data.min(axis=1)  # 5.最小值
    x_p = max(abs(max_[0]), abs(min_[0]))  # 6.峰值
    for i in range(len_):
        x_rms += data.iloc[0, i] ** 2
        absXbar += abs(data.iloc[0, i])
        x_r += math.sqrt(abs(data.iloc[0, i]))
        S += (data.iloc[0, i] - mean_[0]) ** 3
        K += (data.iloc[0, i] - mean_[0]) ** 4
    x_rms = math.sqrt(x_rms / len_)  # 7.均方根值
    absXbar = absXbar / len_  # 8.绝对平均值
    x_r = (x_r / len_) ** 2  # 9.方根幅值
    W = x_rms / mean_[0]  # 10.波形指标
    C = x_p / x_rms  # 11.峰值指标
    I = x_p / mean_[0]  # 12.脉冲指标
    L = x_p / x_r  # 13.裕度指标
    S = S / ((len_ - 1) * std_[0] ** 3)  # 14.偏斜度
    K = K / ((len_ - 1) * std_[0] ** 4)  # 15.峭度

    fea = [mean_[0],absXbar,var_[0],std_[0],x_r,x_rms,x_p,max_[0],min_[0],W,C,I,L,S,K]
    return fea

频域指标

振动信号频域分析首先需要把信号的时域波形借助离散傅里叶变换转化为频谱信息，公式如下：

式中：$x ( k Δ t ) $为振动信号的采样值； N 为采样点数； $\Delta t$为采样间隔； k为时域离散值的序号。
求得频谱信息后，可根据频域统计指标公式计算相应的值，公式如下：

python程序

def get_fre_domain_features(f,y):
    fre_line_num = len(y)
    p1 = y.mean()
    p2 = math.sqrt(sum((y-p1)**2)/fre_line_num)
    p3 = sum((y-p1)**3)/(fre_line_num*p2**3)
    p4 = sum((y-p1)**4)/(fre_line_num*p2**4)
    p5 = sum(f*y)/sum(y)
    p6 = math.sqrt(sum((f-p5)**2*y)/fre_line_num)
    p7 = math.sqrt(sum(f**2*y)/sum(y))
    p8 = math.sqrt(sum(f**4*y)/sum(f**2*y))
    p9 = sum(f**2*y)/math.sqrt(sum(y)*sum(f**4*y))
    p10 = p6/p5
    p11 = sum((f-p5)**3*y)/(p6**3*fre_line_num)
    p12 = sum((f-p5)**4*y)/(p6**4*fre_line_num)
    p13 = sum(abs(f-p5)*y)/(math.sqrt(p6)*fre_line_num)
    p = [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13]
    return p

补充：---------------------------------------------------------------------
程序中的f是傅里叶变换之后的频率轴，y是傅里叶变换之后的幅值，也就是平时画的频谱图的x轴和y轴。

def nextpow2(x):
    if x == 0:
        return 0 
    else:
        return int(np.ceil(np.log2(x)))

def Do_fft(sig,Fs):#输入信号和采样频率
    xlen = len(sig)
    sig = sig - sig.mean()
    NFFT = 2**nextpow2(xlen)
    yf = np.fft.fft(sig,NFFT)/xlen*2
    yf = abs(yf[0:int(NFFT/2+1)])
    f = Fs/2*np.linspace(0,1,int(NFFT/2+1))
    f = f[:]
    return f,yf
    #频域离散值的序号

f,y = Do_fft(sig, Fs)
p = get_fre_domain_features(f,y)

基于谱线数，采样频率的实现方式

import numpy as np

def FFT_demo(x,N,fs):
 	'''
    :param x: 输入信号
    :param N: 谱线数
    :param fs: 信号采样频率
    :return: 频谱值, 频率
    '''
    if N % 2 > 0:
        N -= 1
    if N > len(x):
        xs = np.append(x,np.zeros(N-len(x)))
    else:
        xs = x[:N]
    xf = np.fft.rfft(xs)/N
    freq = np.linspace(0,fs/2,int(N/2+1))
    xf = np.abs(xf)*2
    return xf, freq

示例Demo

import numpy as np
import matplotlib.pyplot as plt

def FFT_demo(x,N,fs):
    if N % 2 > 0:
        N -= 1
    if N > len(x):
        xs = np.append(x,np.zeros(N-len(x)))
    else:
        xs = x[:N]
    xf = np.fft.rfft(xs)/N
    freq = np.linspace(0,fs/2,int(N/2+1))
    xf = np.abs(xf)*2
    return xf, freq
    
w = 5 # fre of rotational speed
z = 30 # gear teeth
fs = 1024 # vib signal sampling rate
fsw = 5 # rotational speed sampling rate
Time = 1 
f = w * z
t = np.linspace(0,Time-1/fs,int(Time*fs))
x = (1+1*np.sin(2*np.pi*20*t))*np.sin(2*np.pi*f*t)
Amplitude, fre = FFT1(x,1024,fs)
plt.subplot(2,1,1)
plt.plot(t,x)
plt.title('original')
plt.ylabel('Amplitude')
plt.subplot(2,1,2)
plt.plot(fre,Amplitude)
plt.title('Frequency')
plt.ylabel('Amplitude')
plt.show()

参考：https://blog.csdn.net/baidu_38963740/article/details/111824612