755 lines
21 KiB
Python
755 lines
21 KiB
Python
|
from math import floor, ceil
|
|||
|
|
import numpy as np
|
|||
|
|
from numpy import real, sqrt, mean, array
|
|||
|
|
from scipy import fftpack
|
|||
|
|
# 频域滤波
|
|||
|
|
from scipy.signal import butter, sosfilt, hilbert, argrelextrema
|
|||
|
|
|
|||
|
|
|
|||
|
|
def fft_spectrum(signal: np.ndarray, fs: float):
|
|||
|
|
"""
|
|||
|
|
傅里叶变换频谱
|
|||
|
|
:param signal:np.ndarray 时序振动信号
|
|||
|
|
:param fs:float 采样频率
|
|||
|
|
:return:
|
|||
|
|
f:np.ndarray, 频率
|
|||
|
|
y:np.ndarray, 幅值
|
|||
|
|
"""
|
|||
|
|
n = len(signal)
|
|||
|
|
f = fs * np.arange(n // 2) / n
|
|||
|
|
fft_data = abs(np.fft.fft(signal))[:int(n / 2)] * 2 / n # 归一化处理，双边频谱反转后需要×2
|
|||
|
|
return f, fft_data
|
|||
|
|
|
|||
|
|
|
|||
|
|
def fft_filter(signal, lpf1, lpf2, fs):
|
|||
|
|
"""
|
|||
|
|
signal: np.ndarray 原始振动信号
|
|||
|
|
fs: float 采样频率
|
|||
|
|
lpf1: float 滤波频率-低频
|
|||
|
|
lpf2: float 滤波频率-高频
|
|||
|
|
:return
|
|||
|
|
y: np.ndarray 滤波信号
|
|||
|
|
"""
|
|||
|
|
yy = fftpack.fft(signal)
|
|||
|
|
m = len(yy)
|
|||
|
|
k = m / fs
|
|||
|
|
for i in range(0, floor(k * lpf1)):
|
|||
|
|
yy[i] = 0
|
|||
|
|
for i in range(ceil(k * lpf2 - 1), m):
|
|||
|
|
yy[i] = 0
|
|||
|
|
y = 2 * real(fftpack.ifft(yy))
|
|||
|
|
return y
|
|||
|
|
# 高通滤波
|
|||
|
|
def fft_filter_high_pass(signal, lpf_high, fs):
|
|||
|
|
"""
|
|||
|
|
signal: np.ndarray 原始振动信号
|
|||
|
|
fs: float 采样频率
|
|||
|
|
lpf_high: float 滤波频率
|
|||
|
|
:return
|
|||
|
|
y: np.ndarray 滤波信号
|
|||
|
|
"""
|
|||
|
|
yy = fftpack.fft(signal)
|
|||
|
|
m = len(yy)
|
|||
|
|
k = m / fs
|
|||
|
|
for i in range(0, floor(k * lpf_high) + 1):
|
|||
|
|
yy[i] = 0
|
|||
|
|
for i in range(ceil(k * (fs - lpf_high) - 1), m):
|
|||
|
|
yy[i] = 0
|
|||
|
|
y = real(fftpack.ifft(yy))
|
|||
|
|
return y
|
|||
|
|
# 低通滤波
|
|||
|
|
def fft_filter_low_pass(signal, lpf_low, fs):
|
|||
|
|
"""
|
|||
|
|
signal: np.ndarray 原始振动信号
|
|||
|
|
fs: float 采样频率
|
|||
|
|
lpf_low: float 滤波频率
|
|||
|
|
:return
|
|||
|
|
y: np.ndarray 滤波信号
|
|||
|
|
"""
|
|||
|
|
yy = fftpack.fft(signal)
|
|||
|
|
m = len(yy)
|
|||
|
|
k = m / fs
|
|||
|
|
for i in range(ceil(k * lpf_low - 1), m):
|
|||
|
|
yy[i] = 0
|
|||
|
|
y = 2 * real(fftpack.ifft(yy))
|
|||
|
|
return y
|
|||
|
|
|
|||
|
|
# 希尔伯特变换
|
|||
|
|
def hilbert_envelop(signal):
|
|||
|
|
"""
|
|||
|
|
输入原始信号，输出包络时域信号
|
|||
|
|
"""
|
|||
|
|
hx = hilbert(signal) # 对信号进行希尔伯特变换。
|
|||
|
|
analytic_signal = signal - hx * 1j # 进行hilbert变换后的解析信号
|
|||
|
|
return np.abs(analytic_signal)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def acc2dis(data: np.ndarray, fs: float):
|
|||
|
|
"""
|
|||
|
|
采用时域积分的方式，将振动加速度信号转化为速度信号和位移信号
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
data: np.ndarray, 振动加速度信号
|
|||
|
|
fs: float, 采样频率
|
|||
|
|
|
|||
|
|
Return
|
|||
|
|
------
|
|||
|
|
s_ifft: np.ndarray, 积分速度信号
|
|||
|
|
d_ifft：np.ndarray, 积分位移信号
|
|||
|
|
"""
|
|||
|
|
n = len(data)
|
|||
|
|
a_mul_dt = data / fs
|
|||
|
|
|
|||
|
|
s = []
|
|||
|
|
total = a_mul_dt[0]
|
|||
|
|
for i in range(n - 1):
|
|||
|
|
total = total + a_mul_dt[i + 1]
|
|||
|
|
s.append(total)
|
|||
|
|
s_fft = np.fft.fft(s)
|
|||
|
|
s_fft[:2] = 0 # 去除直流分量
|
|||
|
|
s_fft[-3:] = 0 # 去除直流分量
|
|||
|
|
s_ifft = np.real(np.fft.ifft(s_fft))
|
|||
|
|
|
|||
|
|
s_mut_dt = s_ifft / fs
|
|||
|
|
d = []
|
|||
|
|
total = s_mut_dt[0]
|
|||
|
|
for i in range(n - 2):
|
|||
|
|
total = total + s_mut_dt[i + 1]
|
|||
|
|
d.append(total)
|
|||
|
|
d_fft = np.fft.fft(d)
|
|||
|
|
d_fft[:2] = 0
|
|||
|
|
d_fft[-3:] = 0
|
|||
|
|
d_ifft = np.real(np.fft.ifft(d_fft))
|
|||
|
|
return s_ifft * 1000, d_ifft * 1000000 # 单位转换
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_vel_pass_rms(signal, fs):
|
|||
|
|
"""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
:return: 通频速度有效值（10-1000Hz）
|
|||
|
|
"""
|
|||
|
|
#x1 = fft_filter(signal, 10, 1000, fs)
|
|||
|
|
x2, x3 = acc2dis(signal, fs)
|
|||
|
|
#x4 = fft_filter(x2, 10, 1000, fs)
|
|||
|
|
vel_pass_rms = np.sqrt(np.mean(x2 ** 2))
|
|||
|
|
return vel_pass_rms
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 低频速度有效值（3-1000Hz）
|
|||
|
|
def calc_vel_low_rms(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
v_rms: 速度有效值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 3, 1000, fs)
|
|||
|
|
x2, x3 = acc2dis(x1, fs)
|
|||
|
|
vel_low_rms = np.sqrt(np.mean(x2 ** 2))
|
|||
|
|
return vel_low_rms
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 加速度有效值（3-10KHz）
|
|||
|
|
def calc_acc_rms(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
a_rms: 加速度有效值
|
|||
|
|
"""""
|
|||
|
|
#x1 = fft_filter(signal, 3,10000, fs)
|
|||
|
|
acc_rms = np.sqrt(np.mean(signal ** 2))
|
|||
|
|
return acc_rms
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 加速度峰值（3-10KHz）
|
|||
|
|
def calc_acc_p(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
a_p: 加速度有效值
|
|||
|
|
"""""
|
|||
|
|
#x1 = fft_filter(signal, 3,10000, fs)
|
|||
|
|
x2 = sorted(abs(signal), reverse=True)
|
|||
|
|
x3 = x2[0:100]
|
|||
|
|
acc_p = np.mean(x3)
|
|||
|
|
return acc_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 振动冲击值（5K~10KHz）
|
|||
|
|
def calc_vibration_impulse(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
impulse: 振动冲击值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5000, 10000, fs)
|
|||
|
|
x2 = hilbert_envelop(x1)
|
|||
|
|
x3 = fft_filter(x2, 3, 500, fs)
|
|||
|
|
impulse = np.sqrt(np.mean(x3 ** 2))
|
|||
|
|
return impulse
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 加速度峭度指标(3~10KHz)
|
|||
|
|
def calc_acc_kurtosis(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
acc_k: 峭度指标
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 3, 10000, fs)
|
|||
|
|
K = sum(x1 ** 4) / len(x1)
|
|||
|
|
a = sqrt(sum(x1 ** 2) / len(x1))
|
|||
|
|
acc_kurtosis = K / (a ** 4)
|
|||
|
|
return acc_kurtosis
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 加速度歪度指标(3-10KHz)
|
|||
|
|
def calc_acc_skew(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
acc_s: 歪度指标
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 3, 10000, fs)
|
|||
|
|
S = sum(x1 ** 3) / len(x1)
|
|||
|
|
a = sqrt(sum(x1 ** 2) / len(x1))
|
|||
|
|
acc_skew = S / (a ** 3)
|
|||
|
|
return acc_skew
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 速度峰值(3-1KHz)
|
|||
|
|
def calc_vel_p(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动加速度信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
vel_p: 速度峰值
|
|||
|
|
"""""
|
|||
|
|
#x1 = fft_filter(signal, 10, 1000, fs)
|
|||
|
|
x2, x3 = acc2dis(signal, fs)
|
|||
|
|
x4 = sorted(abs(x2), reverse=True)
|
|||
|
|
x5 = x4[0:100]
|
|||
|
|
vel_p = np.mean(x5)
|
|||
|
|
return vel_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_acc(signal, fs):
|
|||
|
|
"""
|
|||
|
|
加速度传感器特征值计算
|
|||
|
|
"""
|
|||
|
|
# 通频速度有效值
|
|||
|
|
vel_pass_rms = calc_vel_pass_rms(signal, fs)
|
|||
|
|
# 低频速度有效值
|
|||
|
|
vel_low_rms = calc_vel_low_rms(signal, fs)
|
|||
|
|
# 加速度有效值
|
|||
|
|
acc_rms = calc_acc_rms(signal, fs)
|
|||
|
|
# 加速度峰值
|
|||
|
|
acc_p = calc_acc_p(signal, fs)
|
|||
|
|
# 振动冲击值
|
|||
|
|
vibration_impulse = calc_vibration_impulse(signal, fs)
|
|||
|
|
# 加速度峭度指标
|
|||
|
|
acc_kurtosis = calc_acc_kurtosis(signal, fs)
|
|||
|
|
# 加速度歪度指标
|
|||
|
|
acc_skew = calc_acc_skew(signal, fs)
|
|||
|
|
# 速度峰值
|
|||
|
|
vel_p = calc_vel_p(signal, fs)
|
|||
|
|
return vel_pass_rms, vel_low_rms, acc_rms, acc_p, vibration_impulse, acc_kurtosis, acc_skew, vel_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 最大正向峰值(5-500Hz)
|
|||
|
|
def calc_max_positive_p(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
max_positive_p: 最大正向峰值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
x2 = []
|
|||
|
|
for i in range(1, len(x1)):
|
|||
|
|
if x1[i - 1] > 0:
|
|||
|
|
k = x1[i - 1]
|
|||
|
|
x2.append(k)
|
|||
|
|
max_positive_p = abs(max(x2) - np.mean(x1))
|
|||
|
|
return max_positive_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 最大负向峰值(5-500Hz)
|
|||
|
|
def calc_max_negative_p(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
max_negative_p: 最大负向峰值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
x2 = []
|
|||
|
|
for i in range(1, len(x1)):
|
|||
|
|
if x1[i - 1] < 0:
|
|||
|
|
k = x1[i - 1]
|
|||
|
|
x2.append(k)
|
|||
|
|
max_negative_p = abs(min(x2) - np.mean(x1))
|
|||
|
|
return max_negative_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 监测峰峰值(5-500Hz)
|
|||
|
|
def calc_monitor_pp(signal, fs, L):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
L: 传感器量程
|
|||
|
|
Return
|
|||
|
|
mon_pp: 监测峰峰值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
for i in range(1, len(x1)):
|
|||
|
|
if x1[i - 1] < x1[i]:
|
|||
|
|
x1[i] = x1[i - 1] + (1000 / fs) * (L * 0.05)
|
|||
|
|
elif x1[i - 1] > x1[i]:
|
|||
|
|
x1[i] = x1[i - 1] - x1[i - 1] * 0.63 * (1 / fs)
|
|||
|
|
mon_pp = max(x1) - min(x1)
|
|||
|
|
return mon_pp
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 诊断峰峰值(5-500Hz)
|
|||
|
|
def calc_diagnosis_pp(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
dia_pp: 诊断峰峰值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
dia_pp = max(x1) - min(x1)
|
|||
|
|
return dia_pp
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 峰值(5-500Hz)
|
|||
|
|
def calc_peak(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
Peak: 峰值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
peak = (max(x1) - min(x1)) / 2
|
|||
|
|
return peak
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 峰值因子(5-500Hz)
|
|||
|
|
def calc_peaking_factor(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
Cf: 峰值因子
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
PP = (max(x1) - min(x1)) / 2
|
|||
|
|
rms = np.sqrt(np.mean(x1 ** 2))
|
|||
|
|
peaking_factor = PP / rms
|
|||
|
|
return peaking_factor
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 推导峰值(5-500Hz)
|
|||
|
|
def calc_derive_p(signal, fs):
|
|||
|
|
"""""
|
|||
|
|
signal: 振动位移信号
|
|||
|
|
fs: 采样频率
|
|||
|
|
Return
|
|||
|
|
de_p: 推导峰值
|
|||
|
|
"""""
|
|||
|
|
x1 = fft_filter(signal, 5, 500, fs)
|
|||
|
|
rms = np.sqrt(np.mean(x1 ** 2))
|
|||
|
|
de_p = 1.414 * rms
|
|||
|
|
return de_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_displacement(signal, fs, L):
|
|||
|
|
"""
|
|||
|
|
位移传感器特征值计算
|
|||
|
|
"""
|
|||
|
|
# 最大正向峰值
|
|||
|
|
max_positive_p = calc_max_positive_p(signal, fs)
|
|||
|
|
# 最大负向峰值
|
|||
|
|
max_negative_p = calc_max_negative_p(signal, fs)
|
|||
|
|
# 监测峰峰值
|
|||
|
|
mon_pp = calc_monitor_pp(signal, fs, L)
|
|||
|
|
# 诊断峰峰值
|
|||
|
|
dia_pp = calc_diagnosis_pp(signal, fs)
|
|||
|
|
# 峰值
|
|||
|
|
peak = calc_peak(signal, fs)
|
|||
|
|
# 峰值因子
|
|||
|
|
peaking_factor = calc_peaking_factor(signal, fs)
|
|||
|
|
# 推导峰值
|
|||
|
|
de_p = calc_derive_p(signal, fs)
|
|||
|
|
return max_positive_p, max_negative_p, mon_pp, dia_pp, peak, peaking_factor, de_p
|
|||
|
|
|
|||
|
|
|
|||
|
|
def filter_wave(x: np.ndarray, order: int, wn: int or list, type: str, fs: float, output: str = 'sos'):
|
|||
|
|
"""
|
|||
|
|
振动信号过滤,实现高通，低通，旁通等功能。
|
|||
|
|
:param x: np.ndarray
|
|||
|
|
:param order: int,阶次
|
|||
|
|
:param wn: int or list, if type is lowpass or highpass,the value is int,otherwise the value is list 截止频率
|
|||
|
|
,top and bottom limit
|
|||
|
|
:param type:str,lowpass,highpass,bandpass,bands-top
|
|||
|
|
:param fs:float,frequency
|
|||
|
|
:param output:str,three types: ba,sos,zpk
|
|||
|
|
:return:np.array,过滤的信号
|
|||
|
|
"""
|
|||
|
|
sos = butter(order, wn, type, False, fs=fs, output=output)
|
|||
|
|
return sosfilt(sos, x)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def envelope_detection(x: np.ndarray):
|
|||
|
|
"""
|
|||
|
|
提取振动信号的包络曲线
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
x: np.ndarray, 原始振动信号
|
|||
|
|
y： np.ndarray, 希尔伯特变换变换后的解析信号
|
|||
|
|
Return
|
|||
|
|
------
|
|||
|
|
z: np.ndarray, 包络曲线
|
|||
|
|
"""
|
|||
|
|
x1 = hilbert(x)
|
|||
|
|
y = x - x1 * 1j
|
|||
|
|
z = abs(y)
|
|||
|
|
return z
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_ylb_SWE(signal, fs):
|
|||
|
|
"""
|
|||
|
|
:return ylb_SWE: 应力波能量
|
|||
|
|
"""
|
|||
|
|
f_min = 32000 # 低频截至频率
|
|||
|
|
f_max = 40000 # 高频截至频率
|
|||
|
|
signal_filter_data = filter_wave(signal, 8, [f_min, f_max], "bandpass", fs)
|
|||
|
|
signal_envelope = envelope_detection(signal_filter_data)
|
|||
|
|
ylb_SWE = sum(signal_envelope[signal_envelope > 0])
|
|||
|
|
return ylb_SWE
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_ylb_SWPE(signal, fs):
|
|||
|
|
"""
|
|||
|
|
:return ylb_SWPE: 脉冲能量
|
|||
|
|
"""
|
|||
|
|
f_min = 32000 # 低频截至频率
|
|||
|
|
f_max = 48000 # 高频截至频率
|
|||
|
|
L_times = 0.1 # 最小排序法时有效:人工设置极限阈值
|
|||
|
|
signal_filter_data = filter_wave(signal, 8, [f_min, f_max], "bandpass", fs)
|
|||
|
|
signal_envelope = envelope_detection(signal_filter_data)
|
|||
|
|
signal_envelope_less = signal_envelope[argrelextrema(signal_envelope, np.less)] # 求局部极小值点,np.greatest是极大值点
|
|||
|
|
L_sort = np.sort(signal_envelope_less)
|
|||
|
|
L_Limit = L_sort[round(len(L_sort) * L_times)] # 设置最小阀值
|
|||
|
|
signal_ylb = signal_envelope.copy()
|
|||
|
|
signal_ylb[signal_ylb <= L_Limit] = 0 # 低于阀值设置为0
|
|||
|
|
signal_ylb[0] = 0 # 首元素设置为0
|
|||
|
|
signal_ylb[len(signal_envelope) - 1] = 0 # 尾元素设置为0
|
|||
|
|
signal_ylb[signal_ylb > L_Limit] = 1
|
|||
|
|
signal_dt = np.diff(signal_ylb)
|
|||
|
|
dt_s = [i for i in range(len(signal_dt)) if signal_dt[i] == 1]
|
|||
|
|
dt_x = [j for j in range(len(signal_dt)) if signal_dt[j] == -1]
|
|||
|
|
signal_ylb_SWPE = np.linspace(0, 0, len(dt_s)) # 应力波脉冲能量值(SWPE)赋初值0
|
|||
|
|
for i in range(len(dt_s)):
|
|||
|
|
signal_ylb_SWPE[i] = sum(signal_envelope[dt_s[i] + 1:dt_x[i] + 1] - L_Limit)
|
|||
|
|
ylb_SWPE = sum(signal_ylb_SWPE)
|
|||
|
|
return ylb_SWPE
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_ylb_SWPA(signal, fs):
|
|||
|
|
"""
|
|||
|
|
:return ylb_SWPA: 最大峰高
|
|||
|
|
"""
|
|||
|
|
f_min = 32000 # 低频截至频率
|
|||
|
|
f_max = 40000 # 高频截至频率
|
|||
|
|
L_times = 0.1 # 最小排序法时有效:人工设置极限阈值
|
|||
|
|
signal_filter_data = filter_wave(signal, 8, [f_min, f_max], "bandpass", fs)
|
|||
|
|
signal_envelope = envelope_detection(signal_filter_data)
|
|||
|
|
signal_envelope_less = signal_envelope[argrelextrema(signal_envelope, np.less)]
|
|||
|
|
L_sort = np.sort(signal_envelope_less)
|
|||
|
|
L_Limit = L_sort[round(len(L_sort) * L_times)] # 设置最小阀值
|
|||
|
|
signal_ylb = signal_envelope.copy()
|
|||
|
|
signal_ylb[signal_ylb <= L_Limit] = 0 # 低于阀值设置为0
|
|||
|
|
signal_ylb[0] = 0 # 首元素设置为0
|
|||
|
|
signal_ylb[len(signal_envelope) - 1] = 0 # 尾元素设置为0
|
|||
|
|
signal_ylb[signal_ylb > L_Limit] = 1
|
|||
|
|
signal_dt = np.diff(signal_ylb)
|
|||
|
|
dt_s = [i for i in range(len(signal_dt)) if signal_dt[i] == 1]
|
|||
|
|
dt_x = [j for j in range(len(signal_dt)) if signal_dt[j] == -1]
|
|||
|
|
signal_ylb_SWPA = np.linspace(0, 0, len(dt_s)) # 应力波脉冲能量峰值(SWPA)赋初值0
|
|||
|
|
for i in range(len(dt_s)):
|
|||
|
|
signal_ylb_SWPA[i] = max(signal_envelope[dt_s[i]+1:dt_x[i]+1])
|
|||
|
|
ylb_SWPA = max(signal_ylb_SWPA) # 最大峰高 SWPA
|
|||
|
|
return ylb_SWPA
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_ylb(signal):
|
|||
|
|
"""
|
|||
|
|
应力波传感器特征值计算
|
|||
|
|
:return ylb_SWE, ylb_SWPE, ylb_SWPA: 应力波能量 脉冲能量 最大峰高
|
|||
|
|
"""
|
|||
|
|
fs_ylb = 131072
|
|||
|
|
f_min = 32000 # 低频截至频率
|
|||
|
|
f_max = 40000 # 高频截至频率
|
|||
|
|
L_times = 0.1 # 最小排序法时有效:人工设置极限阈值
|
|||
|
|
signal_filter_data = filter_wave(signal, 8, [f_min, f_max], "bandpass", fs_ylb)
|
|||
|
|
signal_envelope = envelope_detection(signal_filter_data)
|
|||
|
|
signal_envelope_less = signal_envelope[argrelextrema(signal_envelope, np.less)]
|
|||
|
|
L_sort = np.sort(signal_envelope_less)
|
|||
|
|
L_Limit = L_sort[round(len(L_sort) * L_times)] # 设置最小阀值
|
|||
|
|
signal_ylb = signal_envelope.copy()
|
|||
|
|
signal_ylb[signal_ylb <= L_Limit] = 0 # 低于阀值设置为0
|
|||
|
|
signal_ylb[0] = 0 # 首元素设置为0
|
|||
|
|
signal_ylb[len(signal_envelope) - 1] = 0 # 尾元素设置为0
|
|||
|
|
signal_ylb[signal_ylb > L_Limit] = 1
|
|||
|
|
signal_dt = np.diff(signal_ylb)
|
|||
|
|
dt_s = [i for i in range(len(signal_dt)) if signal_dt[i] == 1]
|
|||
|
|
dt_x = [j for j in range(len(signal_dt)) if signal_dt[j] == -1]
|
|||
|
|
signal_ylb_SWPA = np.linspace(0, 0, len(dt_s)) # 应力波脉冲能量峰值(SWPA)赋初值0
|
|||
|
|
signal_ylb_SWPE = np.linspace(0, 0, len(dt_s)) # 应力波脉冲能量值(SWPE)赋初值0
|
|||
|
|
for i in range(len(dt_s)):
|
|||
|
|
signal_ylb_SWPA[i] = max(signal_envelope[dt_s[i]+1:dt_x[i]+1])
|
|||
|
|
signal_ylb_SWPE[i] = sum(signal_envelope[dt_s[i] + 1:dt_x[i] + 1] - L_Limit)
|
|||
|
|
ylb_SWE = sum(signal_envelope[signal_envelope > 0])
|
|||
|
|
ylb_SWPE = sum(signal_ylb_SWPE)
|
|||
|
|
ylb_SWPA = max(signal_ylb_SWPA) # 最大峰高 SWPA
|
|||
|
|
return ylb_SWE, ylb_SWPE, ylb_SWPA
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_multiple_frequency(fft_data, ift, sp):
|
|||
|
|
"""
|
|||
|
|
获取给定频率的幅值
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
fft_data: 频谱数据 [频率, 幅值]
|
|||
|
|
ift: 给定频率
|
|||
|
|
sp: 采样频率除以采样点数
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
f_iX:i倍频幅值
|
|||
|
|
"""
|
|||
|
|
try:
|
|||
|
|
f_iX = max(fft_data[floor((ift - 0.5) / sp):ceil((ift + 0.5) / sp)])
|
|||
|
|
except ValueError:
|
|||
|
|
f_iX = 0
|
|||
|
|
return f_iX
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_HS(fft_data, fea, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
:param fft_data: 频谱数据
|
|||
|
|
:param fea: 特征值
|
|||
|
|
:param sp: 采样频率除以采样点数
|
|||
|
|
:param f_st: 起始阶次
|
|||
|
|
:param f_ord: 计算阶次
|
|||
|
|
:return: fea_HS: f_st~(f_st+f_ord-1)倍频幅值
|
|||
|
|
"""
|
|||
|
|
fea_HS = []
|
|||
|
|
for i in range(f_st, f_st + f_ord):
|
|||
|
|
fea_HS.append(calc_multiple_frequency(fft_data, i * fea, sp))
|
|||
|
|
return fea_HS
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_HS(fft_data, fea, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
计算特征值的HS函数
|
|||
|
|
:param fft_data: 频谱数据
|
|||
|
|
:param fea: 特征值
|
|||
|
|
:param sp: 采样频率除以采样点数
|
|||
|
|
:param f_st: 起始阶次
|
|||
|
|
:param f_ord: 计算阶次
|
|||
|
|
:return: HS:
|
|||
|
|
"""
|
|||
|
|
mul = 0.707
|
|||
|
|
fea_HS = calc_fea_HS(fft_data, fea, sp, f_st, f_ord)
|
|||
|
|
HS = sqrt(sum(array(fea_HS) ** 2)) * mul
|
|||
|
|
return HS
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_HRS(fft_data, fea_min, fea_max, sp):
|
|||
|
|
"""
|
|||
|
|
:param fft_data: 频谱数据
|
|||
|
|
:param fea_min: 频率下限
|
|||
|
|
:param fea_max: 频率上限
|
|||
|
|
:param sp: 采样频率除以采样点数
|
|||
|
|
:return: HRS: 能量和
|
|||
|
|
"""
|
|||
|
|
mul = 0.707
|
|||
|
|
fea_X = fft_data[floor((fea_min - 0.5) / sp):ceil((fea_max + 0.5) / sp)]
|
|||
|
|
HRS = sqrt(sum(fea_X ** 2)) * mul
|
|||
|
|
return HRS
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_HCS_up(fft_data, fc, fb, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
获取给定中心频率和边带的上边带能量和
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
fft_data: 频谱数据 [频率, 幅值]
|
|||
|
|
fc: 中心频率
|
|||
|
|
fb: 边带频率
|
|||
|
|
sp: 采样频率除以采样点数
|
|||
|
|
f_st: 起始阶次
|
|||
|
|
f_ord: 计算阶次
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
HRS_up:f_st~(f_st+f_ord-1)倍上边带能量和
|
|||
|
|
"""
|
|||
|
|
mul = 0.707
|
|||
|
|
HCS_up = []
|
|||
|
|
for i in range(f_st, f_st + f_ord):
|
|||
|
|
xb = []
|
|||
|
|
for j in range(1, 6):
|
|||
|
|
xb.append(calc_multiple_frequency(fft_data, i * fc + j * fb, sp))
|
|||
|
|
xb.append(calc_multiple_frequency(fft_data, i * fc, sp))
|
|||
|
|
HCS_up.append(sqrt(sum([x ** 2 for x in xb])) * mul)
|
|||
|
|
return HCS_up
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_HCS_low(fft_data, fc, fb, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
获取给定中心频率和边带的下边带能量和
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
fft_data: 频谱数据 [频率, 幅值]
|
|||
|
|
fc: 中心频率
|
|||
|
|
fb: 边带频率
|
|||
|
|
sp: 采样频率除以采样点数
|
|||
|
|
f_st: 起始阶次
|
|||
|
|
f_ord: 计算阶次
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
HRS_low:f_st~(f_st+f_ord-1)倍下边带能量和
|
|||
|
|
"""
|
|||
|
|
mul = 0.707
|
|||
|
|
HCS_low = []
|
|||
|
|
for i in range(f_st, f_st + f_ord):
|
|||
|
|
xb = []
|
|||
|
|
for j in range(1, 6):
|
|||
|
|
xb.append(calc_multiple_frequency(fft_data, i * fc - j * fb, sp))
|
|||
|
|
xb.append(calc_multiple_frequency(fft_data, i * fc, sp))
|
|||
|
|
HCS_low.append(sqrt(sum([x ** 2 for x in xb])) * mul)
|
|||
|
|
return HCS_low
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_HDS(fft_data, fea, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
:param fft_data: 频谱数据
|
|||
|
|
:param fea: 特征值
|
|||
|
|
:param sp: 采样频率除以采样点数
|
|||
|
|
:param f_st: 起始阶次
|
|||
|
|
:param f_ord: 计算阶次
|
|||
|
|
:return: fea_HS: 1/f_st~1/(f_st+f_ord-1)倍频幅值
|
|||
|
|
"""
|
|||
|
|
fea_HDS = []
|
|||
|
|
for i in range(f_st, f_st + f_ord):
|
|||
|
|
fea_HDS.append(calc_multiple_frequency(fft_data, 1 / i * fea, sp))
|
|||
|
|
return fea_HDS
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_HDS(fft_data, fea, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
:param fft_data: 频谱数据
|
|||
|
|
:param fea: 特征值
|
|||
|
|
:param sp: 采样频率除以采样点数
|
|||
|
|
:param f_st: 起始阶次
|
|||
|
|
:param f_ord: 计算阶次
|
|||
|
|
:return: HDS:
|
|||
|
|
"""
|
|||
|
|
mul = 0.707
|
|||
|
|
fea_HDS = calc_fea_HDS(fft_data, fea, sp, f_st, f_ord)
|
|||
|
|
HDS = sqrt(sum(array(fea_HDS) ** 2)) * mul
|
|||
|
|
return HDS
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_fea_HCR(fft_data, cf, sf, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
计算给定中心频率和边带频率的边带能量比
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
fft_data: 频谱数据 [频率, 幅值]
|
|||
|
|
cf: 中心频率
|
|||
|
|
sf: 边带频率
|
|||
|
|
sp: 采样频率除以采样点数
|
|||
|
|
f_st: 起始阶次
|
|||
|
|
f_ord: 计算阶次
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
fea_HCR:f_st~(f_st+f_ord-1)倍能量比
|
|||
|
|
"""
|
|||
|
|
fea_HCR = []
|
|||
|
|
for i in range(f_st, f_st + f_ord):
|
|||
|
|
xb1 = []
|
|||
|
|
xb2 = []
|
|||
|
|
for j in range(1, 6):
|
|||
|
|
xb1.append(calc_multiple_frequency(fft_data, i * cf - j * sf, sp))
|
|||
|
|
xb2.append(calc_multiple_frequency(fft_data, i * cf + j * sf, sp))
|
|||
|
|
try:
|
|||
|
|
Hb1 = sum([x ** 2 for x in xb1]) # 下边带能量和
|
|||
|
|
except ValueError:
|
|||
|
|
Hb1 = 0
|
|||
|
|
try:
|
|||
|
|
Hb2 = sum([x ** 2 for x in xb2]) # 上边带能量和
|
|||
|
|
except ValueError:
|
|||
|
|
Hb2 = 0
|
|||
|
|
xc = calc_multiple_frequency(fft_data, i * cf, sp) # 给定中心频率幅值
|
|||
|
|
if xc == 0:
|
|||
|
|
H = 0
|
|||
|
|
else:
|
|||
|
|
H = sqrt(Hb1 + Hb2) / xc
|
|||
|
|
fea_HCR.append(H)
|
|||
|
|
return fea_HCR
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calc_HCR(fft_data, cf, sf, sp, f_st, f_ord):
|
|||
|
|
"""
|
|||
|
|
计算给定中心频率和边带频率的边带能量比
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
fft_data: 频谱数据 [频率, 幅值]
|
|||
|
|
cf: 中心频率
|
|||
|
|
sf: 边带频率
|
|||
|
|
sp: 采样频率除以采样点数
|
|||
|
|
f_st: 起始阶次
|
|||
|
|
f_ord: 计算阶次
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
HCR:f_st~(f_st+f_ord-1)倍能量比的平方和开根号
|
|||
|
|
"""
|
|||
|
|
fea_HCR = calc_fea_HCR(fft_data, cf, sf, sp, f_st, f_ord)
|
|||
|
|
HCR = sqrt(sum(array(fea_HCR) ** 2))
|
|||
|
|
return HCR
|
|||
|
|
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
freq = 8192
|
|||
|
|
data = np.loadtxt(r"E:\Workspace\Software\Python\无线传感器算法\特征值计算\09f07f1800158d00-Y.csv", delimiter=',', usecols=(0))
|
|||
|
|
data_list = data.tolist()
|
|||
|
|
data_list = data_list[1:]
|
|||
|
|
|
|||
|
|
print("数据长度：", len(data_list))
|
|||
|
|
|
|||
|
|
#速度有效值
|
|||
|
|
speed_rms_value = calc_vel_pass_rms(data_list, freq)
|
|||
|
|
print("速度有效值：", speed_rms_value)
|
|||
|
|
# 速度峰值
|
|||
|
|
speed_vel_p = calc_vel_p(data_list, freq)
|
|||
|
|
print("速度有峰值：", speed_vel_p)
|
|||
|
|
# 加速度有效值
|
|||
|
|
acc_rms = calc_acc_rms(data_list, freq)
|
|||
|
|
print("加速度有效值：", acc_rms)
|
|||
|
|
# 加速度峰值
|
|||
|
|
acc_p = calc_acc_p(data_list, freq)
|
|||
|
|
print("加速度峰值：", acc_p)
|