SignalProcessing
SavitzkyGolayFilter
apply Savitzky-Golay filter to a signal
Calling Sequence
Parameters
Options
Description
Examples
References
Compatibility
SavitzkyGolayFilter( signal, degree, frameradius, options )
SavitzkyGolayFilter( signal, degree, weights, options )
signal
-
1-D rtable or list of real-valued data
degree
nonnegative integer for the degree of the polynomials used for filtering
frameradius
positive integer for the radius of the frame used for filtering
weights
2-D rtable or list of positive real-valued weights for the least-squares fitting
extend: Keyword true or false. Specifies if the signal is to be extended on the left and right by frameradius in order to improve the quality of the filtered signal. The default is false.
extrapolation: Keyword periodic, polynomial, or spline. Specifies how the signal is to be extended when extend=true. The default is polynomial.
filteredplotoptions: List of additional plot options used to create the plot of the filtered signal. The default is [].
plotoptions: List of additional plot options used to create the combined plot of the unfiltered and filtered signals. The default is [].
recurrencevariables: List of three unique names for the index, old value, and new value used when constructing the recurrence relation for the filter. The default is [n,x,y].
splinepoints: Positive integer other than 1 or the keyword infinity. Specifies the maximum number of data points used to perform spline extrapolation. The default is 10.
timerange: Range of the form t1..t2, where t1 and t2 are numeric values. Specifies the time range for the plot and Vector of times. The default is 1..numelems(signal).
unfilteredplotoptions: List of additional plot options used to create the plot of the unfiltered signal. The default is [].
output: The type of output. Let m be the signal size, d be the degree, r be the frame radius, and l=2⁢r+1 be the frame length. The supported options are:
coefficients: The float[8] Vector, of size l, storing the coefficients of the filter.
convolution: The float[8] Vector, of size m−2⁢r when extend=false and m when extend=true, storing the convolution of the filter.
derivatives: The float[8] Matrix, of size m by d, storing the derivatives of the signal up to degree d.
differentiator: The float[8] Matrix, of size l by d+1, storing the differentiator of the filter.
filteredsignal: The float[8] Vector, of size m, storing the filtered signal.
hankel: The float[8] Matrix, of size d+1 by d+1, storing the Hankel Matrix of the filter.
nrr: The numeric value corresponding to the noise-reduction ratio of the filtered signal.
plot: The plot of the unfiltered and filtered signals.
projection: The float[8] Matrix, of size l by l, storing the projection of the filter.
recurrence: The linear recurrence relation corresponding to the filter.
times: The float[8] Vector with the time values determined by timerange.
unfilteredsignal: The float[8] Vector, of size m, storing the unfiltered signal.
vandermonde: The float[8] Matrix, of size l by d+1, storing the Vandermonde Matrix of the filter.
weights: The float[8] Vector, of size l, storing the weights.
record: Returns a record with the previous options.
list of any of the above options: Returns an expression sequence with the corresponding outputs, in the same order. The default is filteredsignal.
The SavitzkyGolayFilter command applies the Savitzky-Golay Filter to a 1-D signal with real-valued data, which is useful for smoothing noisy data and estimating derivatives. The method is also known as the Least-Squares Smoothing, Polynomial Smoothing Filter, and Locally Weighted Scatterplot Smoothing (LOWESS).
The filter works by applying polynomial fitting to a succession of frames/windows containing a fixed number of points. More specifically, for each interior data point of the signal, say xk which corresponds to time tk, the frame consists of the points xk−r,..,xk+r, where r is the frame radius. The total number of points in the frame is the frame length, l=2⁢r+1, and the index k satisfies r+1≤k and k≤m+r, where m is the size of the signal. The filtered value yk of xk is determined by the value of the best-fit polynomial of degree d through the points in the frame at time tk.
To determine the filtered signal, Maple does the following:
Suppose the frame consists of r nodes ni spaced uniformly between −1 and 1.
Let p be the polynomial of degree d having coefficients c with respect to the basis defined by the columns of a Vandermonde Matrix V, with dimensions l by d+1. The elements of V are such that the value of p at node ni is element i of y=V·c. The Vandermonde Matrix used by Maple has elements defined by Vi,j=i−r−1rj−1.
It can be shown that the optimal choice of c is given by the the middle row of the Projection Matrix P=F·VT, where the Differentiator Matrix is given by F=V·H−1 and the Hankel Matrix is given by H=VT·W·V. More specifically, if x is the signal over the frame, the objective function to minimize is given by φ=x−V·cT·W·x−V·c.
For each point of the unfiltered signal, say xk with r+1≤k and k≤m+r, the filtered value yk is found using a frame centered at xk and computing the value of the best-fit polynomial at the center node using the above method. All the filtered values together between indices r+1 and m+r, namely the steady-state signal, are computed as the valid convolution of the unfiltered signal and the coefficients Vector.
To find the filtered signal on the left and right of the steady-state signal, namely the input-on transient and input-off transient, respectively, let u be the first r elements of the unfiltered signal x in reverse order, v be the last r elements of x in reverse order, Q be the first r rows of P in reverse order, and R be the last r rows of P in reverse order. The input-on transient is given by R·u and the input-off transient is given by Q·v.
When extend=false, the original signal is extended on the left and right by r so that the overall signal has size m+2⁢r. This is done so that the resulting steady-state of the filtered signal is of size m, which will be used as the filtered version of the original signal.
The noise-reduction ratio is computed as the sum of the squares of the coefficients. Note that the sum of the coefficients is equal to 1.
The Derivatives Matrix is found by convolving the unfiltered signal with the appropriately weighted columns of the Differentiator Matrix. When extend=false the original unfiltered signal is used with the convolution shape being same. On the other hand, when extend=true, the extended unfiltered signal is used with the convolution shape being valid.
The frame length l and radius r satisfy l=2⁢r+1. Moreover, the size of the signal m must satisfy l≤m, and the degree of the filtering polynomials d must satisfy d<2⁢r.
When a weights container is passed, its size corresponds to l and must be an odd integer larger than 1. Internally, the weights will be normalized with respect to root mean square.
The signal and weights rtables cannot have an indexing function and must use Fortran ordering and rectangular storage.
Maple attempts to coerce signal and weights to an rtable of float[8] datatype, and an error is thrown if this is not possible. For this reason, it is most efficient for the passed containers to be rtables of this datatype.
The SavitzkyGolayFilter command is not thread safe.
with⁡SignalProcessing:
Example 1
X≔GenerateSignal⁡exp⁡−0.5⁢t⁢sin⁡2⁢t,t=0..2⁢π,200,noisedeviation=0.05
d≔3:
r≔25:
SavitzkyGolayFilter⁡X,d,r
SavitzkyGolayFilter⁡X,d,r,output=plot
Example 2
Consider the following signal:
f≔sin⁡t+15⁢cos⁡3⁢t
f≔sin⁡t+cos⁡3⁢t5
t1≔0:
t2≔2⁢π:
m≔200:
T,X≔GenerateSignal⁡f,t=t1..t2,m,noisedeviation=0.2,output=times,signal:
Use of non-constant weights can improve the smoothing:
d≔2:
r≔10:
W≔seq⁡1..r+1,seq⁡r..1,−1:
SavitzkyGolayFilter⁡X,d,r,timerange=t1..t2,plotoptions=title=Unweighted Savitzky-Golay Filtered Signal,output=plot
SavitzkyGolayFilter⁡X,d,W,timerange=t1..t2,plotoptions=title=Weighted Savitzky-Golay Filtered Signal,output=plot
Example 3
The Savitzky-Golay Filter can be used to determine the derivatives of a signal. For example, consider the following signal:
g≔t↦2⋅cos⁡t+sin⁡3⋅t
t1≔−2⁢π:
m≔500:
T,X≔GenerateSignal⁡g,t1..t2,m,includefinishtime=false,output=times,signal:
Here, we will compare the unextended and extended cases with the actual derivatives:
r≔20:
Y1≔SavitzkyGolayFilter⁡X,d,r,timerange=t1..t2,extend=false,output=derivatives:
Y2≔SavitzkyGolayFilter⁡X,d,r,timerange=t1..t2,extend=true,extrapolation=periodic,output=derivatives:
dg≔Vector⁡d,k↦Dk⁡g:
Z≔Matrix⁡m,d,i,j↦dgj⁡Ti,datatype=float8:
As we can see, both versions are accurate but the extended one is more accurate at the ends of the time interval:
dataplot⁡T,Y1,style=line,legend=seq⁡sprintf⁡%s derivative,convert⁡i,english,ordinal,i=1..d,title=Savitzky-Golay Derivatives without Extension,view=DEFAULT,min⁡Z..max⁡Z
dataplot⁡T,Y2,style=line,legend=seq⁡sprintf⁡%s derivative,convert⁡i,english,ordinal,i=1..d,title=Savitzky-Golay Derivatives with Extension,view=DEFAULT,min⁡Z..max⁡Z
Using relative root mean square error, we can quantify the accuracy:
seq⁡RelativeRootMeanSquareError⁡Y1..,j,Z..,j,j=1..d
0.343537779779587826,0.393736670852186055,0.940391637995258156
seq⁡RelativeRootMeanSquareError⁡Y2..,j,Z..,j,j=1..d
0.0104018750855233832,0.158908529185620678,0.130136176920078700
As expected, though, the unextended version is very accurate away from the endpoints:
seq⁡RelativeRootMeanSquareError⁡Y1r+1..−r+1,j,Zr+1..−r+1,j,j=1..d
0.0102865451746982136,0.159079567157564050,0.130110621919648450
Example 4
The linear recurrence relation corresponding to the filter can be returned:
req1≔SavitzkyGolayFilter⁡seq⁡7+−1k,k=1..25,2,2,recurrencevariables=i,u,v,output=recurrence
req1≔vi=−0.0857142857142857⁢ui−2+0.342857142857143⁢ui−1+0.485714285714286⁢ui+0.342857142857143⁢ui+1−0.0857142857142857⁢ui+2
In terms of fractions:
req2≔subsindets⁡req1,float,identify
req2≔vi=−3⁢ui−235+12⁢ui−135+17⁢ui35+12⁢ui+135−3⁢ui+235
Example 5
And finally, an example included merely because it looks cool:
SavitzkyGolayFilter⁡Array⁡1..100,k↦−1k,2,15,output=plot
Orfanidis, Sophocles J. Introduction to Signal Processing. Pearson Education, Inc. (2009).
The SignalProcessing[SavitzkyGolayFilter] command was introduced in Maple 2023.
For more information on Maple 2023 changes, see Updates in Maple 2023.
See Also
CurveFitting[Lowess]
SignalProcessing[Convolution]
SignalProcessing[DifferentiateData]
SignalProcessing[Welch]
Statistics[Lowess]
Download Help Document