Episode #125 of the Stack Overflow podcast is here. We talk Tilde Club and mechanical keyboards. Listen now

3 of 4 added 10 characters in body

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

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[].