Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I do have noisy data and want to smooth them by a Savitzky-Golay filter because I want to keep the magnitude of the signal.

a) Is there a ready-to-use Filter available for that?

b) what are appropriate values for m (the half width) and for the coefficients for 3000-4000 data points?

Example data: data

share|improve this question
3  
There is one here library.wolfram.com/infocenter/MathSource/789 –  Nasser Nov 19 '13 at 22:55
    
Smoothing data depends on lots of things - for example if you are recording the data in cm^-1, but you actually need the spectrum in nm. Describe the specific applications and also the sampling rate as there are parameters that depend on it. Also have you considered using other tools to smooth your data ? –  Sektor Nov 19 '13 at 22:55
    
@Nasser: found that link, but the m-file can not be downloaded. –  Shukoff Nov 19 '13 at 23:01
2  
the m-file can not be downloaded It worked for me, just downloaded it the m file and the .nb file. The trick is to right-click->Save Link As... and not to click on it. –  Nasser Nov 19 '13 at 23:11
5  
Here is a link to some code posted on mathgroup that I have used before and found ok: forums.wolfram.com/mathgroup/archive/2012/Feb/msg00036.html –  Mike Honeychurch Nov 19 '13 at 23:46

3 Answers 3

up vote 15 down vote accepted

The following code will filter noisy data…

SGKernel[left_?NonNegative, right_?NonNegative, degree_?NonNegative, derivative_? NonNegative] :=
Module[{i, j, k, l, matrix, vector},
    matrix = Table[  (* matrix is symmetric *)
        l = i + j;
        If[l == 0,
            left + right + 1,
        (*Else*)
            Sum[k^l, {k, -left, right}]
        ],
        {i, 0, degree},
        {j, 0, degree}
    ];
    vector = LinearSolve[
        matrix, 
        MapAt[1&, Table[0, {degree+1}], derivative+1]
    ];
    (* vector = Inverse[matrix][[derivative + 1]]; *)
    Table[
        vector.Table[If[i == 0, 1, k^i], {i, 0, degree}],
        {k, -left, right}
    ]
] /; derivative <= degree <= left+right


SGSmooth[list_?VectorQ, window_, degree_, derivative_:0]:=
    Module[{pairs},
        pairs = MapThread[List, {Range[Length[list]], list}];
        Map[Last, SGSmooth[pairs, window, degree, derivative]]
    ] 


SGSmooth[list_?MatrixQ, window_, degree_, derivative_:0]:=
    Module[{kernel, list1, list2, margin, space, smoothData},

        (* determine a symmetric margin at the ends of the raw dataset.
        The window width is split in half to make a symmetric window 
        around a data point of interest *)
        margin = Floor[window/2];

        (* take only the 1st column of data in the list to be smoothed (the
        independant Values) and extract the data from the list by removing
        half the window width 'i.e., margin' from the ends of the list *)
        list1 = Take[Map[First, list], {margin + 1, Length[list] - margin}];

        (* take only the 2nd column of data in the list to be smoothed 
        (the dependent Values) and Map them into list2 *)
        list2 = Map[Last, list];

        (* get the kernel coefficients for the left and right margins, the 
           degree, and the requested derivative *)
        kernel = SGKernel[margin, margin, degree, derivative];

        (* correlation of the kernel with the list of dependent values *)
        list2 = ListCorrelate[kernel, list2];

        (* Data _should_ be equally spaced, but... calculate spacing anyway by getting
           the minimum of all the differences in the truncated list1, remove the first 
           and last points of list1 *)
        space = Min[Drop[list1, 1] - Drop[list1, -1]];

        (* condition the dependant values based on spacing and the derivative *)
        list2 = list2*(derivative!/space^derivative);

        (* recombine the correlated (x-y) data pairs (that is list1 and list2), 
        put these values back together again to form the smooth data list *)
        smoothData=Transpose[{list1, list2}]


    ] /; derivative <= degree <= 2*Floor[window/2] && $VersionNumber >= 4.0

I did not apply this to your data, but you can do that later. This example is applied to noisy random data.

Using a noisy sine data function…

