采用欧式距离进行波形匹配
做项目涉及波形相似性比对,并找到偏移量,重合显示两个波形
红色波形为2号波形,为待匹配波形,采集时间上滞后一定的时间,
蓝色波形为1号波形,现在需要将红色波形移动匹配到蓝色波形上
python代码如下:
1 #!/usr/bin/python3 2 # -*- coding: utf-8 -*- 3 4 import numpy as np 5 import matplotlib.pyplot as plt 6 import sys 7 import scipy.signal as signal 8 from scipy.fftpack import fft,ifft 9 from collections import defaultdict 10 from fastdtw import fastdtw 11 from scipy.spatial.distance import euclidean 12 import struct 13 14 15 16 plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签 17 plt.rcParams['axes.unicode_minus']=False #用来正常显示负号 18 19 ''' 20 返回cwf或者crw文件的波形数组 21 ''' 22 23 ''' 24 局放波形数据(.pwf) 25 头文件说明 26 PDWaveform_HeadFile 27 { 28 //total 28bytes 29 double PD_Test_Voltage;//测试放电电压 8bytes 30 double WaveSpeed;//波速 8bytes 31 int PD_pCRange;//局放范围pC为单位 4bytes 32 int PD_Max_pC;//PD最大pC值 4bytes 33 int WaveformLength;//波形文件长度 4bytes 34 } 35 ''' 36 def getDataFromWaveFile(datafile): 37 with open(datafile, 'rb') as f: 38 ss = f.read(28) 39 varss = struct.unpack('<ddIII',ss) 40 PD_pCRange = varss[2] 41 WaveformLength = varss[4] 42 ss = f.read() 43 hexbytes = struct.unpack(str(WaveformLength)+"b",ss) 44 print('hexbytes:' , hexbytes[1:100]) 45 list1= np.array(hexbytes) / 128 * PD_pCRange 46 return (list1,varss) 47 48 def generateWaveByPwf(datafile, cc): 49 wavedata,_ = getDataFromWaveFile(datafile) 50 list1 = wavedata 51 #list1 = butter_lowpass_filter(list1 ,5e6) 52 list1 = list1.tolist() 53 xx=range(len(list1) ) 54 mIndex = list1.index(max(list1)) 55 plt.plot(xx,list1 , linewidth=0.5,c=cc , linestyle="-", label=datafile)# 折线 1 x 2 y 3 color 56 57 plt.scatter([mIndex, ], [list1[mIndex], ], s=50, color='g') #在这点加个蓝色的原点 原点大小50 58 return list1,mIndex 59 60 def getFastLine(segList, targetList): 61 listret=[] 62 len1 = len(segList) 63 len2 = len(targetList) 64 xx1 = range(len2-len1) 65 for ii in xx1: 66 dist = euclidean(segList , targetList[ii:ii+len1]) 67 listret.append(dist) 68 xindex = listret.index(min(listret)) 69 return listret ,xindex 70 71 def plot(list1 , cc='b'): 72 xx1 = range(len(list1)) 73 plt.plot(xx1,list1 , linewidth=0.5,c=cc)# 折线 1 x 2 y 3 color 74 75 def plot1(list1 , cc='b'): 76 xx1 = range(len(list1)) 77 plt.plot(xx1,list1 , linewidth=0.5,c=cc)# 折线 1 x 2 y 3 color 78 xindex = list1.index(min(list1)) 79 plt.scatter([xindex,],[list1[xindex],], 50, color ='blue') 80 81 82 def main(): 83 84 85 plt.figure(1) # 第一张图 86 list1,maxIndex1 = generateWaveByPwf("L1_12.3kV_1nC_401m_20181122160048.pwf" , 'b') 87 list2,maxIndex2 = generateWaveByPwf("L1_12.3kV_1nC_401m_20181122160208F.pwf" , 'r') 88 plt.legend(loc='upper right') 89 print(maxIndex1 , ':' , maxIndex2) 90 lrnum = 50 # 第一个最值的-20,20个点 91 offset = 500 #2×offset个单元进行fastdtw比较 92 list_seg = list1[maxIndex1-offset :maxIndex1+offset ] 93 list_target = list2[maxIndex2-offset-lrnum:maxIndex2+offset+lrnum] 94 95 96 listret , iindex = getFastLine(list_seg , list_target) 97 print( '窗口偏移量:' , iindex) 98 99 diffx = maxIndex2 - maxIndex1 100 101 list2 = list2[(diffx + lrnum - iindex):] 102 print( '最终偏移量:' , diffx + lrnum - iindex) 103 plt.figure(2) # 第二张图 104 plt.subplot(211) 105 plt.title(U'匹配图') 106 plot(list1 , 'r') 107 plot(list2 , 'b') 108 109 plt.subplot(212) 110 plt.title(U"滑动窗欧式距离曲线") 111 plot1(listret , cc='b') 112 113 plt.grid(axis='y') 114 plt.show() 115 print('done.') 116 117 if __name__ == '__main__': 118 main()
lrnum = 50 # 第一个最值的-50,50个点
offset = 500 #2×offset个单元进行比较,滑动窗大小为1000个点
匹配后最终偏移量: 1246192
滑动33个点为欧式距离最小点,为最佳匹配点。如下图:
放大后的匹配效果
达到预期效果!
github地址
https://github.com/shift0ogg/wavematch