2 of 4 finally!

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___] /; 
  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[].