General Least-Square Smoothing and Differentiation
V. P. Stokes (vs@it.uu.se), 2009-04-16, 2009-04-12

Introduction

The main purpose of this notebook is to show how exact values for all coefficients in a sliding window (convolution kernel) over all points in sequence of equally spaced data values can be generated for all derivative orders for least-square (LS) smoothing. I have implemented a single function GramPolys that can be used to obtain the convolution kernel under these rather general conditions.These are often referred to as the Savitzky-Golay (SG) smoothers, since Savitzky and Golay [3] were the first to show how these polynomial could be applied to the central points of equally spaced data values.
 
The earliest reference that I have found for the application of the recursive form of the Gram polynomials to LS approximation was by Gorry [1]. Later, these same recursions were given by Jakubowska and Kubiak [2]. Note, that I have derived recursive forms based on the Gram polynomials that also include the derivatives [4] which for convenience I will refer to as SR  (Stokes recursions). The following is taken from [2],

(Subscript[P, 0]^m)(i)=1, (Subscript[P, -1]^m)(i)=0

(Subscript[P, k]^m)(i)=((2(2k-1))/(k(2m-k+1)))i (Subscript[P, k-1]^m)(i)-(((k-1)(2m+k))/(k(2m-k+1))) (Subscript[P, k-2]^m)(i), k=1,2,...,n

i = point of evaluation in the window ( i∈[-m,m] )
k = degree of polynomial
n = maximum degree to be used
2m+1 = size (number of points) in sliding window

The coefficients of the convolution kernel for these complete LS approximation is given by,

(Subscript[c, i])=(∑_(j=-m)^m)(∑_(k=0)^n)(((2k+1)(((2m))^(k))))/((2m+k+1))) (Subscript[P, k]^m)(j) (Subscript[P, k]^m)(i), i =-m,-m+1,...,-1,0,1,...,m-1,m
Note, that this includes the left-end (front) points which uses (Subscript[c, i]), i=-m,-m+1,...-1 and the right-end (trailing) points of the data sequence which uses (Subscript[c, i]), i=1,2,...m. I first  show one way to implement these equations in Mathmatica and then how they can be applied to uniformly sampled data.

Gram Polynomial Function

Let's first generate the Gram polynomials in symbolic form (in terms of m and i ) as the degree (k) varies using these recursions. Note, this low-level function will be used later for the generation of the coefficients in the convolution kernel.

