2010-10-28 21 views
2

雖然每個週期的持續時間和內容彼此不同,但我有一個重複信號,其大小每秒重複一次,一些參數內的一點點。我的每個信號數據都有一千個x,y座標。每個週期內的一小部分但重要的數據段已損壞,並且我想用向上的拋物線替換每個已損壞的段。使用numpy在重複信號的一部分內繪製拋物線

對於需要被拋物線替換的每個數據段,我有三個點的x,y座標。頂點/最小值就是其中的一個點。另外兩點是向上的U形的拋物線的左右頂部。換句話說,左上方是該函數域中x值最低的x,y座標對,右上方是該函數域中x值最高的x,y座標對。左上角和右上角的y座標彼此相等,並且是數據段中兩個最高的y值。

如何編寫代碼來繪製這個向上拋物線中的剩餘數據點?請記住,對於每一分鐘的數據,這個函數需要被調用60或70次,並且每次調用這個函數時拋物線的形狀/公式都需要改變,以便說明這三對之間的不同關係在每個拋物線中的x,y座標。

def ReplaceCorruptedDataWithParabola(Xarray, Yarray, LeftTopX, LeftTopY 
            , LeftTopIndex, MinX, MinY, MinIndex 
            , RightTopX, RightTopY, RightTopIndex): 

    # Step One: Derive the formula for the upward-facing parabola using 
    # the following data from the three points: 
     LeftTopX,LeftTopY,LeftTopIndex 
     MinX,MinY,MinIndex 
     RightTopX,RightTopY,RightTopIndex 

    # Step Two: Use the formula derived in step one to plot the parabola in 
    # the places where the corrupted data used to reside: 
    for n in Xarray[LeftTopX:RightTopX]: 
     Yarray[n]=[_**The formula goes here**_] 

    return Yarray 

注:Xarray和Yarray與每個兩個陣列鏈接作爲套X,y座標的索引處數據的每個單列向量。它們都是numpy數組。 Xarray包含時間信息並且不會改變,但Yarray包含信號數據,包括將被替換爲需要由此函數計算的拋物線數據的損壞段。

回答

6

所以,據我瞭解,你有3點,你想適應拋物線。

通常情況下,使用numpy.polyfit最簡單,但如果您真的擔心速度問題,並且恰好適合3個點,則使用最小二乘法擬合沒有意義。相反,我們有一個確定系統(擬合拋物線到3 x,y點),我們可以用簡單的線性代數得到一個精確的解。

所以,這一切的一切,你會做這樣的事情(大部分是繪製數據):

import numpy as np                    
import matplotlib.pyplot as plt                 

def main(): 
    # Generate some random data 
    x = np.linspace(0, 10, 100) 
    y = np.cumsum(np.random.random(100) - 0.5) 

    # Just selecting these arbitrarly 
    left_idx, right_idx = 20, 50  
    # Using the mininum y-value within the arbitrary range 
    min_idx = np.argmin(y[left_idx:right_idx]) + left_idx 

    # Replace the data within the range with a fitted parabola 
    new_y = replace_data(x, y, left_idx, right_idx, min_idx) 

    # Plot the data 
    fig = plt.figure() 
    indicies = [left_idx, min_idx, right_idx] 

    ax1 = fig.add_subplot(2, 1, 1) 
    ax1.axvspan(x[left_idx], x[right_idx], facecolor='red', alpha=0.5) 
    ax1.plot(x, y)              
    ax1.plot(x[indicies], y[indicies], 'ro')       

    ax2 = fig.add_subplot(2, 1, 2) 
    ax2.axvspan(x[left_idx], x[right_idx], facecolor='red', alpha=0.5) 
    ax2.plot(x,new_y)             
    ax2.plot(x[indicies], y[indicies], 'ro') 

    plt.show() 

def fit_parabola(x, y): 
    """Fits the equation "y = ax^2 + bx + c" given exactly 3 points as two 
    lists or arrays of x & y coordinates""" 
    A = np.zeros((3,3), dtype=np.float) 
    A[:,0] = x**2 
    A[:,1] = x 
    A[:,2] = 1 
    a, b, c = np.linalg.solve(A, y) 
    return a, b, c 

def replace_data(x, y, left_idx, right_idx, min_idx): 
    """Replace the section of "y" between the indicies "left_idx" and 
    "right_idx" with a parabola fitted to the three x,y points represented 
    by "left_idx", "min_idx", and "right_idx".""" 
    x_fit = x[[left_idx, min_idx, right_idx]] 
    y_fit = y[[left_idx, min_idx, right_idx]] 
    a, b, c = fit_parabola(x_fit, y_fit) 

    new_x = x[left_idx:right_idx] 
    new_y = a * new_x**2 + b * new_x + c 

    y = y.copy() # Remove this if you want to modify y in-place 
    y[left_idx:right_idx] = new_y 
    return y 

if __name__ == '__main__': 
    main() 

Example plot

希望幫助有點...

+0

+ 1個不錯的答案。 – zellus 2010-10-30 18:34:11