1
$\begingroup$

I'm an astronomer trying to smooth variable star data, and one way I'm doing this is using a 7-point, second-order Savitzky-Golay filter. I iteratively apply the filter 51 times (i.e. I apply the filter, then apply the filter again to the result of the first application, etc.) I've found that the equation the filter uses to modify the data for 7-point quadratic SG filter is:

$y_t = \frac{-2x_{t-3} + 3x_{t-2} + 6x_{t-1} + 7x_t + 6x_{t+1} + 3x_{t+2} - 2x_{t+3}}{21}\tag{1}$

Error propagating Equation 1 gives:

$\sigma_{y_t} = \frac{1}{21}\sqrt{4\sigma_{x_{t-3}}^2 + 9\sigma_{x_{t-2}}^2 + 36\sigma_{x_{t-1}}^2 + 49\sigma_{x_t}^2 + 36\sigma_{x_{t+1}}^2 + 9\sigma_{x_{t+2}}^2 + 4\sigma_{x_{t+3}}^2}\tag{2}$

However, I seem to be making some sort of mistake in doing this. After iteratively applying this 51 times to my observational errors, it reduces the errors by 13 orders of magnitude! This is clearly absurd and isn't an error I can seriously report in my paper, but I can't figure out my error. I had the same problem with my convolution method, the error of which I worked out to be:

$\sigma_{y_t} = \frac{1}{w}\sqrt{\sum_{i=0}^w{\sigma_{x_i}^2}}\tag{3}$

where $w$ is the window size, for which I also used 7, giving an equation similar to Equation 2 without the coefficients in front of the errors. This also reduced the error by about 13 orders of magnitude, and through guesswork, I discovered that:

$\sigma_{y_t} = \sqrt{\frac{1}{w}\sum_{i=0}^w{\sigma_{x_i}^2}}\tag{4}$

gave errors of the proper order of magnitude. A similar method doesn't work for Equation 2, however. Since I don't understand why Equation 4 is right and Equation 3 is wrong (I taught error propagation to first-year physics undergrads as part of my GSA), I can neither discover the proper method to error propagate nor explain it in my dissertation.

I think the unknown factor is that I'm iterating the procedure (the 51 times I apply the filter), and I think that is altering the error propagation method somehow. A Google search on 'error propagation for iterated functions' was unhelpful, and I wasn't able to find a similar question on any StackExchange. Some guidance would be appreciated.

Improve this question
$\endgroup$
12
  • 2
    $\begingroup$ You can't hope for a crude approximation to remain any good after so many iterations. You can compute iterated filters directly through convolution. Indeed, that would be a much more efficient way to carry out your process, because you need to compute the coefficients only once and apply that filter just once, rather than filtering the series 51 times. You can further simplify this by noting that the result is going to be extremely well approximated by a Gaussian filter. $\endgroup$
    – whuber ♦
    Commented Jun 2, 2021 at 23:04
  • 1
    $\begingroup$ Would you mind sharing how the data looks like? SG may not be the right choice if it is producing crazy results. I have never heard of 51 iterations of SG on the same data set. May be you can increase the polynomial order of even increase the number of odd points. $\endgroup$
    – ACR
    Commented Jun 3, 2021 at 2:01
  • $\begingroup$ I wouldn't mind sharing it, as it's pretty publicly available anyway. I'd like to clarify, though, that SG itself (Equation 1) isn't producing crazy results. The 51 iterations (and actually, I think it's 50 upon further reflection) produces a nice, smooth curve. What's producing crazy results is the error propagation through the SG filter (Equation 2). $\endgroup$
    – mknote
    Commented Jun 3, 2021 at 12:19
  • $\begingroup$ @whuber I'm a bit confused. Which part of my method is a crude approximation? $\endgroup$
    – mknote
    Commented Jun 3, 2021 at 13:38
  • 1
    $\begingroup$ @whuber Numerical Recipes was what I needed. I now understand how to do as you suggested: find the DFT of a vector of Eqn. 1's coefficients and multiply it by the DFT of what I'm filtering. To do n iterations, I raise the coefficient vector to the power of n. Then I use the inverse DFT to get the filtered values. Turning back to error propagation, the first answer of this question implies that I do the same thing, but with the errors of the values instead of the values themselves. Is that correct? $\endgroup$
    – mknote
    Commented Jun 13, 2021 at 23:41

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.