The Stack Overflow podcast is back! Listen to an interview with our new CEO.

4 of 4 faster Gram polynomial; added notes on built-in SavitzkyGolayMatrix

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.

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:

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

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
 answered Jun 28 '15 at 7:38 J. M. will be back soon♦ 101k1010 gold badges319319 silver badges479479 bronze badges