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>