# sgolay

Savitzky-Golay filter design

## Syntax

``b = sgolay(order,framelen)``
``b = sgolay(order,framelen,weights)``
``[b,g] = sgolay(___)``

## Description

example

````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.```

example

````[b,g] = sgolay(___)` returns the matrix `g` of differentiation filters. You can use these output arguments with any of the previous input syntaxes.```

## Examples

collapse all

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 fourth 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')```

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 fifth 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

collapse all

Polynomial order, specified as a positive integer. The value of `order` must be smaller than `framelen`. If `order` = `framelen` – 1, then the designed filter produces no smoothing.

Frame length, specified as a positive odd integer. The value of `framelen` must be greater than `order`.

Weighting vector, specified as a real positive vector. The weighting vector has the same length as `framelen` and is used to perform least-squares minimization.

## Output Arguments

collapse all

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.

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 filters are used to smooth out noisy signals with a large frequency span. Savitzky-Golay smoothing filters tend to filter out less of the signal's high-frequency content than standard averaging FIR filters. However, they are less successful at rejecting noise when noise levels are particularly high.

In general, filtering consists of replacing each point of a signal by some combination of the signal values contained in a moving window centered at the point, on the assumption that nearby points measure nearly the same underlying value. For example, moving average filters replace each data point with the local average of the surrounding data points. If a given data point has k points to the left and k points to the right, for a total window length of L = 2k + 1, the moving average filter makes the replacement

`${x}_{s}\to {\stackrel{^}{x}}_{s}=\frac{1}{L}\sum _{r=-k}^{k}{x}_{s+r}.$`

Savitzky-Golay filters generalize this idea by least-squares fitting an nth-order polynomial through the signal values in the window and taking the calculated central point of the fitted polynomial curve as the new smoothed data point. For a given point, xs,

`$\left[\begin{array}{c}{x}_{s-k}\\ ⋮\\ {x}_{s-1}\\ {x}_{s}\\ {x}_{s+1}\\ ⋮\\ {x}_{s+k}\end{array}\right]=\left[\begin{array}{c}{b}_{0}+{b}_{1}\left({t}_{s}-k\Delta t\right)+{b}_{2}{\left({t}_{s}-k\Delta t\right)}^{2}+\cdots +{b}_{n}{\left({t}_{s}-k\Delta t\right)}^{n}\\ ⋮\\ {b}_{0}+{b}_{1}\left({t}_{s}-1\Delta t\right)+{b}_{2}{\left({t}_{s}-1\Delta t\right)}^{2}+\cdots +{b}_{n}{\left({t}_{s}-1\Delta t\right)}^{n}\\ {b}_{0}+{b}_{1}\left({t}_{s}-0\Delta t\right)+{b}_{2}{\left({t}_{s}-0\Delta t\right)}^{2}+\cdots +{b}_{n}{\left({t}_{s}-0\Delta t\right)}^{n}\\ {b}_{0}+{b}_{1}\left({t}_{s}+1\Delta t\right)+{b}_{2}{\left({t}_{s}+1\Delta t\right)}^{2}+\cdots +{b}_{n}{\left({t}_{s}+1\Delta t\right)}^{n}\\ ⋮\\ {b}_{0}+{b}_{1}\left({t}_{s}+k\Delta t\right)+{b}_{2}{\left({t}_{s}+k\Delta t\right)}^{2}+\cdots +{b}_{n}{\left({t}_{s}+k\Delta t\right)}^{n}\end{array}\right]=\left[\begin{array}{c}{a}_{0}+{a}_{1}\left(-k\right)+{a}_{2}{\left(-k\right)}^{2}+\cdots +{a}_{n}{\left(-k\right)}^{n}\\ ⋮\\ {a}_{0}+{a}_{1}\left(-1\right)+{a}_{2}{\left(-1\right)}^{2}+\cdots +{a}_{n}{\left(-1\right)}^{n}\\ {a}_{0}+{a}_{1}\left(0\right)+{a}_{2}{\left(0\right)}^{2}+\cdots +{a}_{n}{\left(0\right)}^{n}\\ {a}_{0}+{a}_{1}\left(1\right)+{a}_{2}{\left(1\right)}^{2}+\cdots +{a}_{n}{\left(1\right)}^{n}\\ ⋮\\ {a}_{0}+{a}_{1}\left(k\right)+{a}_{2}{\left(k\right)}^{2}+\cdots +{a}_{n}{\left(k\right)}^{n}\end{array}\right]$`

or, in terms of matrices,

`$x=\left[\begin{array}{ccccc}1& -k& {\left(-k\right)}^{2}& \cdots & {\left(-k\right)}^{n}\\ 1& ⋮& ⋮& ⋰& ⋮\\ 1& -2& {\left(-2\right)}^{2}& \cdots & {\left(-2\right)}^{n}\\ 1& -1& {\left(-1\right)}^{2}& \cdots & {\left(-1\right)}^{n}\\ 1& 0& 0& \cdots & 0\\ 1& 1& {1}^{2}& \cdots & {1}^{n}\\ 1& 2& {2}^{2}& \cdots & {2}^{n}\\ 1& ⋮& ⋮& \ddots & ⋮\\ 1& k& {k}^{2}& \cdots & {k}^{n}\end{array}\right]\left[\begin{array}{c}{a}_{0}\\ ⋮\\ {a}_{n}\end{array}\right]\equiv Ha.$`

To find the Savitzky-Golay estimates, use the pseudoinverse of H to compute a and then premultiply by H:

`$\stackrel{^}{x}=H{\left({H}^{T}H\right)}^{-1}{H}^{T}x=Bx.$`

To avoid ill-conditioning, `sgolay` uses the `qr` function to compute an economy-size decomposition of H as H = QR, in terms of which B = QQT.

It is necessary to compute B only once. The Savitzky-Golay estimates for most signal points result from convolving the signal with the center row of B. The result is the steady-state portion of the filtered signal. The first k rows of B yield the initial transient, and the final k rows of B yield the final transient. See `sgolayfilt` for an example. It is possible to improve noise suppression by increasing the window length, but this introduces ringing analogous to the Gibbs phenomenon near any transients.

## References

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

[2] Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. New York: Cambridge University Press, 1992.

[3] Schafer, Ronald W. “What Is a Savitzky-Golay Filter? [Lecture Notes].” IEEE Signal Processing Magazine Vol. 28, Number 4, July 2011, pp. 111–117. https://doi.org/10.1109/MSP.2011.941097.