采用欧式距离进行波形匹配

做项目涉及波形相似性比对，并找到偏移量，重合显示两个波形

红色波形为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

posted @ 2018-11-30 16:18  探索无止境  阅读(1500)  评论(0)    收藏  举报