使用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) 

输出：