通过降低采样点数快速绘制曲线

5874阅读 2评论2012-04-10 HyryStudio
分类:Python/Ruby

找到所有的峰值点

计算峰值点的方法很简单:

  1. 将曲线上的点分为N组,找到每组中的最小值和最大值,并按照X轴的顺序排序。
  2. 将所有的峰值点连在一起就是所需绘制的曲线的点。

下面是实现这个算法的程序:

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(),在空闲时重回整个图表。

下面是程序的显示效果:

peak_matplotlib.png
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会记录下拖动工具的拖动历史,长时间拖动会降低程序的反应时间,因此这里将这个历史记录列表清空。

下面是程序的显示效果:

peak_chaco.png
上一篇:活用Spyder中的界面控件
下一篇:用guiqwt制作实时波形绘图程序

文章评论