4 faster Gram polynomial; added notes on built-in SavitzkyGolayMatrix
source | link

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 Internal`DCHypergeometricPFQ[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 Internal`DCHypergeometricPFQ[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
source | link

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}] // 
    Developer`ToPackedArray];

SavitzkyGolay[n_Integer?Positive, m_Integer?Positive, opts___]opts___?OptionQ] /; 
  1 < n < 2 m + 1 := 
 Developer`ToPackedArray[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}] // 
    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}]]

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}] // 
    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}]]
2 finally!
source | link

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

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

1
source | link