Signal Processing Toolbox™

sgolay - Savitzky-Golay filter design

Syntax

b = sgolay(k,f) b = sgolay(k,f,w) [b,g] = sgolay(...)

Description

b = sgolay(k,f) designs a Savitzky-Golay FIR smoothing filter b. The polynomial order k must be less than the frame size, f, which must be odd. If k = f-1, the designed filter produces no smoothing. The output, b, is an f-by-f matrix whose rows represent the time-varying FIR filter coefficients. In a smoothing filter implementation (for example, sgolayfilt), the last (f-1)/2 rows (each an FIR filter) are applied to the signal during the startup transient, and the first (f-1)/2 rows are applied to the signal during the terminal transient. The center row is applied to the signal in the steady state.

b = sgolay(k,f,w) specifies a weighting vector w with length f, which contains the real, positive-valued weights to be used during the least-squares minimization.

[b,g] = sgolay(...) returns the matrix g of differentiation filters. Each column of g is a differentiation filter for derivatives of order p-1 where p is the column index. Given a signal x of length f, you can find an estimate of the p^th order derivative, xp, of its middle value from:

xp((f+1)/2) = (factorial(p)) * g(:,p+1)' * x

Remarks

Savitzky-Golay smoothing filters (also called digital smoothing polynomial filters or least squares smoothing filters) are typically used to "smooth out" a noisy signal whose frequency span (without noise) is large. In this type of application, Savitzky-Golay smoothing filters perform much better than standard averaging FIR filters, which tend to filter out a significant portion of the signal's high frequency content along with the noise. Although Savitzky-Golay filters are more effective at preserving the pertinent high frequency components of the signal, they are less successful than standard averaging FIR filters at rejecting noise when noise levels are particularly high. The particular formulation of Savitzky-Golay filters preserves various moment orders better than other smoothing methods, which tend to preserve peak widths and heights better than Savitzky-Golay.

Savitzky-Golay filters are optimal in the sense that they minimize the least-squares error in fitting a polynomial to each frame of noisy data.

Examples

Use sgolay to smooth a noisy sinusoid and compare the resulting first and second derivatives to the first and second derivatives computed using diff. Notice how using diff amplifies the noise and generates useless results.

N = 4;                 % Order of polynomial fit
F = 21;                % Window length
[b,g] = sgolay(N,F);   % Calculate S-G coefficients

dx = .2;
xLim = 200;
x = 0:dx:xLim-1;

y = 5*sin(0.4*pi*x)+randn(size(x));  % Sinusoid with noise

HalfWin  = ((F+1)/2) -1;
for n = (F+1)/2:996-(F+1)/2,
  % Zero-th derivative (smoothing only)
  SG0(n) =   dot(g(:,1), y(n - HalfWin: n + HalfWin));
  
  % 1st differential
  SG1(n) =   dot(g(:,2), y(n - HalfWin: n + HalfWin));
  
  % 2nd differential
  SG2(n) = 2*dot(g(:,3)', y(n - HalfWin: n + HalfWin))';
end

SG1 = SG1/dx;         % Turn differential into derivative
SG2 = SG2/(dx*dx);    % and into 2nd derivative

% Scale the "diff" results
DiffD1 = (diff(y(1:length(SG0)+1)))/ dx;	
DiffD2 = (diff(diff(y(1:length(SG0)+2)))) / (dx*dx);
    
subplot(3,1,1);
plot([y(1:length(SG0))', SG0'])
legend('Noisy Sinusoid','S-G Smoothed sinusoid')

subplot(3, 1, 2);
plot([DiffD1',SG1'])
legend('Diff-generated 1st-derivative', ...
'S-G Smoothed 1st-derivative')

subplot(3, 1, 3);
plot([DiffD2',SG2'])
legend('Diff-generated 2nd-derivative',...
'S-G Smoothed 2nd-derivative')

Note The figures below are zoomed in each figure window panel to show more detail.

References