dataFunction[x_] := Sin[x] + Random[Real, {-\[Pi], \[Pi]}];

Build a table of noisy tabular data from 0 to 2[Pi]…

dataTable = Table[{x, dataFunction[x]}, {x, 0, 2 \[Pi], .01}];

Animate the smoothing operations. Notice the smoothed dataset shrinks with increasing 'window width'. This is an artifact of the ListCorrelate function used in the SGSmooth function. ListCorrelate uses an end buffered dataset.

NOTE: The red line is the filtered data set…

 Manipulate[
 If[showRawData,
  Show[
   ListPlot[dataTable, PlotRange -> {{0, 2 \[Pi]}, {-5.0, 5.0}}],
   ListPlot[
    {
     SGSmooth[dataTable, win, order]
     },
    PlotRange -> {{0, 2 \[Pi]}, {-5.0, 5.0}},
    PlotStyle -> {{Red, Thick}, {Green, Thick}},
    Joined -> True]
   ], (* ELSE just plot smooth data *)
  ListPlot[
   {
    SGSmooth[dataTable, win, order]
    },
   PlotRange -> {{0, 2 \[Pi]}, {-5.0, 5.0}},
   PlotStyle -> {{Red, Thick}, {Green, Thick}},
   Joined -> True]
  ],

 {{win, 100, "window width"}, 2, 300, 1, 
  Appearance -> "Labeled"}, {{order, 1, "order of polynomial"}, 1, 9, 
  1, Appearance -> "Labeled"},
 {{showRawData, True, "Raw Data"}, {True, False}},
 SaveDefinitions -> True
 ]

This will create a Manipulate similar to the following:

enter image description here

Hope this helps!

share|improve this answer
    
Great!. Thanks for sharing the filter. The last Definition of SGSmooth has a ; too much preventing any output. –  Markus Roellig Nov 20 '13 at 19:47
    
@MarkusRoellig thank you finding this little typo. I fixed the above code to plot the smoothed data. I had originally used SGSmooth[] with additional embedded code and did not remove the ';'. My apologies... –  Joseph Nov 21 '13 at 1:43
    
@Joseph: Thank you, that answer is great! Just a question: It is described in literature that the calculation of the derivatives can be directly implemented in the SG-filter (see: pubs.acs.org/doi/abs/10.1021/ac50031a048). I guess this would be an expansion of your solution, right? If so, I open another question and let you know? –  Shukoff Nov 25 '13 at 15:06
1  
@Shukoff: You can use the existing code as is to calculate the data derivative. That is the beauty of the SG filter. This implementation is very similar to the Numerical Recipes implemetation (wire.tu-bs.de/OLDWEB/mameyer/cmr/savgol.pdf). Just use: SGSmooth[dataTable, win, order, 1] for 1st derivative, SGSmooth[dataTable, win, order, 2] for 2nd derivative, etc... –  Joseph Nov 25 '13 at 15:14
    
@Joseph: your implementation is superb! Thanks! –  Shukoff Nov 25 '13 at 15:16

Several years ago in related MathGroups thread Virgil P. Stokes suggested:

A few years back I wrote a Mathematica notebook that shows how one can obtain the SG smoother from Gram polynomials. The code is not very elegant; but, it is a rather general implementation that should be easy to understand. Contact me if you are interested and I will be glad to forward the notebook to you.

I contacted him and received the notebook. I find his implementation of the Savitzky-Golay filter quite stable and working pretty well. Here I publish it with his permission:

Clear[m, i]; (* m, i are global variables !! *)
Clear[GramPolys, LSCoeffs, SGSmooth];

GramPolys[mm_, nmax_] := 
  Module[{k, m = mm},  (* equations (1a), (1b) *)
      (* 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}]    
   ]; 

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}]
   ];

SGSmooth[cc_, data_] := Module[{m, y, datal, 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}]]
   ];
Usage: SGOutput = SGSmooth[LSCoeffs[m,n,d], data]
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; ...)
   data = list of uniformly sampled (spaced) data values to be smoothed (length(data) >=2m+1)  
