一、为什么使用移动面积算法
解：常规波峰判定是采用高低阈值的方法进行筛除，但会出现如图情况。左边噪声高于实际波峰（绿色）高度，甚至高于阈值（红色），会造成波峰高度的误判等。
二、移动面积算法的雏形与原理
选定矩形（mask），此处我设其宽为波峰的1/2，高为波峰最高，面积为S2。通过mask在I-V图中，从左到右移动，计算出每一次移动后mask捕捉到的波峰面积（黄色），黄色面积为S1。
通过ratio = S1 / S2，将得到每一个 mask 位置不一样的黄色面积占比，我称之为 ratio。如图，明显可见最右边的 mask 中，黄色占比大于50%（或0.5）。则阈值考量可以设置为0.5~1之间，具体须看实际数据。
ratio > 0.5 时，则 mask 捕捉到了波峰位置。
而且值得注意的是，图二和图一的峰的横坐标位置不一致，由此可见移动面积算法存在高灵活性的巨大优势。
三、代码注释0~7的细解：
0. 限制范围0.1V~1.1V （排除0.1V附近的干扰噪声）
1. 计算峰最高点位置，坐标（V, Imax）
2. 确定移动矩形（称之为mask） 的长和宽，delta_I 和delta_V
3. 确定 mask 移动的初始电压位置（X轴）为0.1V， 限制范围：0A<I<Imax（峰最高点纵坐标）
面积：
area = (Vn - Vn-1) * (In + In-1)/2 ；
比值： # mask从横坐标左往右移动，分别计算并保存每次的ratio值
ratio = 峰波形在mask 所占面积 / mask 总面积 ；
4. 从保存的所有ratio 值，求最大ratio（最大ratio值大于0.5，则可说明波峰存在）
5.五点法确定峰顶点位置
（1）从保存的所有ratio 值，求ratio>0.7 的数量 # 标准峰型通常为30个，此数量体现峰的宽度，筛除脉冲形状的峰
（2）五点法，前三点上升沿与后三点下降沿，判断峰是否为峰，代码中满足①②时皆为波峰顶点，③存在情况过多且不符。
6. 计算单峰面积（不包括其余噪声，即从峰宽两侧算起） # 面积算法（包括噪声）
此计算公式存在些许不足，有一定偏差，可忽视不用。
7. fitter 滤波。
拟合波形，多项式拟合法。作参考。 # 目前屏蔽，暂时不用
四、全篇代码：
def get_moverectangleratiovalue(volts, diffcrnts, weno, movevspos=0.1, detwidth=0.1, noise_line=0, ratiothreshold=0.40, ratiocntlimit=0): #--Mei
"""Get moving rectangle mask ratio value.(Area of interest is called mask.)
The area occupancy ratio
Args:
volts (list) -- List of voltages (float) in volts
diffcrnts (list) -- List of hybrid-to-base current differences (float) in amps
weno(int) -- we number
detwidth (float) -- The width of the interval used to detect the calculated mask voltage range(Usually 0~0.2V)
movevspos (float) -- The starting voltage value used to move the mask (Usually 0~0.5V)
noise_line (float) -- Below the apparent noise line is noise, and the ampere value of the noise line is used as
the baseline for mask calculation
ratiothreshold (float) -- In the mask, the threshold of the ratio of the waveform area to the total area of the mask
rectangle (threshold size: 0~1) If there is a peak, usually 0.5~1
ratiocntlimit (int) -- Determines the peak width.If the ratio count(>0.5)is greater than the limit,
Returns: result (bool) -- Determine whether there is a peak by this ratio and return the result.
If the ratio is higher than "ratiothreshold" return True,
and lower than "ratiothreshold" return False.
ratio (float*1.00) -- The ratio of the area of the waveform to the area of the mask
ampere_max (float) -- Maximum current value
peakmainarea (float) -- Peak area except for noise range(not the peak range)
"""
# init parameter
result = False
result_code = 0 # Tracing the cause of NEG results (0 ~ 5) 0 means POS result. Other means NEG.
limit_volts = [] # V>0
limit_amps = [] # when V>0
peak_pos = (0, 0) # max peak pos (V, A)
mask_wavearea = 0
mask_waveareamovelist = []
ratio_list = []
list_area = []
ratio_cnt = 0
peakposlis = []
def find_cmax(volts, diffcrnts, detv1=0.1, detv2=1.0):
# Find max diff current value
maxcrnt = -1e20
maxv = 9 # Impossibly large value
for x, volt in enumerate(volts):
if volt < detv1:
continue
if volt > detv2:
break
if diffcrnts[x] > maxcrnt:
maxcrnt = diffcrnts[x]
maxv = volts[x]
return (maxcrnt, maxv)
# Normally we will set the noise to 0 or to 2.5e-8 A(real noise) in order to count.
for index, volt in enumerate(volts):
if volt >= movevspos: # Usually the point before 0.1V is not used as a reference.
limit_volts.append(volt)
limit_amps.append(diffcrnts[index])
if len(limit_amps) <= 0:
limit_amps.append(0)
# current_max = max(limit_amps)
# current_max = get_peakdetectionvalue
# 0. Ensure that 'limit_volts' and 'ratio_list' list have the same number of indexes
current_max, volt_max = find_cmax(volts, diffcrnts)
if current_max < noise_line: # noise line is 0 to 2.5e-8(A)
print("return False: Max ampere is lower than noise line")
result_code = 1
return False, 0.05, 0, 0
# print("===================================================================")
# 1.The first "for" loop is for finding the peak pos
for x, volt in enumerate(volts):
# Now I limit the first voltage start value (volt > 0).
if diffcrnts[x] == current_max and volt >= 0:
# find max peak point pos (volt, ampere)
peak_volt = volts[x]
peak_pos = (peak_volt, current_max)
# 2.The second "for" loop is for finding the valley pos.
# Now that the valley pos is not used, if there is a better way to judge the valley pos ??
for x, volt in enumerate(volts):
if volt > 0 and volt < peak_pos[0] and diffcrnts[x-1] < noise_line and diffcrnts[x] >= noise_line:
valley_pos = (volt, diffcrnts[x])
# Mask is the valley pos(x_min, y_min) between peak pos(x_max,y_max) of rectangle area.
delta_V = detwidth #(peak_pos[0] - valley_pos[0])
# Don't pay attention to the impact of noise line for now.(Reserved for subsequent development.)
delta_A = (peak_pos[1] - noise_line) # eliminate noise about 0.25e-7 or lower.
mask_totalarea = delta_V * delta_A
# print("mask total area = %s" % mask_totalarea)
# 3. The third "for" loop is for calculating every moving wave area on mask.
def mask_algorithm(volts, diffcrnts, noise_line, movevspos):
for x, volt in enumerate(volts):
# Determine the starting point pos of the mask
# checking mask range
mask_Vmax = volt + delta_V # new range from max volt of mask
mask_rect = (volt, noise_line, mask_Vmax, peak_pos[1]) # (Vmin,Amin,Vmax,Amax)
if volt >= movevspos: # movevspos or the first volt.This is defined mask moving start pos.
# count the wave area on mask
mask_wavearea = 0 # reset new wave area is 0
for index, v in enumerate(volts):
# Use this FOR loop to count waveform surface base in rectangular area.
# Only count the waveform area of the rectangular area.
if v >= volt and v <= mask_Vmax:
# Each cycle calculates a slice in the mask,
# and each slice is determined by two adjacent points of VOLT.
ampere2 = diffcrnts[index]
ampere1 = diffcrnts[index - 1]
# Ensure that diffcrnts are larger than the noise line.Otherwise equal to the noise line.
if diffcrnts[index] < noise_line:
ampere2 = noise_line
if diffcrnts[index - 1] < noise_line:
ampere1 = noise_line
if diffcrnts[index] > peak_pos[1]:
ampere2 = peak_pos[1]
if diffcrnts[index - 1] > peak_pos[1]:
ampere1 = peak_pos[1]
# Only count ampere more than noise_line 0.25(e-7)
# area = (V1-V2)*(A2+A1-3*A0)/2; A0 = 0
area = (v - volts[index - 1]) * (ampere2 + ampere1 - (3 * noise_line)) / 2
# print("area = (%s - %s) * (%s + %s -3(%s)) / 2 = %s" % (volts[index], volts[index-1], ampere2, ampere1, noise_line, area))
mask_wavearea += area
list_area.append(area)
# print("every mask area = %s" % mask_wavearea)
# print("********************* every wave area is equaled all slice area ************************")
mask_waveareamovelist.append(mask_wavearea)
# print("mask wave area = %s" % mask_wavearea)
# 3.1 Save the waveform to the mask area ratio
for singlearea in mask_waveareamovelist:
ratio_list.append(round(singlearea / mask_totalarea, 5))
# print(ratio_list)
# print("++++++++++++++++++++++++++ all for loop end +++++++++++++++++++++++++++++++++ ")
# 3. Calculating every moving wave area on mask.
mask_algorithm(volts, diffcrnts, noise_line, movevspos)
# 4.When there is more than one peak or no peak, the judgment result is False.
ratio = max(ratio_list)
def multipeaks_screen(ratio, ratio_cnt, result, result_code):
maxratio = 0
if ratio >= ratiothreshold:
# if someone scale more than ratio threshold(usually 0.5) 计数：比值大于阈值的情况
for rat in ratio_list:
if rat >= ratiothreshold:
ratio_cnt += 1
if ratio_cnt > 30: # ratio_cnt 体现的是峰的宽度，ratio_cnt 越大峰越宽。
result_code = 3
result = False
maxratio = 0
# print("weno = %s, ratio count = %s, max-ratio = %s" % (weno, ratio_cnt, ratio))
# Normally it is 30. You can also use the PRINT above to check the ratio count of the normal peak shape.
if ratio_cnt > ratiocntlimit and ratio_cnt <= 30:
# If the ratio count is greater than the limit, the waveform is too wide.
result = False
maxratio = ratio
# If the waveform is too wide, we need confirm the waveform is whether the waveform is multiple peaks.
if ratio_cnt <= 25 and ratio_cnt > 0:
peak_cnt = 0
volts_xlist = []
for i, v in enumerate(volts):
volts_xlist.append(i)
for x, ratio_x in enumerate(ratio_list):
# If there is not noly one peak(maybe more), I set a limit for cheaking a thing called'AS A PEAK'.
if ratio_x > ratiothreshold: # 五点法确定峰顶点位置
if x > 2 and x+2 <= len(ratio_list)-1 :#and x % 2 == 0:
if (ratio_x - ratio_list[x-1]) >= 0 and (ratio_x - ratio_list[x+1]) >= 0:
if (ratio_list[x-2] - ratio_list[x-1]) < 0 and (ratio_list[x+1] - ratio_list[x+2]) > 0:
# Limit the peak appearing before 0.2V, which is regarded as invalid.
if limit_volts[x] > 0.4 and limit_volts[x] < 1.1:
peak_cnt += 1
maxratio = ratio_x
# peakposlis.append((volts[x], diffcrnts[x])) # log peak pos(Shaped like a peak)
# print("Peak pos = (%s, %s); index = %s" % (volts[x], diffcrnts[x], x))
# print("Peak volt pos= %s" % (limit_volts[x]))
# print("peak cnt = %s" % peak_cnt)
if peak_cnt > 1 or peak_cnt == 0:
result = False
maxratio = 0
result_code = 4
if peak_cnt == 1:
result = True
# else:
# result = True
else:
result_code = 2
return maxratio, ratio_cnt, result, result_code
# 5. Multi peaks screening
ratio, ratio_cnt, result, result_code = multipeaks_screen(ratio, ratio_cnt, result, result_code)
def calculate_area(x, noiseline, volts):
# Area calculation formula
a2 = diffcrnts[x]
a1 = diffcrnts[x - 1]
area = (volts[x] - volts[x - 1]) * (a2 + a1 - (3 * noiseline)) / 2
return area
# 6. Calculate the only one peak area
if result: # 多峰筛除之后，Ture结果进行
# calculate peak area 计算单峰面积
peakposlis.append(peak_pos)
# print(peak_pos)
peakarealis = []
temp_min = peak_pos[1]/2
log_index = 0
if peakposlis != []:
for m, (x, y) in enumerate(peakposlis):
for n, volt in enumerate(volts):
if volt > 0.3:
value = abs(diffcrnts[n] - (y/2)) # half peak high 半峰高度与
if value <= temp_min:
temp_min = value
log_index = n
# print(peakposlis[m][0])
# print(volts[log_index])
min_halfwidth = round(abs(peakposlis[m][0] - volts[log_index]), 4)
# print("min_halfwidth=%s" % min_halfwidth)
x_min = peakposlis[m][0] - (min_halfwidth * 2)
x_max = peakposlis[m][0] + (min_halfwidth * 2)
# print("x_min=%s, x_max=%s" % (x_min,x_max))
peakarea = 0
for index, v in enumerate(volts):
if v >= x_min and v <= x_max:
peakarea += calculate_area(index, noise_line, volts)
# print("peak area = %s" % peakarea)
peakarealis.append(peakarea)
peakmainarea = 0
maxv = 0
for peakarea in peakarealis:
if maxv < peakarea:
maxv = peakarea
peakmainarea = maxv
if peakmainarea == 0:
result = False
result_code = 4
ra_index = ratio_list.index(max(ratio_list, key=abs))
# print("ratio max index = %s " % ra_index)
# print("===========================================================================")
else:
peakmainarea = 0
# 7.
if not result and dglobs.fiter_enable and False: # debug fitting and sum res
# Before fitting debug
befo_maxratio = ratio
befo_res = result
calculate_y = fitter.polynomial_fittodraw(volts, diffcrnts, 12, weno) # only a fiter and draw to show
# Fit the original data and output fitted current data
diffs_fited = fitter.polynomial_fitting(volts, diffcrnts, 15) # complex degree = 15
mask_algorithm(volts, diffs_fited, noise_line, movevspos)
after_maxratio = max(ratio_list)
after_res = None
if after_maxratio > ratiothreshold:
after_res = True
else:
after_res = False
result_code = 5
result = after_res
# print("enable fitter and adjust result")
# getamps_leastsquarefitting(volts, diffcrnts)
# print("-----------------Area-Mei----------------")
# print("we: %s, res: %s, ratio: %s, A max: %s, peak area: %s"
# % (weno, result, ratio, current_max, peakmainarea))
# debug+
# calculate_y = fitter.polynomial_fittodraw(volts, diffcrnts, 12, weno) # only a fiter and draw to show
# time.sleep(1)
# Save volt (value/pos) of the max ampere
# for index, volt in enumerate(volts):
# if diffcrnts[index] == current_max:
# dglobs.allwes_peakpos.append(volt)
# debug+
# print("weno = %s, ratio count = %s, max-ratio = %s" % (weno, ratio_cnt, ratio))
# print("result code = %s" % result_code)
return result, ratio, current_max, peakmainarea, result_code
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
