I have used the Savitzky-Golay filter provided by cookbook, but the result is different from the one given by matlab function <tt>sgolayfilt(x,2,15). I attached the numpy code below:<br><br>#! /usr/bin/env python<br># -*- coding:utf-8 -*-<br>
import numpy as np<br><br>def savitzky_golay(y, window_size, order, deriv=0):<br>    r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.<br>    The Savitzky-Golay filter removes high frequency noise from data.<br>
    It has the advantage of preserving the original shape and<br>    features of the signal better than other types of filtering<br>    approaches, such as moving averages techhniques.<br>    Parameters<br>    ----------<br>
    y : array_like, shape (N,)<br>        the values of the time history of the signal.<br>    window_size : int<br>        the length of the window. Must be an odd integer number.<br>    order : int<br>        the order of the polynomial used in the filtering.<br>
        Must be less then `window_size` - 1.<br>    deriv: int<br>        the order of the derivative to compute (default = 0 means only smoothing)<br>    Returns<br>    -------<br>    ys : ndarray, shape (N)<br>        the smoothed signal (or it's n-th derivative).<br>
    Notes<br>    -----<br>    The Savitzky-Golay is a type of low-pass filter, particularly<br>    suited for smoothing noisy data. The main idea behind this<br>    approach is to make for each point a least-square fit with a<br>
    polynomial of high order over a odd-sized window centered at<br>    the point.<br>    Examples<br>    --------<br>    t = np.linspace(-4, 4, 500)<br>    y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)<br>    ysg = savitzky_golay(y, window_size=31, order=4)<br>
    import matplotlib.pyplot as plt<br>    plt.plot(t, y, label='Noisy signal')<br>    plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')<br>    plt.plot(t, ysg, 'r', label='Filtered signal')<br>
    plt.legend()<br>    plt.show()<br>    References<br>    ----------<br>    .. [1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of<br>       Data by Simplified Least Squares Procedures. Analytical<br>       Chemistry, 1964, 36 (8), pp 1627-1639.<br>
    .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing<br>       W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery<br>       Cambridge University Press ISBN-13: 9780521880688<br>    """<br>
    try:<br>        window_size = np.abs(<a href="http://np.int">np.int</a>(window_size))<br>        order = np.abs(<a href="http://np.int">np.int</a>(order))<br>    except ValueError, msg:<br>        raise ValueError("window_size and order have to be of type int")<br>
    if window_size % 2 != 1 or window_size < 1:<br>        raise TypeError("window_size size must be a positive odd number")<br>    if window_size < order + 2:<br>        raise TypeError("window_size is too small for the polynomials order")<br>
    order_range = range(order+1)<br>    half_window = (window_size -1) // 2<br>    # precompute coefficients<br>    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])<br>    m = np.linalg.pinv(b).A[deriv]<br>
    # pad the signal at the extremes with<br>    # values taken from the signal itself<br>    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )<br>    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])<br>
    y = np.concatenate((firstvals, y, lastvals))<br>    return np.convolve( m, y, mode='valid')<br><br>if __name__ == "__main__":<br>    rrstype = np.dtype({'names':['wavelenth', 'rrs'], 'formats':['f','f']})<br>
<br>    record_list = []<br>    record_file = open("200607NO1.txt")<br>    for line in record_file:<br>        item_list = line.strip().split('\t')<br>        record_list.append((item_list[0], item_list[1]))<br>
<br>    record_file.close()<br><br>    record_array = np.array(record_list,dtype=rrstype)<br>    rrs_array = record_array["rrs"]<br>    rrs_array_2 = savitzky_golay(rrs_array, 15, 2)<br>    rrs_array_2.tofile("savitzky-golay-result.txt", "\n")<br>
</tt>