Python环境下信号的包络谱分析

本文主要是介绍Python环境下信号的包络谱分析，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

通过信号的频域分析可以获得轴承振动信号在频域内的频率成分分布，即信号的频谱，从中能够提取与轴承故障相关的频率成分幅值与相位，是信号分析中最根本的方法之一。常见的频域分析方法一般有：傅里叶谱分析、解调谱分析和倒频谱分析等。

傅里叶谱分析是可以直观展示信号在整个谱频率域分布状态的方法。通过傅里叶变换后，将信号从时域转换到频域，获得信号在频域上的分布特点，从而观察信号中的频率特征。在对轴承进行故障诊断时，比较轴承振动信号中不同频率成分在傅里叶谱中的分布情况及幅值大小，能够判断轴承是否发生故障及其故障水平。

解调谱分析，即对调制信号的包络解调。当轴承发生局部损伤故障时，系统固有频率会对故障产生的宽带冲击行进调制，此时利用解调分析处理轴承振动信号，可以从信号高频共振频率中把包含轴承故障的低频成分提取出来，从而从包络谱中观察轴承故障特征。常用的解调方法主要包括：Hilbert解调、能量算子解调、平方算子解调等。

倒频谱分析，有时称它为时谱分析，可以用来提取信号频谱中的周期成分，它定义为信号功率谱对数的功率谱。发生局部损伤的轴承元件转动时因相互碰撞而产生周期冲击，此冲击激发轴承系统产生响应，所以此时获得的轴承振动信号是由上述周期冲击与冲击激发的轴承系统响应卷积而来，使谱中产生谐波分量，倒频谱分析则是通过提取上述谐波分量的距离，进而运用到轴承故障诊断中。

注：Hilbert变换通常用来得到解析信号，可以用来对窄带信号进行解包络，并求解信号的瞬时频率。对信号进行Hilbert变换时，会使信号产生一个90°的相位移，并与原信号构成一个解析信号，即为包络信号。Hilbert变换的实质上相当于把原信号通过了一个原始信号和一个信号做卷积的滤波器。可以看成是将原始信号通过一个滤波器。

基于Hilbert变换的包络谱分析简单明了，可解释性较强，有数学理论作为保证，在诸如管道泄漏检测、机械故障诊断、舰船噪声分析、结构损伤识别、油中水分含量检测、电池故障分析、辐射源个体识别、负荷状态识别、污染预测等方面有重要应用。

本项目采用包络谱对轴承振动信号进行分析，运行环境为Python，部分代码如下：

#信号的包络谱分析
#首先导入相关模块
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, fftpack, stats
#设置绘图参数
from matplotlib import rcParams
config = {"font.family": 'serif', # 衬线字体"font.size": 10, # 相当于小四大小"font.serif": ['SimSun'], # 宋体"mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大'axes.unicode_minus': False # 处理负号，即-号
}
rcParams.update(config)
#加载轴承内圈故障数据
def data_acquision(FilePath):"""fun: 从cwru mat文件读取加速度数据param file_path: mat文件绝对路径return accl_data: 加速度数据，array类型"""data = scio.loadmat(file_path)  # 加载mat数据data_key_list = list(data.keys())  # mat文件为字典类型，获取字典所有的键并转换为list类型accl_key = data_key_list[3]  # 获取'X108_DE_time'accl_data = data[accl_key].flatten()  # 获取'X108_DE_time'所对应的值，即为振动加速度信号,并将二维数组展成一维数组return accl_data
#绘制时域波形
file_path = '1730_12k_0.007-InnerRace.mat'
xt = data_acquision(file_path)
plt.figure(figsize=(12,3))
plt.plot(xt)
print(xt)
#做Hilbert变换
ht = fftpack.hilbert(xt)
print(ht)
#计算信号的包络
at = np.sqrt(ht**2+xt**2)   # at = sqrt(xt^2 + ht^2)
#对包络信号做fft变换即为信号的包络谱
sampling_rate = 12000
am = np.fft.fft(at)   # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am)       # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)]  # 获取正频率
#绘制包络谱
plt.plot(freq, am)
#去直流分量
#在0Hz的幅值比较高，使得其它频率幅值较低，不便观察。这种现象叫直流分量，去直流分量方法，y = y-mean(y)
sampling_rate = 12000
at = at - np.mean(at)  # 去直流分量
am = np.fft.fft(at)    # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am)        # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)]  # 获取正频率
plt.plot(freq, am)
#使用包络谱在低频部分观察
sampling_rate = 12000
at = at - np.mean(at)  # 去直流分量
am = np.fft.fft(at)    # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am)        # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)]  # 获取正频率
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
print('内圈故障特征频率',bpfi)
print('外圈故障特征频率',bpfo)
print('滚动体故障特征频率',bsf)
print(ftf)
print(fr)
#理论故障特征频率与实际故障特征频率验证
sampling_rate = 12000
at = at - np.mean(at)  # 去直流分量
am = np.fft.fft(at)    # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am)        # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)]  # 获取正频率
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
plt.vlines(x=156.13, ymin=0, ymax=0.2, colors='r')  # 一倍频
plt.vlines(x=156.13*2, ymin=0, ymax=0.2, colors='r')  # 二倍频
#与FFT进行对比分析
sampling_rate = 12000
am = np.fft.fft(xt)   # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am)       # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate)  # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)]  # 获取正频率
plt.plot(freq, am)plt.plot(freq, am)
plt.xlim(0, 500)
plt.vlines(x=156.13, ymin=0, ymax=0.05, colors='r')  # 一倍频
plt.vlines(x=156.13*2, ymin=0, ymax=0.05, colors='r')  # 二倍频#外圈故障数据测试
file_path = '1730_12k_0.007-OuterRace3.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfo)
#滚动体故障数据测试分析
file_path = '1730_12k_0.014-Ball.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bsf)