I'm just posting this to record for posterity [something I posted in the chatroom](http://chat.stackexchange.com/transcript/message/4829800#4829800) 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](http://pubs.acs.org/doi/abs/10.1021/ac00205a007) (though I have traced the spirit of the algorithm going as far back as [Hildebrand's book](http://books.google.com/books?id=hqucruPBheQC&pg=PA357)):

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