sgolay

Savitzky-Golay filter design

Syntax

b = sgolay(order,framelen)

b = sgolay(order,framelen,weights)

[b,g] = sgolay(___)

Description

b = sgolay(order,framelen) designs a Savitzky-Golay FIR smoothing filter with polynomial order order and frame length framelen.

b = sgolay(order,framelen,weights) specifies a weighting vector, weights, 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. You can use these output arguments with any of the previous input syntaxes.

Examples

Savitzky-Golay Smoothing of Noisy Sinusoid

Generate a signal that consists of a 0.2 Hz sinusoid embedded in white Gaussian noise and sampled five times a second for 200 seconds.

dt = 1/5;
t = (0:dt:200-dt)';

x = 5*sin(2*pi*0.2*t) + randn(size(t));

Use sgolay to smooth the signal. Use 21-sample frames and 4th-order polynomials.

order = 4;
framelen = 21;

b = sgolay(order,framelen);

Compute the steady-state portion of the signal by convolving it with the center row of b.

ycenter = conv(x,b((framelen+1)/2,:),'valid');

Compute the transients. Use the last rows of b for the startup and the first rows of b for the terminal.

ybegin = b(end:-1:(framelen+3)/2,:) * x(framelen:-1:1);
yend = b((framelen-1)/2:-1:1,:) * x(end:-1:end-(framelen-1));

Concatenate the transients and the steady-state portion to generate the complete smoothed signal. Plot the original signal and the Savitzky-Golay estimate.

y = [ybegin; ycenter; yend];
plot([x y])
legend('Noisy Sinusoid','S-G smoothed sinusoid')

Savitzky-Golay Differentiation

Generate a signal that consists of a 0.2 Hz sinusoid embedded in white Gaussian noise and sampled four times a second for 20 seconds.

dt = 0.25;
t = (0:dt:20-1)';

x = 5*sin(2*pi*0.2*t)+0.5*randn(size(t));

Estimate the first three derivatives of the sinusoid using the Savitzky-Golay method. Use 25-sample frames and 5th-order polynomials. Divide the columns by powers of dt to scale the derivatives correctly.

[b,g] = sgolay(5,25);

dx = zeros(length(x),4);
for p = 0:3
  dx(:,p+1) = conv(x, factorial(p)/(-dt)^p * g(:,p+1), 'same');
end

Plot the original signal, the smoothed sequence, and the derivative estimates.

plot(x,'.-')
hold on
plot(dx)
hold off

legend('x','x (smoothed)','x''','x''''', 'x''''''')
title('Savitzky-Golay Derivative Estimates')

Input Arguments

order — Polynomial order

positive integer

Polynomial order, specified as a positive integer.

framelen — Frame length

positive odd integer

Frame length, specified as a positive odd integer. The value of framelen must be greater than order. If order=framelen-1, then the designed filter produces no smoothing.

weights — Weighing vector

real positive vector

Weighing vector, specified as a real positive vector. The weighing vector has the same length as framelen and are used to perform least-squares minimization.

Output Arguments

b — Time-varying FIR filter coefficients

matrix

Time-varying FIR filter coefficients, specified as a framelen-by-framelen matrix. In a smoothing filter implementation (for example, sgolayfilt), the last (framelen-1)/2 rows (each an FIR filter) are applied to the signal during the startup transient, and the first (framelen-1)/2 rows are applied to the signal during the terminal transient. The center row is applied to the signal in the steady state.

g — Matrix of differentiation filters

matrix

Matrix of differentiation filters, specified as a matrix. 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 framelen, you can find an estimate of the p^th order derivative, xp, of its middle value from xp((framelen+1)/2) = (factorial(p)) * g(:,p+1)' * x.

Algorithms

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.

You can implement data smoothing to measure a variable that is both slowly varying and also corrupted by random noise. Since nearby points measure nearly the same underlying value, you can replace each data point by a local average of the surrounding data points. 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.

References

[1] Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1996.

[2] Press, William. H, Teukolsky, S. A, Vetterling, W. A, and Flannery, B. P . Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, New York, NY, USA 1992 .

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.