Outputs:
   SGOutput =  list of smoothed data values
share|improve this answer
    
I really should try looking for that old notebook I wrote for these coefficients; I remember using an algorithm that did not need to generate the Gram polynomials themselves, and worked with the recurrence coefficients. –  Guess who it is. Jun 8 at 13:53
1  
@Guesswhoitis. Here you can download the original Notebook sent me by Virgil. –  Alexey Popkov Jun 8 at 13:57

I'm just posting this to record for posterity something I posted in the chatroom a not-so-long time ago. As I noted there, the following routine will only do smoothing; I had a more general routine for generating the differentiation coefficients, but I still have not been able to find it. As with Virgil's method (the one in Alexey's answer), this is based on Gorry's procedure (though I have traced the spirit of the algorithm going as far back as Hildebrand's book):

GramP[k_Integer, m_Integer, t_Integer] :=
     (-1)^k HypergeometricPFQ[{-k, 1 + k, -m - t}, {1, -2 m}, 1]

SavitzkyGolay[n_Integer, m_Integer, t_Integer] := 
        Table[Sum[(Binomial[2 m, k]/Binomial[2 m + k + 1, k + 1])
                  GramP[k, m, i] GramP[k, m, t] (1 + k/(k + 1)), {k, 0, n}, 
                  Method -> "Procedural"], {i, -m, m}]
SavitzkyGolay[n_Integer, m_Integer] := Table[SavitzkyGolay[n, m, t], {t, -m, m}]

Usage is pretty straightforward: n is the order of the polynomial smoothing; 2 m + 1 is the window size, and t tells how much to shift the window.


I managed to finally recover the routine I once wrote through the kind assistance of a friend. To share my joy, I now release this to you:

Options[SavitzkyGolay] = {Derivative -> 0, WorkingPrecision -> Infinity};

SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, t_Integer,
              OptionsPattern[]] /; 1 < n < 2 m + 1 && -m <= t <= m := 
Module[{o = OptionValue[Derivative], c, s, h, p, q, u, v, w},
       u = UnitVector[o + 1, 1]; v = ConstantArray[0, o + 1]; 
       c = 1/(2 m + 1); s = Join[{Boole[o == 0] c},
       Table[h = 0;
             {p, q} = {2 (2 k - 1), (k - 1) (2 m + k)}/(k (2 m - k + 1));
             Do[w = u[[j]]; (* evaluate Gram polynomial and derivatives *)
                u[[j]] = p (t w + (j - 1) h) - q v[[j]];
                v[[j]] = h = w,
                {j, Min[k, o] + 1}];
             c *= (2 m - k + 1) (1 + 1/k)/(2 m + k + 1);
             c (1 + k/(k + 1)) u[[o + 1]],
             {k, n}]];
       Table[h = s[[n]] + 2 (2 n - 1) (p = s[[n + 1]]) j/(n (2 m - n + 1));
             Do[q = p; p = h; (* Clenshaw's recurrence *)
                h = s[[k]] + 2 (2 k - 1) p j/(k (2 m - k + 1)) -
                    k (2 m + k + 1) q/((k + 1) (2 m - k)),
                {k, n - 1, 1, -1}];
                N[h, OptionValue[WorkingPrecision]], {j, -m, m}] // 
    Developer`ToPackedArray];

SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___?OptionQ] /; 
  1 < n < 2 m + 1 := 
 Developer`ToPackedArray[Table[SavitzkyGolay[n, m, t, opts], {t, -m, m}]]

As advertised, it uses no matrices, and instead uses the recurrence relation of the Gram polynomial. If need be, the guts of the routine can be embedded within a Compile[].

share|improve this answer
    
Please add an example how it should be applied to a data (as an interesting data I suggest IR spectrum of gaseous water). What is the meaning of the coefficients which SavitzkyGolay returns? –  Alexey Popkov Jul 4 at 6:58
    
@Alexey, I'll add a usage example later. For now: SavitzkyGolay[n, m, Derivative -> d] corresponds to LSCoeffs[m, n, d] in Virgil's implementation in your answer. –  Guess who it is. Jul 4 at 7:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.