Clear[m, i]; (* m, i are global variables !! *)
Remove[GramPolys];
GramPolys[mm_,nmax_] := Module[{k,m=mm}, 
   (* Define recursive equation for Gram polynomials as a function of m,i for degrees 0,1,...,nmax *)
   p[m, 0, i] = 1; 
   p[m, -1, i] = 0;    
   p[m_, k_, i_] := p[m, k, i] = 2*(2*k - 1)/(k*(2*m - k + 1))*i*p[m, k - 1, i] - 
                                   (k - 1)*(2*m + k)/(k*(2*m -k+ 1))*p[m, k - 2, i];
                                   
   (* Return coefficients for degrees 0,1,...,nmax in a list *) 
   Table[p[mm,k,i]//FullSimplify,{k,0,nmax}]    
];

Test GramPolys Function

Generate the Gram polynomials for degee 0,1,2,...,5

n=5;
polynomial = GramPolys[3,n]
polynomial = GramPolys[m,n]

Let's look at some of these polynomials and their deriviatives wrt to i (used to define derivatives for SG smoothing)

polynomial[[1]]
polynomial[[2]]
polynomial[[3]]
polynomial[[4]]
D[polynomial[[4]],i]
polynomial[[5]]
D[polynomial[[5]],{i,3}]
polynomial[[2]] /. {i -> 1}
polynomial[[3]] /. {i ->-2}
polynomial[[4]]/.{i->-1}
polynomial[[5]]/.{m->11,i->-21}
D[polynomial[[4]],i]
polynomial[[6]]

LS Smoothing Coefficients Examples

LS Coefficients in Terms of i

m=3; (* 2m+1 = 7 point sliding window *)
n = 3; (* degree of polynomials = 3, cubic kernel *)
polynomial=GramPolys[m,n];

clist ={};
Do[(* points in each sliding window *)
sum = 0;
Do[ (* degree loop *)
    num = (2 k+1)FactorialPower[2 m,k];
    den = FactorialPower[2 m+k+1,k+1];
    t1 = polynomial[[k+1]]/.{i->j};
    t2 = polynomial[[k+1]]; 
    sum= sum+(num/den)*t1*t2//FullSimplify;
    (*Print["k,polynomial[[k+1]]: ",k,", ",polynomial[[k+1]]];*)
,{k,0,n}];
   clist = Append[clist,sum];
,{j,-m,m}]; clist

Table of LS Coefficients for Derivatives

Note, derivative order = 0, corresponds to the LS smoother.

dorder = 3; (* Derivative order desired *)
Do[
Print["derivative order = ",dorder];
Do[
  klist = D[clist,{i,dorder}]/.{i->j};
If[j≥0,str = "m= ",str="m="];
  Print[str,j,": ",klist];
,{j,-m,m}];
,{dorder,0,4}];

I now introduce a function to generate all the LS coefficients for any deriviative order (d=0,1,...),

Usage: cc = LSCoeffs[m,n,d)
Inputs: 
   m  =   half-width of smoothing window; i.e., 2m+1 points in smoothing kernel
   n =  degree of LS polynomial (n < 2m+1)
   d  =  order of derivative (d =0, smoother; d = 1, 1st derivative; ...)
Otputs:
  cc = list of coefficients
     cc[[1]],cc[[2]],...,c[[m]],   cc[[m+1]],  cc[[m+2]],...,cc[[2m]],cc[[2m+1]]
    -m,    -m+1,        1,         0,          1,           m-1,       m
     |<---------------- left-end -------------->|<-------- central ------->|<-------------------- right-end --------------------->|
                             1,2,...,m,                            m+1,m+2,...,n-m,                              n-m+1,n-m+2,...,n
     Each cc[[i]] is a list with 2m+1 coefficients in it.

Remmove[LSCoeffs];
LSCoeffs[m_,n_,d_] := Module[{k,j,sum,clist,polynomial,cclist},
polynomial=GramPolys[m,n];
clist ={};
Do[(* points in each sliding window *)
sum = 0;
Do[ (* degree loop *)
    num = (2 k+1)FactorialPower[2 m,k];
    den = FactorialPower[2 m+k+1,k+1];
    t1 = polynomial[[k+1]]/.{i->j};
    t2 = polynomial[[k+1]]; 
    sum= sum+(num/den)*t1*t2//FullSimplify;
    (*Print["k,polynomial[[k+1]]: ",k,", ",polynomial[[k+1]]];*)
,{k,0,n}];
   clist = Append[clist,sum];
,{j,-m,m}];
Table[D[clist,{i,d}]/.{i->j},{j,-m,m}]
];

A quick test of this function.

(*m=3;*)
m=3; n=3;
cc=LSCoeffs[m,n,0]
cc[[m+1]] (* coefficients for central points *)
cc=LSCoeffs[m,n,1] (* First derivative *)
cc[[m+1]]

Finally, we define a  function that uses ListConvolve to "smooth" a data sequence with our LS convolution kernel.

Usage: SGOutput = SGSmooth[cc, data]
Inputs: 
   cc      =   list of LS kernel coefficients (from cc = LSCoeffs[m,n,d) )
   data =   list of uniformly sampled (spaced) data values to be smoothed (length(data) ≥2m+1)  
Otputs:
   SGOutput =  list of smoothed data values

Remove[SGSmooth];
SGSmooth[cc_,data_]:= Module[{m,y,data1,datar,k,kk,n,yy},

n = Length[data];
m = (Length[cc]-1)/2;

(* Left end  --- first 2*m+1 points used *)
datal = Take[data,2*m+1];
(* Smooth first m points (1,2,...,m-1,m) *)
kk = 0;
Table[(kk=kk+1;y[k]=ListConvolve[Reverse[cc[[kk]]],datal][[1]]),{k,-m,-1}];

(* Smooth central points (m+1,m+2,...n-m-1) *)
y[0]=ListConvolve[Reverse[cc[[m+1]]],data];

(* Right end --- last 2*m+1 points used *)
datar = Take[data,{n-(2*m+1)+1,n}];
(* Smooth last m points (n-m,n-m+1,...,n) *)
kk = m+1;
Table[(kk=kk+1;y[k]=ListConvolve[Reverse[cc[[kk]]],datar][[1]]),{k,1,m}];

(* And now we concatenate the front-end, central, and back-end estimated data values *)
yy = Join[Table[y[k],{k,-m,-1}],y[0],Table[y[k],{k,1,m}]]
];

Smoothing Examples

Let's apply our smoother (0-order derivative) to some simulated data. First we get the LS polynomial coefficients for the smoother.

(*m=20; (* Must be declared global *)*)
Clear[m];
m=20; n=4;
cc=LSCoeffs[20,4,0];
cc[[m+1]]

Example-1

Linear function without noise (constant slope)

ndata = 51; (* n ≥ 2*m+1 *)
data = Table[k,{k,1,ndata}]

(*SGSmooth[LSCoeffs[[1]],data] (* smoothed data *)*)
SGSmooth[cc,data]
(*SGSmooth[LSCoeffs[[2]],data] (* 1st derivative of data *)
SGSmooth[LSCoeffs[[3]],data] (* 2nd derivative of data *)*)

Example-2

Suppose the following represents a function that we wish to measure,

Remove[epsp];
epsp[t_] := If[t > 50., -65. -E^(-(t - 50.) / 10.) + E^(-(t - 50.) / 100.),-65.]

We plot it and consider this to be the "true" underlying signal that we would like to extract from a measured sample of this signal.

ndata = 301;
epspPlot = Plot[epsp[t], {t, 0, ndata-1}, PlotRange -> {-66, -63},PlotStyle -> RGBColor[1, 0, 0]]

Add some Gaussian noise to simulate the measured data. Thus, the measured signal is the underlying signal + Gaussian noise (additive noise).

SeedRandom[1234];
measurement =
Table[epsp[t] + Random[NormalDistribution[0, 0.10]],{t, 0, ndata-1}];

Ok, we now have our measured signal (data) --- let's plot it.

measurementPlot = ListPlot[measurement, PlotRange -> {-66, -63},PlotJoined -> True,
                           PlotStyle -> {GrayLevel[0.5]}]

We plot the coefficients of our LS polynomial over the interval {-m,-m+1,...,-1,0,1,...,m-1,m} which will be convolved with the central points of the measured data. Note, the left and right ends of the data also are convolved with non-central LS coefficients (these are not plotted for brevity).

ListPlot[Transpose[{Range[-m, m], cc[[m+1]]}],AxesLabel -> {"n", "c[n]"}]

We now smooth these measured data via convolution of the kernel (LS polynomial coefficients) with the measured data.

smoothOutput=SGSmooth[cc,measurement];

This is what the smoothed measured data looks like (we plot each smoothed point).

smoothOutputPlot = ListPlot[smoothOutput, PlotRange -> {-66, -63}]

And finally we superimpose the actual signal (red), the measured data (gray), and the smoothed data (blue).

Show[measurementPlot, smoothOutputPlot, epspPlot]

It seems clear from this plot that the smoother is able to remove much of the additive noise and thus we can obtain a reasonable approximation to the underlying signal.

Example-3

We consider the following Sin function to be the "true" underlying signal

                   Sin(((2π)/400) t)
                   
 that we would like to extract from a measured sample of this signal.  We plot one cycle of this sinusoid,

ndata = 401;
sigPlot =
Plot[Sin[(2*Pi/400)*t], {t, 0, ndata-1}, PlotRange -> {-1.5, 1.5},PlotStyle -> RGBColor[1, 0, 0]]

Add some Gaussian noise to simulate a signal that was measured. Thus, the measured signal is the underlying signal + Gaussian noise (additive noise).

SeedRandom[234191731]; (* use a seed for reproducibility *)
sdnoise = 0.07; (* standard deviation of noise *)
measurement = Table[Sin[(2*Pi/400)*t] + Random[NormalDistribution[0, sdnoise]],{t, 0, ndata-1}];

Ok, we now have our measured signal --- let's plot it.

measurementPlot = ListPlot[measurement, PlotRange -> {-1.5, 1.5},PlotJoined -> True,
                           PlotStyle -> {GrayLevel[0.5]}]

We also plot the coefficients of our LS polynomial for the central points over the interval {-m,-m+1,...,-1,0,1,...,m-1,m} which will be convolved with the all points except the first m points and the last m points. Note, the left and right ends of the data also are convolved with non-central LS coefficients; but, for brevity these are not plotted.

ListPlot[Transpose[{Range[-m, m], cc[[m+1]]}],AxesLabel -> {"n", "c[n]"}]

Apply this smoother to the measured data via convolution of the kernel (LS polynomial coefficients) with the measured data.

smoothOutput=SGSmooth[cc,measurement];

This is what the smoothed measured data looks like.

smoothOutputPlot = ListPlot[smoothOutput, PlotRange -> {-1.5, 1.5}]

And finally we superimpose the actual signal (red), the measured data (gray), and the smoothed data (blue).

Show[measurementPlot, smoothOutputPlot, sigPlot]

It seems clear from this plot that the smoother is able to remove much of the additive noise and thus we can obtain a reasonable approximation to the underlying signal.

Example-4

Let's look at the estimation of the 1st deriviate of this noisy sinusoid. We first, reduce the degree of the LS polynomial to reduce the high-frequence variations introduced by differentiation.

m=20; n=2;
cc=LSCoeffs[20,4,0];

We plot one cycle of the "true" underlying signal that we would like to extract from the noisy measured sample, which is just

                       ((2π)/400) Cos(((2π)/400)t)
 
 That is, we would like to obtain a good estimate of the true 1st derivative of the underlying signal that contains additive Gaussian noise.

ndata = 401;
amp = 2*Pi/400;
derPlot =
Plot[amp*Cos[(2*Pi/400)*t], {t, 0, ndata-1}, PlotRange -> {-0.04, 0.04},PlotStyle -> RGBColor[1, 0, 0]]

We also plot the coefficients of our LS polynomial for the central points over the interval {-m,-m+1,...,-1,0,1,...,m-1,m} which will be convolved with all points except the first m points and the last m points. Note, the left and right ends of the data also are convolved with non-central LS coefficients; but, for brevity these are not plotted.

cc=LSCoeffs[m,n,1];

ListPlot[Transpose[{Range[-m, m], cc[[m+1]]}], AxesLabel -> {"n", "c[n]"}]

Apply this 1st deriviative smoother to the measured data via convolution of the kernel (LS polynomial coefficients) with the measured data.

dsmoothOutput=SGSmooth[cc,measurement];

This is what the smoothed estimate for the 1st derivative of the noisy measured data looks like.

dsmoothOutputPlot = ListPlot[dsmoothOutput, PlotRange -> {-0.04, 0.04}]

And finally we superimpose the actual 1st derivative of the underlying signal (red) and the smoothed estimate for the 1st derivative of the noisy measured data (blue).

Show[dsmoothOutputPlot, derPlot]

It should be clear that we have obtained a rather good estimate for the 1st derivative our noisy measured data.

References

[1] Gorry, P. A. (1990) General Least-Squares Smoothing and Differentiation by the Convolution (Savitzky-Golay) Method. Analytical Chemistry, 62, no. 6 (March), pp. 570-573.
[2] Jakubowska, M. and Kubiak, W. W. (2004) Adaptive-degree polynomial filter for voltammetric signals. Analytica Chimica Acta, 512,  pp. 241-250.
[3] Savitzky, A. and Golay M.J.E. (1964) Smoothing and Differentialtion of Data by Simplified Least Squares Procedures. Analytical Chemistry, 36, no. 8 (July), pp. 1627-1639.
[4] Stokes, V.P. (1997) Methods for Analysis of 3D Movement of the Trunk During Treadmill Locomotion. Ph.D. Thesis. Karolinska Institute. Stockholm, Sweden.