使用python進行階躍函數分析

Question

我一直在反彈這個問題的一些想法，但認爲我會諮詢在線社區，看看是否有更好的選擇。

所以我有階梯狀的函數圖形是這樣的：

而與此我想計算步驟之間的y位移。

正如人們所看到的，這些步驟並不是完全水平的，而是在升級之前取一小部分y值。

所以，問題是：

（1）什麼是「正確」的方式（如果有的話），每次取「水平」的平均y值？我不確定我應該在哪裏使用我最左邊的點和最右邊的每個點 - 所以我可以將這些點之間的數值取平均值，然後達到每個級別的「平均值」，希望這可以使得感。正如人們所看到的，它們並不全都在x中位移相同。最終目標是獲得水平之間的y位移，一旦我獲得每個「水平」的「平均」值，取得差異是微不足道的。

我可能正在考慮採用曲線的導數，並且在每個關卡的左右兩個點上看到其等於零的位置，但我不確定它會起作用，因爲每個關卡都包含點（ dy/dx = 0） - 所以我可以使用一些見解。

謝謝:)噢 - 這必須在Python中完成 - 它不僅僅是這些圖形，而且它們中的很多具有類似的風格，所以代碼必須足夠通用以處理其他類似步驟圖表。

日期文件圖1：http://textuploader.com/5nwsh

數據文件圖2：http://textuploader.com/5nwsv

數據文件圖3：http://textuploader.com/5nwsj

散點圖Python代碼：

import numpy as np 
import matplotlib.pyplot as plt 
import pylab as pl 


data=np.loadtxt('data-file') 
x= data[:,0] 
y=data[:,1] 


pl.plot(data[:,0],data[:,1],'r') 

pl.xlabel('x') 
pl.ylabel('y') 
plt.show() 

Answer 1

讓我們假設您的數據是晶體表面的AFM測量結果（所有單位均爲m），並且您希望獲得臺階高度的晶體。以下將幫助你達成目標。

from __future__ import division 
from ipywidgets import * 
import numpy as np 
import matplotlib.pyplot as p 
%matplotlib inline 

def rotatedata(x,y,a): 
    cosa=np.cos(a) 
    sina=np.sin(a) 
    x = x*cosa-y*sina 
    y = x*sina+y*cosa 
    return x,y 

data=np.loadtxt('plot3.txt') 
data=data.T 
x,y=data[0],data[1] 

def workit(a2): 
    fig=p.figure(num=None, figsize=(18, 16), dpi= 80, facecolor='w', edgecolor='k') 
    p.subplot(511) # , aspect='equal') 
    p.plot(x,y) 

    #what is the slope? 
    m,b = np.polyfit(x, y, 1) 

    x1,y1=rotatedata(x,y, -np.arctan(m)) # rotate data to have flat surface, 
              # interesting for surface roughness 
    p.subplot(512) 
    p.plot(x1,y1) 

    x2,y2=rotatedata(x,y, a2) # rotate data to 'sharpen' histogram 
    p.subplot(513) 
    p.plot(x2,y2) 

    p.subplot(514) 
    p.hist(y2,bins=130)   

    y3=np.diff(y2) 
    p.subplot(515) 
    p.plot(y3) 

    return HTML() 

interact(workit,a2=[-0.002,0.002,0.00001]) 

第一個情節是原始數據。在第二個圖中，我刪除了數據的斜率，以顯示如果您關心計算表面粗糙度的情況下將使用的數據。

第三個圖顯示相同的數據，如果它們旋轉使得所有斜坡都是水平的。

這是如何完成的（交互式滑塊）顯示在第4個圖中，它是旋轉數據的直方圖。您可以簡單地移動滑塊（旋轉數據），直到直方圖具有最大銳度（所有最大值都是最小寬度）。 我用手做了這個（滑塊）沒有一個函數，但有一個簡單的autofocus例程，最大化相鄰值之間絕對差值的總和可以使用。

在最後一個圖中，我顯示了一階導數（現在是多餘的，但水平區域的平均斜率應該以零爲中心作爲一個雙重檢查）。

水晶層的寬度（您要求的高度）現在由直方圖中最大值的距離給出。

我將每個步驟的實際測定結果作爲練習（例如閾值直方圖，然後計算每個峯的重心，然後獲得相鄰峯重心的差異）。

Answer 2

這裏的第二步是自動調整直方圖。我定義了一個函數focus，它告訴我sharp直方圖是如何使用該函數來查找給出最清晰直方圖的角度的。然後我選擇該角度，對直方圖進行閾值並找出每個直方圖峯值的重心。這些地點之間的差異是步驟。如果這些數字以米爲單位，那麼步長約爲4埃。事實證明，在二維平面AFM或白光干涉儀數據中可以使用同樣的想法，並且非常精確地確定階梯高度，不僅僅是晶體，而是納米級塗層或蝕刻深度。

from __future__ import division 
from ipywidgets import * 
import numpy as np 
import matplotlib.pyplot as p 
%matplotlib inline 

def rotatedata(x,y,a): 
    cosa=np.cos(a) 
    sina=np.sin(a) 
    x = x*cosa-y*sina 
    y = x*sina+y*cosa 
    return x,y 


data=np.loadtxt('plot3.txt') 
data=data.T 
x,y=data[0],data[1] 

def rotateAndCheck(a2): 
    x2,y2=rotatedata(x,y, a2) 
    vals,edges=np.histogram(y2,bins=230) 
    focus=np.sqrt(np.sum((np.diff(vals))**2)) 
    return focus 

focus=[] 
amin,amax,astep=-0.01,0.01,0.0001 
for i in np.arange(amin,amax,astep): 
    focus.append(rotateAndCheck(i)) 


fig=p.figure(num=None, figsize=(18, 16), dpi= 80, facecolor='w', edgecolor='k') 


p.subplot(311)  

p.plot(focus,'.-') 
nm=np.argmax(focus) 
angle=amin+astep*nm 


p.subplot(312) 
x2,y2=rotatedata(x,y, angle) 
vals,edges,_=p.hist(y2,bins=230) 


#now threshold 
p.subplot(313) 
vals[vals<3]=0 
#print len(edges),len(vals) 
deltaedge=edges[1]-edges[0] 
#print deltaedge 
#p.bar(edges[:-1],vals,0.05e-10) 
p.bar(np.arange(len(vals)),vals,0.05e-10) 
p.show() 

# now you go through the histogram from left to right, identify each group and compute the center of gravity for each group 
# this could get trickier if the bin size is not well chosen. 

from scipy.ndimage.measurements import center_of_mass 
levels=[] 

for i in range(1,len(edges)-2): 
    if vals[i-1]==0 and vals[i]>0: 
     istart=i 
     #print 'istart: ',istart 
    if vals[i]>0 and vals[i+1]==0: 
     istop=i 
     #print 'istop', istop 
     sum=np.sum(vals[istart:istop+1]) 
     c= center_of_mass(vals[istart:istop+1])[0] 
     level= edges[istart]+c*deltaedge 
     levels.append(level) 

     #print i, sum 

print 'levels: ',levels   
print 
print 'steps: ' ,np.diff(levels) 

輸出：