一、為什麼使用移動面積演算法
解：常規波峰判定是採用高低閾值的方法進行篩除，但會出現如圖情況。左邊噪聲高於實際波峰（綠色）高度，甚至高於閾值（紅色），會造成波峰高度的誤判等。
二、移動面積演算法的雛形與原理
選定矩形（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