Pythonによるステップ関数解析

Question

私はこの問題についていくつかのアイデアを跳ね返ってきましたが、オンラインコミュニティと相談して、誰かがより良い選択肢を持っているかどうかを調べると考えました。

だから私はこのような階段状の関数グラフいる：

そして、これでは私がステップ間のy変位を計算したいです。

ステップが完全に水平ではなく、むしろステップアップする前に小さな範囲のy値を取ることがわかります。

そこで質問です：

（1）各「レベル」の平均y値を取るために（1がある場合）、「適切」な方法は何ですか？私は左のポイントと右のポイントをそれぞれのレベルで使うべきであるかどうかわからないので、これらのポイント間の値を取って平均し、各レベルの平均を達成することができますセンス。見ることができるように、それらはすべてxの同じ変位にまたがるわけではありません。究極の目標はレベル間にy変位を得ることです。各レベルの「平均」値を取ったら、その差を取るのは簡単です。

私は、カーブの派生を取って、各レベルの左と右のほとんどのポイントでゼロと等しいところを見ていたかもしれませんが、各レベルがポイントを含んでいるので、 dy/dx = 0） - 私はいくつかの洞察力を使用することができます。

ありがとうございます：ああ - これはPythonで行う必要があります - そして、これらのグラフだけではなく、同様のスタイルを持つ多くのものですので、コードは他の階段状グラフ。グラフ3のためのhttp://textuploader.com/5nwsv

データファイル：グラフ2用http://textuploader.com/5nwsh

データファイル：グラフ1のため

日ファイル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

ましょうあなたのデータが結晶表面（すべての単位はm）の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 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]) 

最初のプロットは、あなたの生データです。 2番目のプロットでは、データの傾きを削除して、表面粗さの計算に気をつけた場合に使用するデータを表示しました。

3つ目のプロットは、すべての傾きが水平になるように回転させると同じデータを示しています。

回転したデータのヒストグラムである第4プロットには、これがどのように行われたか（対話スライダー）が示されています。ヒストグラムが最大のシャープネス（すべての最大値が最小幅です）になるまで、単純にスライダを動かす（データを回転させます）。 私はこれを手動で（スライダー）関数ではなく単純なautofocusルーチンで行いました。隣り合った値の絶対差の合計を最大にすることはOKです。

最後のプロットでは、1次導関数を示します（現在は冗長ですが、水平領域の平均スロープは、ダブルチェックとしてゼロ付近を中心にする必要があります）。

ヒストグラムの最大値の距離によって、結晶層の幅（求めていたステップの高さ）が表示されるようになりました。

各ステップの実際の決定は練習問題として残しておきます（ヒストグラムをスレッシュホールドし、各ピークの重心を計算してから、隣接するピーク重心の差を求めるなど）。

Answer 2

第2ステップは、ヒストグラムの自動レベリングです。私は関数focusを教えてくれるのですが、sharpヒストグラムがどのようになっているのですか？その関数を使って、最も鮮明なヒストグラムを与える角度を見つけることができます。私はその角度を選択し、ヒストグラムを閾値化し、各ヒストグラムのピークの重心を見つける。これらの場所の違いは手順です。数字が実際のメートルであれば、階段は約4オングストロームになります。

2Dからレベルの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) 

出力：我々は傾いかどうroadrunner66 @