2.1.22.4.4 ocmath_savitsky_golay

Description

Savitzky-Golay smoothing filter.

Syntax

int ocmath_savitsky_golay( const double * pY, double * pYs, uint nSize, int nLeft, int nRight = -1, int nPolyDeg = 2, int nDervOrder = 0, int nPadding = EDGEPAD_EXTRAPOLATE )

Parameters

pY
[input] pointer to Y vector data
pYs
[output] pointer to smoothed data, or derivatives, depending on nDervOrder
nSize
[input] vector size of both pY and pYs
nLeft
[input] number of data points to the left to be used in filter convolution,
nRight
[input] number of data points to the right to be used in filter convolution, default (-1) will assume nLeft.
Total number of points used in the polynomial fits are (nLeft + nRight + 1) and it must be odd, namely nLeft + nRight must be even.
nPolyDeg
[input] The polynomial order. Higher order will preserve sharper features. nPolyDeg must be less then (nLeft + nRight + 1).
nDervOrder
[input] order of derivative desired (0 = smoothing). To generate a fourth derivative, a minimum quartic (order 4) smoothing must be used.
For a third derivative, a minimum cubic (order 3) smoothing is needed. Similarly, a second derivative requires a minimum quadratic (order 2) smoothing.
nPadding
[input] data are padded with nLeft and nRight points on both ends. Possible values are
EDGEPAD_NONE no padding
EDGEPAD_REFLECT pad reflect, end points are repeated such that on the left, [-1] = [0], [-2] = [1], [-3] = [2] and etc
EDGEPAD_REPEAT pad with [0] values the left and with [nSize-1] on the right
EDGEPAD_EXTRAPOLATE linear extrapolation using nLeft points on the left and nRight points on the right
EDGEPAD_PERIODIC pad periodic, [-1] = [nSize -1], [-2] = [nSize - 2]

Return

OE_NOERROR for success

Examples

EX1

//Assume in the current graph, curve's XY data is in the first data plot. This piece
//of code get the XY data of the curve from the first data plot and Savitzky-Golay
//smoothing filter on it. The result is output in a new worksheet and the new curve
//will plot in the original data plot with color red.
void ocmath_savitsky_golay_ex1()
{
GraphLayer gl = Project.ActiveLayer();
if (!gl)
{
out_str("Active layer is not a graph.");
return;
}

//get XY data from the first dataplot
DataPlot dp = gl.DataPlots(0);
DataRange dr;
vector vx, vy;
if(dp.GetDataRange(dr))
{
DWORD dwPlotID;
if(dr.GetData(DRR_GET_DEPENDENT | DRR_NO_FACTORS, 0, &dwPlotID, NULL, &vy, &vx) < 0)
{
printf("get data failed GetData");
return;
}
}
vector vSmooth;
vSmooth.SetSize(vy.GetSize());

ocmath_savitsky_golay(vy, vSmooth, vy.GetSize(), 3, 3, 2, 1); // Left=Right=7, quadratic, 1st order derivative

//new a worksheet to output the result
Worksheet wks;
wks.Create("Smooth");
wks.SetSize(-1, 2);
wks.SetColDesignations("XY");

DataRange drOut;
drOut.Add("X", wks, 0, 0, -1, 0);
drOut.Add("Y", wks, 0, 1, -1, 1);
drOut.SetData(&vSmooth, &vx);

//plot the curve after smoothing filter
int nPlot = gl.AddPlot(drOut, IDM_PLOT_LINE);
dp = gl.DataPlots(nPlot);
dp.SetColor(1);
}

Remark

Savitzky-Golay smoothing filter.

This time-domain smoothing is based on least squares polynomial fitting across a moving window within the data.

This smoothing method preserves higher moments in the data, thus reducing the distortion of essential features of the data like peak heights and line widths in a spectrum. It can be used to just do data smoothing or to calculate derivatives.

Please note that this routine requires uniform spacing in the time domain for the data.

This function does not handle missing values so work with a copy of data with missing values removed if needed.

cf. Numerical Recipes (Press et al. 1992, Sec.14.8)

origin.h