找到所有的峰值点计算峰值点的方法很简单：
- 将曲线上的点分为N组，找到每组中的最小值和最大值，并按照X轴的顺序排序。
- 将所有的峰值点连在一起就是所需绘制的曲线的点。
下面是实现这个算法的程序：
def get_peaks(x, y, n, x0=None, x1=None):
if x0 is None:
x0 = x[0]
if x1 is None:
x1 = x[-1]
index0 = np.searchsorted(x, x0) ❶
index1 = np.searchsorted(x, x1, side="right")
step = (index1 - index0) // n
if step == 0:
step = 1
index1 += 2 * step
if index0 < 0:
index0 = 0
if index1 > len(x) - 1:
index1 = len(x) - 1
x = x[index0:index1+1] ❷
y = y[index0:index1+1]
y = y[:len(y)//step*step] ❸
yy = y.reshape(-1, step) ❹
index = np.c_[np.argmin(yy, axis=1), np.argmax(yy, axis=1)] ❺
index.sort(axis=1) ❻
index += np.arange(0, len(y), step).reshape(-1, 1) ❼
index = index.reshape(-1) ❽
return x[index], y[index]
get_peaks()有5个参数，其中x和y是曲线上各点的X轴和Y轴坐标，要求x中的数据是递增的。n为等分的组数，x0和x1为X轴上的显示范围。❶由于x为递增序列，因此可以使用searchsorted()进行二分查找，快速找到x0和x1对应的下标index0、index1。为了让曲线充满X轴的显示范围，稍微增大一点index1的值，并将index0和index1控制在x的下标范围之内。
❷获得x和y中在显示范围之内的部分，❸并将y的长度调整为step的整数倍，step为每个分组的点数。❹最后将y改为二维数组，第1轴的长度为step。❺找到每组的最小值和最大值的下标，❻并对其进行排序。
❼将index的下标加上每个分组的偏移量得到峰值点在x和y中的下标，❽将下标数组改为一维，并用它获得所有峰值点的X轴坐标和Y轴坐标。
我们用下面的程序产生一组测试数据，在数据中存在许多随机的脉冲噪声，我们希望绘图程序不会漏掉这些脉冲。
def get_sample_data(N):
x = np.linspace(-14, 14, N)
y = np.sin(x) * x**3 * np.sin(100*x)
np.random.seed(1)
y[np.random.randint(0, N, 100)] += np.random.randint(100, 300, 100)
return x, y
matplotlib的绘图程序下面我们用get_peak()加速matplotlib的曲线绘制：
import numpy as np
import matplotlib.pyplot as plt
from get_peaks import get_peaks, get_sample_data
x, y = get_sample_data(1000000)
ax = plt.subplot(111)
line, = plt.plot(*get_peaks(x, y, 500))
ax = plt.gca()
def update_data(ax):
x0, x1 = ax.get_xlim()
line.set_data(*get_peaks(x, y, 500, x0, x1)) ❶
ax.figure.canvas.draw_idle() ❷
ax.callbacks.connect('xlim_changed', update_data) ❸
plt.show()
❸当ax的X轴范围改变时，会产生’xlim_changed’事件，我们用update_data()响应此事件。在update_data()中，首先获得X轴的显示范围，❶并调用line.set_data()将曲线的数据点设置为重新计算的峰值点。❷调用draw_idle()，在空闲时重回整个图表。
下面是程序的显示效果：
Chaco的绘图程序Chaco的程序要稍微复杂一些，不过处理方式类似：
from enthought.traits.api import HasTraits, Instance, on_trait_change
from enthought.traits.ui.api import View, Item
from enthought.chaco.api import Plot, ArrayPlotData
from enthought.enable.component_editor import ComponentEditor
from enthought.chaco.tools.api import PanTool, ZoomTool, DragZoom
from get_peaks import get_peaks, get_sample_data
class LinePlot(HasTraits):
plot = Instance(Plot)
data = Instance(ArrayPlotData)
traits_view = View(
Item('plot',editor=ComponentEditor(), show_label=False),
width=500, height=500, resizable=True, title="Chaco Plot")
def __init__(self, **traits):
super(LinePlot, self).__init__(**traits)
self.x, self.y = get_sample_data(1000000)
peak_x, peak_y = get_peaks(self.x, self.y, 500)
self.data = ArrayPlotData(x=peak_x, y=peak_y)
plot = Plot(self.data)
plot.plot(("x", "y"), type="line", color="blue")
plot.title = "sin(x) * x^3"
plot.tools.append(PanTool(plot))
self.drag_tool = DragZoom(plot, drag_button="right",
maintain_aspect_ratio=False)
plot.tools.append(self.drag_tool)
plot.overlays.append(ZoomTool(plot))
self.plot = plot
@on_trait_change("plot.index_range.updated") ❶
def index_range_changed(self):
x0, x1 = self.plot.index_range.low, self.plot.index_range.high ❷
peak_x, peak_y = get_peaks(self.x, self.y, 500, x0, x1)
self.data["x"] = peak_x ❸
self.data["y"] = peak_y ❸
self.drag_tool._history = [] ❹
if __name__ == "__main__":
p = LinePlot()
p.configure_traits()
❶plot.index_range.updated是一个类型为Event的Trait属性，当X轴显示范围改变时，将触发此事件。我们在事件处理方法中，❷获得X轴的显示范围，并更新数据集中的数据。数据集中的数据更新之后，Chaco将自动重回图表。
❹Chaco会记录下拖动工具的拖动历史，长时间拖动会降低程序的反应时间，因此这里将这个历史记录列表清空。
下面是程序的显示效果：
阅读(5885) | 评论(2) | 转发(0) |
