Savitzky-Golay filter design
b = sgolay(k,f)
b = sgolay(k,f,w)
[b,g] = sgolay(...)
b = sgolay(k,f) designs
a Savitzky-Golay FIR smoothing filter
b. The polynomial
k must be less than the frame size,
which must be odd. If
the designed filter produces no smoothing. The output,
f matrix whose rows
represent the time-varying FIR filter coefficients. In a smoothing
filter implementation (for example,
(f-1)/2 rows (each an FIR filter) are
applied to the signal during the startup transient, and the first
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
which contains the real, positive-valued weights to be used during
the least-squares minimization.
[b,g] = sgolay(...) returns
g of differentiation filters. Each column
g is a differentiation filter for derivatives
p is the
column index. Given a signal
x of length
you can find an estimate of the
xp, of its middle value from:
xp((f+1)/2) = (factorial(p)) * g(:,p+1)' * x
sgolay to smooth a noisy sinusoid. Compute the resulting first and second derivatives.
N = 4; % Order of polynomial fit F = 21; % Window length [b,g] = sgolay(N,F); % Calculate S-G coefficients dx = 0.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, % Zeroth 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
Find the first and second derivatives using
diff. Compare the 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','S-G Smoothed') legend boxoff ylabel('Sinusoid') subplot(3,1,2) plot([DiffD1',SG1']) legend('Diff','S-G') legend boxoff ylabel('1st derivative') subplot(3,1,3) plot([DiffD2',SG2']) legend('Diff','S-G') legend boxoff ylabel('2nd derivative')
Zoom in to see more detail. Note how
diff amplifies the noise and generates useless results.
for kj = 1:3 subplot(3,1,kj) axis([400 500 -Inf Inf]) end
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.
 Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1996.