[lecture NOTES]
Ronald W. Schafer
What Is a Savitzky-Golay Filter?
ecently, while reading a paper impulse response. The lowpass filters digital filters and a basic knowledge of
on heart rate monitoring obtained by this method are widely matrices.
using an accelerometer [1], I known (in some sectors) as Savitzky-
found myself asking the Golay filters. Savitzky and Golay were LEAST-SQUARES
Rquestion posed by the above interested in smoothing noisy data SMOOTHING OF SIGNALS
title. While searching for the answer, I obtained from chemical spectrum ana- The basic idea behind least-squares poly-
discovered many things that seemed to lyzers, and they demonstrated that least- nomial smoothing is depicted in
be well known to others outside the field squares smoothing reduces noise while Figure 1, which shows a sequence of
3 4
of digital signal processing (DSP) but not maintaining the shape and height of samples x n of a signal as solid dots.
to me. After adding some results of my waveform peaks (in their case, Gaussian- Considering for the moment the group
own, I presented what I had learned at a shaped spectral peaks). In researching of 2M 1 1 samples centered at n 5 0, we
poster session at the 2011 DSP/SPE this topic, I did find some awareness of obtain (by a process to be described) the
Workshop [2]. As I interacted with peo- S-G filters in the signal processing com- coefficients of a polynomial
ple at the poster session, I asked if they munity. Hamming s book [7] has a dis-
N
1 2
had ever heard of Savitzky-Golay (S-G) cussion of the use of least-squares in p n 5 ak nk (1)
a
k50
filters. Given my own ignorance, it was data smoothing, and Orfanidis has a
comforting that only one out of about 20 detailed discussion of S-G filters in his that minimize the mean-squared approx-
had heard of them. What is remarkable book, which is now out of print but imation error for the group of input
about this is that Savitzky and Golay s available for free download [8]. Some samples centered on n 5 0,
paper [3], published in 1964, was researchers have found the peak shape
M
described in 2000 by editors of the jour- preservation property of the S-G filters
1 1 2 3 4 22
EN 5 p n 2 x n
a
nal Analytical Chemistry as number five to be attractive in applications such as
n52M
M N 2
among the top ten papers ever published electrocardiogram processing [1] and the
3 4
5 a ak nk 2 x n b . (2)
a a
in that journal [4]. They stated,  It can basic concept of least-squares polynomi-
n52M k50
be argued that the dawn of computer- al smoothing has been generalized to
controlled analytical chemistry can be two dimensions [5] and applied in pro- The analysis is the same for any other
traced to this article. For this reason, I cessing images such as ultrasound and group of 2M 1 1 input samples. We shall
feel that it could be useful to use the synthetic aperture radar. refer to M as the  half width of the
"Lecture Notes" forum to introduce (or While frequency-domain representa- approximation interval. In Figure 1,
reintroduce) my colleagues in signal tions of S-G filters have been discussed where N 5 2 and M 5 2, the solid curve
processing to the S-G filters. [6], [7], most presentations on S-G filters on the left in Figure 1 is the polynomial
1 2
(e.g., [9], [10]) have emphasized time- p n evaluated on a fine grid between
RELEVANCE domain properties (such as complicated 22 and 12, and the smoothed output
1 2
In their  seminal [4] paper [3], Savitzky closed-form expressions for the impulse value is obtained by evaluating p n at
3 4
and Golay proposed a method of data responses) without reference to such the central point n 5 0. That is, y 0 , the
smoothing based on local least-squares frequency-domain features as passband output at n 5 0, is
polynomial approximation. They showed width or stopband attenuation. Therefore,
3 4 1 2
that fitting a polynomial to a set of input the purpose of this article is to examine y 0 5 p 0 5 a0, (3)
samples and then evaluating the result- S-G filters from the frequency-domain
ing polynomial at a single point within viewpoint and to quantify some of their i.e., the output value is just equal to the
the approximation interval is equivalent frequency-domain properties. 0th polynomial coefficient. In general,
to discrete convolution with a fixed the approximation interval need not be
PREREQUISITES symmetric about the evaluation point.
This article assumes only a familiarity This leads to nonlinear phase filters,
Digital Object Identifier 10.1109/MSP.2011.941097
Date of publication: 15 June 2011 with finite-impulse response (FIR) which can be useful for smoothing at
1053-5888/11/$26.00©2011IEEE IEEE SIGNAL PROCESSING MAGAZINE [111] JULY 2011
lecture NOTES continued
[ ]
The equations in (6) are known as the
normal equations for the least-squares
x (m ) or x (n )
approximation problem. It is impor-
tant to note before proceeding that a
unique solution requires that we have
at least as many data samples as we
have coefficients in the polynomial
approximation. That is, we require
N # 2M. In fact, the equations in (6)
become ill-conditioned if M and N are
large and N is close to 2M. Further-
more, if N 5 2M the polynomial fits
the 2M 1 1 data samples exactly and
0 10 m or n
no smoothing results.
Additional insight can be obtained by
[FIG1] Illustration of least-squares smoothing by locally fitting a second-degree
expressing the equations in (6) in matrix
polynomial (solid line) to five input samples: d denotes the input samples, ~
form. To do this, it is helpful to define a
denotes the least-squares output sample, and 3 denotes the effective impulse
response samples (weighting constants). (The dotted line denotes the polynomial
1 2 1 2 5 6
2M 1 1 by N 1 1 matrix A 5 an,i
approximation to centered unit impulse.)
as the matrix with elements
M an, i 5 ni, 2 M # n # M,
the ends of finite-length input sequenc-
3 4 3 4 3 4
y n 5 h m x n 2 m
i 5 0, 1, c, N.
a
m52M
es. The output at the next sample is ob-
n1M
tained by shifting the analysis interval This matrix is called the design matrix
3 4 3 4
5 h n 2 m x m . (4)
a
to the right by one sample, redefining for the polynomial approximation prob-
m5n2M
the origin to be the position of the mid- lem [10]. The transpose of A is
5 6
dle sample of the new block of 2M 1 1 The values marked with 3 in Figure 1 AT 5 ai,n and the product matrix
1 2 1 2
samples, and repeating the polynomial are the shifted impulse responses B 5 ATA is an N 1 1 3 N 1 1 sym-
3 4 3 4
fitting and evaluation at the central lo- h 0 2 m and h 10 2 m that could be metric matrix with elements
cation. This can be repeated at each used to compute the output samples
M M
i1
k
sample of the input, each time labeled with ~, thus replacing the poly- bi, k 5
aai, nan, k 5 an 5bk, i,
n5 2
2 M n5 M
producing a new polynomial and a new nomial fitting process at each sample
3 4
value of the output sequence y n . with a single evaluation of (4). for i 5 0, 1, c, N and k 5 0, 1, c, N,
Another example is shown on the right To show that we can find a single finite- which we see are the coefficients for the
in Figure 1 where the center of the in- duration impulse response that is equiva- set of equations in (6). Furthermore, if
terval is shifted to sample n 5 10 and lent to least-squares polynomial smoothing we define the vector of input samples as
the new polynomial fit to the samples for all shifts of the 2M 1 1-sample interval,
3 3 4 3 4
x 5 x 2 M , c, x 2 1 ,
8 # n # 12 is shown again by the solid we must first determine the optimal coeffi-
curve and the output at n 5 10 is the cients of the polynomial in (1) by differenti-
3 4 3 4 3 4 4T
x 0 , x 1 , c, x M
value of the new polynomial evaluated ating EN in (2) with respect to each of the
3 4T
at the center location. N 1 1 unknown coefficients and setting and define a 5 a0, a1, c, aN as the
The original paper by Savitzky and the corresponding derivative equal to zero. vector of polynomial coefficients, then it
Golay [3] showed that at each position, This yields, for i 5 0, 1, c, N, follows that the equations in (6) can be
the smoothed output value obtained by represented in matrix form as
' EN M N
sampling the fitted polynomial is identi-
3 4
5 2ni a ak nk 2 x n b 5 0,
a a
' ai n52M k50
cal to a fixed linear combination of the Ba 5 ATAa 5 ATx.
(5)
local set of input samples; i.e., the set of
2M 1 1 input samples within the which, by interchanging the order of the Therefore, the solution for the polynomi-
approximation interval are effectively summations, becomes the set of N 1 1 al coefficients can be written as
combined by a fixed set of weighting equations in N 1 1 unknowns
21
1 2
coefficients that can be computed once a 5 ATA ATx 5 Hx.
N M
for a given polynomial order N and
a ni1kbak
a a
approximation interval of length Now recall that the output for the
k50 n52M
2M 1 1. That is, the output samples can group of samples centered on n 5 0 is
M
3 4
be computed by a discrete convolution y 0 5 a0; i.e., we only need to obtain
3 4
5 nix n i 5 0, 1, c, N. (6)
a
of the form n52M the coefficient a0. Furthermore, we see
IEEE SIGNAL PROCESSING MAGAZINE [112] JULY 2011
that we only need the 0th row of the
1 2 1 2 S-G Impulse Response: N = 6, M = 16
N 1 1 3 2M 1 1 m a t r i x H 5
21
1 2
ATA AT, which by the definition of
0.15
matrix multiplication gives a0 as a linear
0.1
1 2
combination of the 2M 1 1 elements
0.05
1 2
of the 2M 1 1 3 1 column vector x. 0
The important observation is that the -0.05
-20 -15 -10 -5 0 5 10 15 20
matrix H depends only on N and M and
Sample Time n
is independent of the input samples.
Thus, the same weighting coefficients
[FIG2] Impulse response of an S-G filter with M 5 16 and N 5 6. The dashed curve is
will be obtained at each group of 2M 1 1
1 2
the polynomial | n evaluated on a dense grid.
p
input samples, and so we can think of
N
least-squares smoothing as a shift-invari- approximation of the impulse input as |,
a
|1n2 5 |
p ak nk 2M # n # M. (8)
a
ant discrete convolution process. which is given by
k50
One approach to finding the impulse
|
21
1 2
response of the equivalent linear time- a 5 ATA ATd, Therefore, the impulse response of the
invariant (LTI) system is to compute the S-G filter is
3 4T
matrix H. Then, by the definition of where d 5 0, 0, c, 0, 1, 0, c, 0, 0
|1n2.
3 4
h 2n 5 h0, n 5 p
1 2
matrix multiplication, the output will be is a 2M 1 1 3 1 column vector
1 2 3 4
impulse and AT is the N 1 1 3 As before, this equation gives h 2 n
M
3 4 3 4 1 2
y 0 5 a0 5 h0, mx m , 2M 1 1 matrix since the impulse response is flipped
a
m52M
around n 5 0 in evaluating discrete
AT 5
where hi, n denotes the elements of the convolution. Henceforth, we shall refer
c c
1 20 1 20
2M 21 1 10 M0
1 2 1 2 1 2
N 1 1 3 2M 1 1 matrix H and h0, m to | n as the impulse response design
p
c c
1 21 1 21
2M 21 0 11 M1
is an element of the 0th row. Therefore, polynomial. As we will discuss in a
c c
1 22 1 22
E 2M 21 0 12 M2U. (7)
comparing this equation to the second later section, (8) is the basis for a sim-
( ( ( ( ( ( (
term of (4) with n 5 0, we observe that ple method for the design of S-G filters
c c
1 2N 1 2N
2M 21 0 1N MN
using the polynomial fitting functions
3 4
h 2 5h0, m 2M #m #M. in MATLAB.
m
Then for the impulse input d, it follows that
3 4 1 2
Note that this equation gives h 2m ATd is the N 1 1 3 1 column vector PROPERTIES OF S-G FILTERS
since, as shown in (4), the impulse Figure 2 shows the impulse response of
3 4T
response is flipped with respect to the ATd 5 1, 0, c, 0 . an S-G filter with N 5 6 and M 5 16.
input in evaluating discrete convolu- Although this is a specific example, its
tion. Efficient matrix inversion tech- This means that the symmetric matrix time-domain properties are representa-
1 221
niques can be employed to compute ATA must have the form tive of the entire class of symmetric S-G
only this first row rather than the entire filters.
| | |
matrix H [10]. a0 a1 c aN Figure 3 shows the frequency re-
| c
Another approach is to note that since sponse of several S-G filters designed by
a1 . .
1 221
ATA 5 e" Ä„ ,
the same weighting coefficients are MATLAB statements given in the section
( ( ( (
| c
obtained irrespective of the signal vector,  Design of S-G Filters. The impulse
aN . .
1 2
we can set x equal to a unit impulse cen- response lengths are all 2M 1 1 5
#
tered in the interval 2M # n # M, and where the matrix entries denoted " do 2 16 1 1 5 33 with implicit polynomial
|
solve for all the coefficients of the corre- not enter into the computation of a. orders of N 5 0, 2, 4, 6, 12. Figures 2 and
sponding polynomial approximation. Note Now, since AT is as given in (7), it fol- 3 illustrate properties shared by all S-G
that these polynomial coefficients, lows from the definition of matrix mul- filters. These properties, which result
denoted as |, will generally not be equal tiplication that the 0th row of the from the structures of the matrices B and
a
1 221
to those of any of the local approxima- matrix H 5 ATA AT is H, are summarized below:
tions that are implicitly generated for P1 The odd-indexed coefficients
c
3 4
h0, , h0, , c, h0, 0, , h0, M
2 M 2 M11
each group of 2M 1 1 input samples. of the impulse response design poly-
Then, the impulse response can be c c |1M24,
nomial are all zero so that we can
3|1 2 1 2 1 2
5 p 2M , | 2M11 , ,| 0 , , p
p p
1 2
obtained by evaluating the corresponding express | n as
p
|1n2 is the polynomial fit to the
polynomial at locations 2M # n # M. where p
N/2
: ;
To show that this statement is true, unit impulse evaluated at the integers
|1n2 5 |
p a2k n2k, (9)
a
we denote the coefficient vector for 2M # n # M, k50
IEEE SIGNAL PROCESSING MAGAZINE [113] JULY 2011
lecture NOTES continued
[ ]
and M, the next section gives an
Frequency Response of S-G Filters (M = 16) approximate empirical relation for
20
fc as a function of both N and M.
N = 0
P7 The S-G filters have mediocre
10 N = 2
N = 4 attenuation characteristics in their
N = 6
stopband regions (except at the fre-
0
N = 12
quencies corresponding to zeros on the
-10
unit circle). Defining the stopband as
the frequency range from the first zero
-20
up to p radians, we see from Figure 3
1 2
that for the MA filter N 5 0 or 1 , the
-30
minimum attenuation in the stopband
-40
(amplitude of first peak after the first
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
zero) is approximately 13 dB. For
Normalized Frequency /
N $ 2, the minimum attenuation in
the stopband is approximately 11 dB.
[FIG3] Frequency response of S-G filters for M 5 16 and various polynomial orders.
Figure 3 also shows that the peak stop-
band gain tends to increase with
where : ; means rounding down. (9) and the normal equations (6) that increasing N for fixed M and that the
1 2
Among other things, this means that H ejv |v50 5 1 and frequency response decreases slightly
the S-G filters for N and N 1 1 are in gain as frequency increases above
M
1 2
drH ejv
identical for N an even integer. 1 2r a 3 4 the nominal cutoff frequency.
` 5 2j nrh n 5 0,
dvr
v50 n52M
P2 Moving average (MA) filter-
(10)
ing defined as DESIGN OF S-G FILTERS
for r 5 1, 2, c, N. Furthermore, it Recall from (8) that the impulse
n1M
1 can be shown using the product rule response of an S-G filter can be comput-
3 4 3 4
y n 5 x m
a
2M 1 1 from differential calculus and ed as samples of the Nth degree polyno-
m5n2M
Parseval s theorem that (10) guaran- mial fit to the unit impulse sequence.
is identical to S-G smoothing with tees that the first N moments of the This method of computing the S-G
3 4
polynomials of order N 5 0 (con- input signal x n are preserved in the filters is easily implemented using
3 4
stant) and M 5 1 (straight line). output y n ; i.e., MATLAB s polynomial functions as in the
P3 The impulse response is following MATLAB statements:
|1n2
3 4
symmetric since h 2n 5 p 5
3 4 3 4
Ä…
|12n2 5 h3n4. Therefore, the fre- n52`nry n 5 Ä… nrx n
p n52` a=polyfit(-ML:MR,...
r 5 1, 2, c, N. (11)
quency response is purely real. [zeros(1,ML),1,zeros
(The shifted impulse response (1,MR)],N);
3 4
h n 2 M is causal and the corre- P6 The nominal normalized h=fliplr( polyval(a,-ML:MR))
sponding frequency response has cutoff (3 dB down) frequency,
linear phase corresponding to the fc 5vc/p, depends on both the The MATLAB function polyfit()
time delay of M samples.) S-G fil- implicit polynomial order N and computes the coefficients of the impulse
ters are type I FIR lowpass filters the length of the impulse response, response design polynomial and poly-
2
[11] with nominal passband gain of (2M 1 1 . If M is fixed as in Figure val() evaluates the polynomial at a
unity. 3, the passband of the filter gets discrete set of points. Note that these
P4 The zeros of the system func- wider approximately in proportion statements can be used to compute non-
1 2
tion H z of an S-G filter are either to N. Although not illustrated in symmetric S-G filters by setting ML?MR.
on the unit circle of the z-plane or Figure 3, the cutoff frequency also The MATLAB Signal Processing Toolbox
they occur in complex conjugate depends inversely on M. S-G filters has a function sgolayfilt() for
reciprocal groups [11]. The unit circle are often compared against a MA fil- designing and implementing both sym-
zeros are, of course, responsible for ter with the same impulse response metric and nonsymmetric S-G filters.
the sharp dips (high attenuation) in length [10]. Figure 3 shows that There are some important constraints
the stopband of the frequency this is somewhat unfair since a in the use of polynomial fitting in general.
responses in Figure 3. shorter MA filter could have rough- Specifically, the number of data points (in
P5 S-G filters have very flat fre- ly the same cutoff frequency as a this case 2M 1 1) must be strictly greater
quency response in their passbands longer S-G filter with higher value than the number of undetermined
since it can be shown using (8) and of N. To clarify this interaction of N coefficients N 1 1 to achieve smoothing
IEEE SIGNAL PROCESSING MAGAZINE [114] JULY 2011
Log Magnitude (dB)
by the S-G process. Furthermore, if the
order of the polynomial, N, is too large, 1
the approximation problem is badly con-
0.9
17
16
15
14
ditioned and the solution will be of no 13
12
0.8 11
10
9
8
7
value. (The function polyfit() issues
6
18
0.7
5
4 19
an alert when the approximation problem
0.6 3 20
is ill conditioned.) Although these factors
2
0.5
are significant limitations, a wide range of
M = 25
0.4
frequency-domain characteristics can be
0.3
achieved nevertheless by choosing M and
M = 50
0.2
N appropriately.
To quantify the frequency-domain 0.1 M = 100
M = 200
behavior of S-G filters, impulse responses
0
0 5 10 15 20 25 30 35 40
were computed for various values of M
Polynomial Order, N
and N within the constraints mentioned
above, and the corresponding fre-
[FIG4] Relationship between 3 dB cutoff frequency fc 5 vc /p, polynomial length N,
quency responses were computed for
and impulse response half-length M.
0 #v#p. The passband of the filter
was defined by the frequency where
1 2
20 log10|H ejv | is  3 dB down from the The values of fc 5vc/p predicted by of M and N. The formula does not fit
value of 0 dB, the gain of the filter at this equation are marked with a red as well for values of M less than 25.
+
v5 0. The results for measurements on and connected by a red line. Figure 4 However, the dependence of fc on N is
filters with M 5 2, 3, c, 20 and shows that this simple formula fits the still linear except for small M. For
M 5 25, 50, 100, 200 for even polyno- measurements quite well even for the 10 # M , 25 and N suitably restrict-
mial orders N are displayed in Figure 4. case M 5 25 where the measurements ed, a formula similar to (12) with 4.6
The points marked with * and connected deviate only slightly from the straight replaced by 2 gives more accurate pre-
by a blue line are the measured cutoff fre- line over the entire range of N. The dictions. While a more complicated
quencies for a fixed value of M. The values relative error in predicting the mea- functional form based on more mea-
for commonly used short filters sured cutoff frequency is less than 4% surements could provide more accu-
1
2 # M # 6 and N , 2M) are given over the range M 5 25, 50, 100, 200 rate predictions over a wider range of
more precisely in Table 1. These are the and N 5 4, c, 32, and the relative M and N, (12) should be adequate for
cutoff frequencies for all possible sym- error is within 8% for the cases most applications of S-G filters where
metric S-G filters for impulse response M 5 25,50, 100, 200 and N 5 2. As can precise specification of the cutoff fre-
1 2
lengths 5 # 2M 1 1 # 13. Observe be seen, large values of M and small N quency is not required.
that the cutoff frequencies that are lead to extremely narrow passbands, Figure 4 points out another feature
achievable range from 0.165 to 0.681, but which would be of limited usefulness of S-G filters that is often overlooked. A
only a discrete set of values is possible. except when the signal components given desired cutoff frequency fc can be
In all cases in Figure 4, fc varies are greatly oversampled. Even though realized by different combinations of N
almost linearly with N when N V 2M the function polyfit() gave an ill- and M. For example, we can achieve a
with the slope being dependent conditioned warning for the larger val- cutoff frequency fc 5vc/p < 0.4 using
inversely on M, but the curves for ues of N, the resulting filters remained the (N, M) pairs (4,4), (9,10), (14,16),
M , 25 tend to deviate from a straight acceptable for values of N up to about and (22,19). These filters differ in the
line as N increases toward 2M. 40. The formula in (12) becomes sharpness of their transition from flat
However, when M is large, as in the four increasingly accurate for larger values passband to stopband, which just
cases M 5 25, 50, 100, 200, the linear reflects a familiar property of FIR dis-
region of the curve coincides with the crete-time lowpass filters that the
range of usable values of N, so a nearly [TABLE 1] NORMALIZED 3 dB widths of transition regions are typi-
CUTOFF FREQUENCIES fc 5vc /p
linear relation holds over a wide range cally inversely proportional to the
AS A FUNCTION OF M AND N.
of N. A reasonably accurate approxima- length of the impulse response.
POLYNOMIAL ORDER, N
tion to this behavior for the indicated
M 2 4 6 8 10
range of parameters is the equation DISCUSSION
2 0.475    
While there is value in knowing that a
3 0.319 0.561   
N 1 1
single S-G impulse response implicitly
4 0.243 0.406 0.615  
M $ 25 and N , M.
fc 5
5 0.197 0.323 0.467 0.653 
3.2M 2 4.6
achieves local polynomial fitting for
6 0.165 0.269 0.382 0.512 0.681
(12)
every output sample, in many
IEEE SIGNAL PROCESSING MAGAZINE [115] JULY 2011
c
c
3 dB Cutoff,
f
=
/
lecture NOTES continued
[ ]
the transition region and the location of
the first zero of the frequency response
 3 dB Least-Squares Filter (M = 16, N = 6)
0
were approximately in the same loca-
Parks-McClellan Filter (L = 2M + 1 = 33)
tion as those of the corresponding S-G
-20
filter. The measured 3 dB cutoff fre-
-40
quency of the S-G filter was fc 5 0.143
(the formula of (12) predicts fc 5 0.15).
-60
A very flat passband is achieved with the
P-M design algorithm by imposing a
-80
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
10:1 weighting ratio between the pass-
Normalized Frequency /
band equiripple approximation error
(a)
and the stopband approximation error.
Detail of Passbands
Larger ratios will make the passband
1
0 even flatter. In the case of the S-G filter,
-1
the gain at the first local maximum
-2
beyond the first zero of the frequency
-3
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
response is 211.73 dB, while the equi-
Normalized Frequency /
ripple maxima of the P-M filter have
(b)
gains of 219.9 dB. The lower part of the
plot shows that the passband gain of the
[FIG5] Comparison of an S-G filter (N 5 6 and M 5 16) with an equal-length
P-M filter has small ripple about 0 dB,
equiripple filter designed with the PM algorithm. (a) Entire frequency response and
(b) the passband region. and the flat region is in fact wider than
that of the S-G filter. However, due to
applications, signals are not character- that the signal is  smooth enough to the tendency of S-G frequency
ized in terms of their ability to be mod- be represented by a polynomial of  high responses to fall off at high frequencies,
eled by polynomials but rather in terms enough degree. However, S-G filters are the S-G filter has slightly lower peak
of their frequency spectra. Thus, we have often used in applications where a direct stopband gain than the P-M filter after
focused in this article on the frequency- frequency-domain specification is more about v/p50.5.
domain properties of the S-G filters. precise or more easily related to models Given the close similarity of the two
S-G filters are often preferred (even for signal production. In such cases, the frequency responses in Figure 5, it is clear
revered in some circles) because, when empirical relationship in (12) or the plot that for the case of a signal confined to the
they are appropriately designed to match of Figure 4 may be useful. Even in the band |v| , 0.143p with additive white
the waveform of an oversampled signal case of the sampled Gaussian wavelets noise, the performance of the two systems
corrupted by noise, they tend to preserve that model chemical spectrum lines, the will be nearly identical. Experiments with
the width and height of peaks in the sig- corresponding Fourier transform also longer impulse responses show that P-M
nal waveform. While such performance has Gaussian shape, and it is straight- filters can achieve very flat passbands and
features are often explained in terms forward to determine the frequency-do- greater stopband attenuation than an
of matching fitted polynomial slopes main width as a function of the width of equivalent S-G filter. Also, recall that the
to signal slopes or to the preservation the Gaussian wavelet. cutoff frequencies of the S-G filters are
of signal moments, the reason for this From the frequency-domain point of restricted to a finite set while those of the
behavior is more obvious from the fre- view, the question naturally arises as to P-M filter are not.
quency-domain properties of the filters. whether the main desirable property of
Specifically, S-G filters have extremely the S-G filters (very flat passband) could WHAT WE HAVE LEARNED
flat passbands with modest attenuation be achieved with another design This article has attempted to answer the
in their stopbands. Furthermore, the method, and perhaps with greater atten- question  What is a Savitky-Golay fil-
symmetric S-G filters have zero phase uation in the stopband region. Figure 5 ter? in terms that will be familiar to
so that features of the signal are not shows the frequency response of an S-G the DSP community and readers of
shifted. Thus, if the signal has most of its filter with M 5 16 (impulse response IEEE Signal Processing Magazine. We
energy in the filter passband (implying length L 5 2M 1 1 5 33) and N 5 6. reviewed the definition and properties
significant over-sampling), the signal Also shown is the frequency response of of S-G filters and showed how they can
components are undistorted while some a length L 5 33 filter designed by the be designed easily using polynomial
high-frequency noise is reduced but not Parks-McClellan (P-M) algorithm. In approximation of an impulse sequence.
completely eliminated. Of course, as- this example, the passband and stop- In contrast to most discussions of S-G
suming that the signal is lowpass and band cutoff frequencies of the P-M filter filters, we focused on the frequency-
oversampled is equivalent to assuming were adjusted by trial and error so that domain properties, and offered an
IEEE SIGNAL PROCESSING MAGAZINE [116] JULY 2011
Log Magnitude (dB)
Log Magnitude (dB)
approximate formula for the 3-dB cutoff coauthor of several DSP textbooks includ- [5] M. Sühling, M. Arigovindan, P. Hunziker, and M.
Unser,  Multiresolution moment filters: Theory and
frequency as a function of polynomial ing Discrete-Time Signal Processing
applications, IEEE Trans. Image Processing, vol.
13, no. 4, pp. 484 495, Apr. 2004.
order N and impulse response half- (with Oppenheim), Signal Processing
[6] M. U. A. Bromba and H. Ziegler,  Applica-
length M. Engineers with a frequency- First (with McClellan and Yoder), and
tion hints for Savitzky-Golay smoothing filters,
Anal. Chem., vol. 53, no. 11, pp. 1583 1586,
domain mindset (like the author) may Theory and Applications of Digital
Sept. 1981.
find this useful if they choose to use Speech Processing (with Rabiner).
[7] R. W. Hamming, Digital Filters, 3rd ed.
S-G filters in their application.
Englewood Cliffs, NJ: Prentice-Hall, 1989.
REFERENCES
[8] S. J. Orfanidis. (1995 2009). Introduction to
[1] K. Pandia, S. Revindran, R. Cole, G. Kovacs,
signal processing [Online]. Available: www.ece.
and L. Giaovangrandi,  Motion artifact cancella-
AUTHOR
rutgers.edu/~orfanidis/intro2sp
tion to obtain heart sounds from a single chest-
Ronald W. Schafer ([e-mail])
worn accelerometer, in Proc. ICASSP-2010, 2010, [9] P.-O. Persson and G. Strang,  Smoothing by
pp. 590 593. Savitzky-Golay and Legendre filters, IMA Vol. Math.
is an HP Fellow in the Mobile and
Systems Theory Biol., Comm., Comp., and Finance,
[2] R. W. Schafer,  On the frequency-domain
vol. 134, pp. 301 316, 2003.
Immersive Experience Lab at HP Labs,
properties of Savitzky-Golay filter, in Proc. 2011
DSP/SPE Workshop, Sedona, AZ, Jan. 2011, pp.
[10] W. H. Press, S. A . Teukolsk y, W. T.
Palo Alto, California, where he is involved
54 59.
Vertterling, and B. P. Flannery, Numerical
in research on acoustic and audio signal
Recipes, 3rd ed. Cambridge, U. K.: Cambridge
[3] A. Savitzky and M. J. E. Golay,  Soothing and
Univ. Press, 2007.
differentiation of data by simplified least squares
processing. From 1974 to 2004, he was
procedures, Anal. Chem., vol. 36, pp. 1627 1639,
[11] A. V. Oppenheim and R. W. Schafer, Discrete-
John and Marilu McCarty Professor of the 1964.
Time Signal Processing, 3rd ed. Upper Saddle
[4] J. Riordon, E. Zubritsky, and A. Newman,  Top River, NJ: Pearson, 2010.
School of Electrical and Computer
10 articles, Anal. Chem., vol. 72, no. 9, pp. 324A
Engineering at Georgia Tech. He is the 329A, May 2000. [SP]
Kivanc Kose and A. Enis Cetin
Low-Pass Filtering of Irregularly Sampled
Signals Using a Set Theoretic Framework
n this article, the goal is to show samples as a part of the filtering pro- PROBLEM STATEMENT
1 2
that it is possible to filter non- cess. Since there are specifications in Let us assume that samples xc ti ,
u n i f o r m l y s a m p l e d s i g n a l s both time and frequency domains, it i 5 0, 1, 2, c, L 2 1, of a continuous
1 2
according to specs defined in the is possible to iterate between time and time-domain signal xc t are available.
IFourier domain. In many practi- frequency domains using the fast These samples may not be on an uni-
cal applications, it is necessary to fil- Fourier transform (FFT) while impos- form sampling grid. Let us define
3 4 1 2
ter irregularly sampled data including ing the constraints in each domain. xd n 5 xc nTs as the uniformly sam-
seismic signal processing, synthetic pled version of this signal. We assume
aperture radar (SAR) imaging sys- RELEVANCE that the sampling period Ts is sufficient-
tems, three-dimensional (3-D) mesh- The ideas presented here can be used to ly small (below the Nyquist period) for
1 2
es, and digital terrain models [1], [2]. develop filtering algorithms for irregu- the signal xc t . In a typical discrete-
3 4
In almost all of these practical prob- larly sampled one or higher dimension- time filtering problem, we have xd n or
lems, it is possible to define the al data. It can be used as a teaching its noisy version, and we apply a dis-
desired filtering solution in a set the- material in advanced undergraduate crete-time low-pass filter to the uni-
3 4
oretic framework. This lecture note and graduate discrete-time signal pro- formly sampled signal xd n . However,
3 4
presents a new method for filtering cessing, optimization as well as applied xd n is not available in this problem.
1 2
irregularly sampled data by defining mathematics courses. Only nonuniformly sampled data xc ti ,
stopband tolerance regions in the i 50, 1, 2, c, L 21 are available in this
Fourier domain and time-domain PREREQUISITES problem.
upper and lower bounds on the signal The prerequisites for understanding this
article s material are linear algebra, dis- GOAL
crete-time signal processing, and basic Our goal is to low-pass filter the non-
Digital Object Identifier 10.1109/MSP.2011.941098
1 2
Date of publication: 15 June 2011 optimization theory. uniformly sampled data xc ti according
1053-5888/11/$26.00©2011IEEE IEEE SIGNAL PROCESSING MAGAZINE [117] JULY 2011