4 faster Gram polynomial; added notes on built-in SavitzkyGolayMatrix edited Dec 16 '15 at 16:28 J. M. is away♦ 100k1010 gold badges317317 silver badges474474 bronze badges 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 aI had a more general routine for generating the differentiation coefficients, but I still have not been able to find it. (For the more general routine for generating the differentiation coefficients, but I still have not been able to find itless compact version, see below.) 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, 1k + k1, -mt - tm}, {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.Added 12/17/2015 Here is a faster routine for evaluating the Gram polynomial, using some undocumented functionality:GramP[k_Integer, m_Integer, t_Integer] := (-1)^k InternalDCHypergeometricPFQ[k, {-k, k + 1, -m - t}, {1, -2 m}, 1]  I managed to finally recover the general SG routine I once wrote through the kind assistance of a friend. To share my joy, I now release this to you: 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[].Added 12/17/2015 Altho SavitzkyGolayMatrix[] is now built-in in version 10, it is only limited to producing the "central" coefficients, as opposed to the routine SavitzkyGolay[] in this answer that can also generate coefficients for the left and right ends.SavitzkyGolayMatrix[{2}, 3, 1, WorkingPrecision -> ∞] {1/12, -2/3, 0, 2/3, -1/12} SavitzkyGolay[3, 2, 0 (* central *), Derivative -> 1] {1/12, -2/3, 0, 2/3, -1/12}  In general, the result of SavitzkyGolayMatrix[] is built from appropriate outer products of coefficient lists.SavitzkyGolayMatrix[{3, 4}, {2, 3}, WorkingPrecision -> ∞] === Outer[Times, SavitzkyGolay[2, 3, 0], SavitzkyGolay[3, 4, 0]] True  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: 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[]. 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. (For the more general, but less compact version, see below.) 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, k + 1, -t - m}, {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.Added 12/17/2015 Here is a faster routine for evaluating the Gram polynomial, using some undocumented functionality:GramP[k_Integer, m_Integer, t_Integer] := (-1)^k InternalDCHypergeometricPFQ[k, {-k, k + 1, -m - t}, {1, -2 m}, 1]  I managed to finally recover the general SG routine I once wrote through the kind assistance of a friend. To share my joy, I now release this to you: 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[].Added 12/17/2015 Altho SavitzkyGolayMatrix[] is now built-in in version 10, it is only limited to producing the "central" coefficients, as opposed to the routine SavitzkyGolay[] in this answer that can also generate coefficients for the left and right ends.SavitzkyGolayMatrix[{2}, 3, 1, WorkingPrecision -> ∞] {1/12, -2/3, 0, 2/3, -1/12} SavitzkyGolay[3, 2, 0 (* central *), Derivative -> 1] {1/12, -2/3, 0, 2/3, -1/12}  In general, the result of SavitzkyGolayMatrix[] is built from appropriate outer products of coefficient lists.SavitzkyGolayMatrix[{3, 4}, {2, 3}, WorkingPrecision -> ∞] === Outer[Times, SavitzkyGolay[2, 3, 0], SavitzkyGolay[3, 4, 0]] True  3 added 10 characters in body edited Jul 7 '15 at 4:01 J. M. is away♦ 100k1010 gold badges317317 silver badges474474 bronze badges Usage is pretty straightforward: n is the order of the polynomial smoothing; 2 m+1m + 1 is the window size, and t tells how much to shift the window.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}] // DeveloperToPackedArray]; SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___]opts___?OptionQ] /; 1 < n < 2 m + 1 := DeveloperToPackedArray[Table[SavitzkyGolay[n, m, t, opts], {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.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}] // DeveloperToPackedArray]; SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___] /; 1 < n < 2 m + 1 := DeveloperToPackedArray[Table[SavitzkyGolay[n, m, t, opts], {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.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}] // DeveloperToPackedArray]; SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___?OptionQ] /; 1 < n < 2 m + 1 := DeveloperToPackedArray[Table[SavitzkyGolay[n, m, t, opts], {t, -m, m}]]  2 finally! edited Jul 4 '15 at 1:09 J. M. is away♦ 100k1010 gold badges317317 silver badges474474 bronze badges This is a placeholder answer; 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): The complete routine I once had entirely skips having to defineUsage is pretty straightforward: n is the Gram polynomials separately andorder of the polynomial smoothing; 2 m+1 is also ablethe window size, and t tells how much to produceshift the derivative coefficientswindow. I'll post I managed to finally recover the routine as soon as I find that notebookonce 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}] // DeveloperToPackedArray]; SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___] /; 1 < n < 2 m + 1 := DeveloperToPackedArray[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[]. This is a placeholder answer; 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): The complete routine I once had entirely skips having to define the Gram polynomials separately and is also able to produce the derivative coefficients. I'll post the routine as soon as I find that notebook... 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): 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}] // DeveloperToPackedArray]; SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___] /; 1 < n < 2 m + 1 := DeveloperToPackedArray[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[]. 1 answered Jun 28 '15 at 7:38 J. M. is away♦ 100k1010 gold badges317317 silver badges474474 bronze badges