Savitzky-Golay filter design
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));
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
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')
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')
order— Polynomial order
Polynomial order, specified as a positive integer.
framelen— Frame length
Frame length, specified as a positive odd integer. The value of
framelen must be greater than
order=framelen-1, then the designed filter produces no
weights— Weighing 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
b— Time-varying FIR filter coefficients
Time-varying FIR filter coefficients, specified as a
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 of differentiation filters, specified as a matrix. Each column of
g is a differentiation filter for derivatives of order
p is the column index. Given a signal
x of length
framelen, you can find an estimate
pth order derivative,
xp, of its middle value from
(factorial(p)) * g(:,p+1)' * x.
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.
 Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1996.
 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 .