轉載請注明出處：https://www.cnblogs.com/bethansy/p/10560341.html
1、波峰識別
序列數據是按照時間進行采集，其中400個點一個周期，一條數據共10個周期，即4000個點。 首先針對序列數據判斷每個周期是否存在波峰，並在存在波峰的情況下進行波峰形狀識別。 其中波峰識別主要遍歷尋找某個值既大於左邊又大於右邊
1.1、參數：
①peak_thre：獲取時域振動數據95%(20/400)分位數值大小，將此值得0.8倍作為波峰閾值
②peak_width：數據中間隔15個點及其以內的波峰被識別為一個波峰
圖1.1 波峰識別邏輯
1.2、算法原理：
data[i-1]<data[i] & data[i]>data[i+1] & data[i]> peak_thre
OR
data[i-1]=data[i] & data[i]> peak_thre
1.3、步驟
（1）平滑處理序列數據
（2）波峰個數識別
1.4、返回值
（1）N：一條時域數據中波峰個數
（2）peak_index：每個波峰對應的索引值
圖1.2 數據平滑處理
算法處理后，共識別出10個波峰，每個波峰對應的索引分別為[37,77,117,157,197,237,277,317,357,397]
2、波峰形狀識別
這里波峰形狀主要是指等腰三角形和直角三角形，三角形原理
1、參數
①peak_index：上一節獲得的每個波峰的索引
②noise_ratio：為便於計算每個三角形的占空比，默認為每圈最大值的0.2或0.3作為濾噪基數，低於以下的值當做0處理。
2、算法原理
如圖2.2所示的兩種形式會被識別為直角三角形：2.2-a（左邊寬度大於右邊寬度，且占空比小於0.2），2.2-b（右邊寬度大於左邊，且占空比小於0.3）。其中占比0.2和0.3是根據大量數據案例統計得到。
2.2-a 2.2-b
圖2.2 識別為直角三角形的兩種形式
3、步驟
（1）將時域數據每一圈故障波形進行降噪處理
（2）每個波形占空比計算
（3）識別每個波峰三角形形狀，統計一條振動數據出現次數最多的三角形形狀作為每個故障波形三角形形狀輸出結果。
4、返回值
（1）flag：三角形形狀判斷輸出
（2）zanbi：每個三角形的占空比
以圖2.2中時序列數據為例，輸出結果為內部缺陷，每個三角形占空比分別為
[ 0.07, 0.0, 0.2, 0.125, 0.125, 0.175, 0.3, 0.4, 0.1, 0.39]
3、代碼實現
# 波峰個數判別
def find_peak_triangleNum(filter_signal, peak_width):
'''
判斷序列數據波形的三角形個數
:param filter_signal: 平滑后的波形
:param step:連續幾個
:param peak_width:尖峰之間的寬距小於peak_width時划分為一個峰，頻域數據一般定義在20；
時域數據三角形識別一般定義在15，太大會濾掉雙峰
:return:
'''
# 判斷是否有凸起
length_data = len(filter_signal)
thre = 0.7*np.percentile(filter_signal, 95) # 設置閾值高度 95%是前400個點的20個波峰點
# 在整個區域內找極值
l = []
for i in range(1,length_data-1):
if filter_signal[i-1] < filter_signal[i] and filter_signal[i]>filter_signal[i+1] and filter_signal[i]>thre:
l.append(i)
elif filter_signal[i] == filter_signal[i-1] and filter_signal[i]>thre:
l.append(i) # 最高點前后可能有相等的情況
CC = len(l) # 統計極值得個數
cou = 0
ll = l.copy()
for j in range(1, CC):
if l[j]-l[j-1] < peak_width: # 此判斷用於將位於同一個峰內的極值點去掉
if l[j] > l[j-1]: # 同一個峰內的數據，將小的值替換成0
ll[j-1] = 0
else:
ll[j] = 0
cou = cou+1
rcou = CC -cou
ll = [i for i in ll if i > 0] # 去掉0的值
peak_index = []
# 找到每個區間內波峰最大值
# 截斷每個區間再求區間最大值的索引
for i in range(len(ll)):
if i == 0:
index_range = np.array(l)[np.array(l) <= ll[i]]
else:
index_range = np.array(l)[(np.array(l)<=ll[i]) &(np.array(l)>ll[i-1])]
# 找到每個區間最大值得索引
peak_index.append(index_range[np.argmax(filter_signal[index_range],axis=0)])
return [rcou,peak_index]
# 三角形形狀判別
def triangle_flag_meiquan(q_flush, n, ratio,max_index):
'''
判斷三角形的形狀，首先每400個點進行一次過濾，過濾的規則是value-max(value)*0.1
:param q_flush:
:param n: 一條振動數據包含的數據點數
:param ratio: 默認0.2或者0.3
:return:
'''
circle_range = [0]
zhankongbi = []
left = []
right = []
# 生成int_num*step 的矩陣，頭尾需要判斷 ,其中q_flush是list類型。
for i in range(len(max_index)-1):
next_index = int((max_index[i+1]-max_index[i])/2)+max_index[i]
circle_range.append(next_index)
q = q_flush[circle_range[i]:circle_range[i+1]]
# 每圈除噪
newq = q-max(q)*ratio
newq[newq < 0] = 0
# 占空比計算三角形形狀判斷
zhankongbi.append((max(max(np.where(newq > 0))) - min(min(np.where(newq > 0)))) /len(newq))
max_index_now = np.argmax(newq)
left.append(max_index_now - min(min(np.where(newq > 0))))
right.append(max(max(np.where(newq>0))) - max_index_now)
# 進行三角形識別，不考慮左邊為0
right_zb = np.array(zhankongbi)[np.where(np.array(right) >= np.array(left))] # 右邊大於等於左邊對應占空比
right_count = np.sum(right_zb <= 0.3) # 右邊大於左邊占空比小於0.3的個數
left_zb = np.array(zhankongbi)[np.where(np.array(left) > np.array(right))] # 左邊大於等於右邊對應占空比
left_count = np.sum(left_zb <= 0.2) # 左邊大於右邊但是占空比小於0.2
if left_count+right_count >= len(zhankongbi)/2:
flag = 2 # 直角三角形
else:
flag = 1 # 等腰三角形
return [flag, zhankongbi]
# 平滑曲線（濾波）
def smooth(a,WSZ):
'''
:param a: 原始數據，NumPy 1-D array containing the data to be smoothed
必須是1-D的，如果不是，請使用 np.ravel()或者np.squeeze()轉化
:param WSZ: smoothing window size needs, which must be odd number
:return:
'''
out0 = np.convolve(a,np.ones(WSZ,dtype=int),'valid')/WSZ
r = np.arange(1,WSZ-1,2)
start = np.cumsum(a[:WSZ-1])[::2]/r
stop = (np.cumsum(a[:-WSZ:-1])[::2]/r)[::-1]
return np.concatenate((start, out0, stop))