1

I'm currently researching the Kolmogorov-Zurbenko filter and trying to implement it myself as a way to smooth one-dimensional signal strength values.

The basic filter per se is pretty easy to implement by using two loops:

for each iteration do:
    for each point do:
        movingAverage(point)

Since this filter is going to be applied quite often and performance is an issue, I'd like to precalculate the coefficients $a^{m, k}_s$ - see https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Zurbenko_filter#Definition - once so that the iteration loop may be replaced by one simple multiplication (see the last few lines in the definition section).

To do this, $c^{k, m}_s$ has to be calculated: $$a^{m, k}_s = \frac{c^{k, m}_s}{m^k}$$ The problem is that I have trouble understanding how to do that.

The definition section of the linked wiki article states that $c^{k, m}_s$ may be obtained by the equation following, but when I try to remove the sum and the factor $z$ (since I want to calculate $c$ for one specific values of $s$ respectively $r$) from this equation, I end up with $c^{k, m}_s = 1$, regardless of the parameters.

Obviously I'm missing something here - I'd appreciate any hints. Thanks!

justonemorething
  • 113

2 Answers2

1

From that link, you see that the terms $c^{k,m}_s$ are defined as the coefficients of a polynomial: $$ \sum^{k(m-1)}_{r = 0}z^rc^{k, m}_{r - k(m - 1)/2}=(1 + z + ... + z^{m-1})^k $$ which means that to compute them you need to expand the power $(1 + z + \dotsb + z^{m-1})^k$ for appropriate values of $k$ and $m$.

The case $m = 2$ is fairly simple, because then by the binomial theorem you know that $$ (1+z)^k = \sum_{r = 0}^k \binom{k}{r} z^r $$ hence $$ c^{k,2}_{r - k/2} = \binom{k}{r} := \frac{k!}{r! (k-r)!} $$

For a simple example, consider $k = 3$, $m = 2$. Then $$ (1 + z + ... + z^{m-1})^k = (1 + z)^3 = 1 + 3 z + 3 z^2 + z^3 $$ means that $$ c^{3,2}_{\pm 3/2} = 1 \qquad c^{3,2}_{\pm 1/2} = 3 $$

For the general case, it may be useful to observe that $c^{k,m}_s = c^{k,m}_{-s}$.

Also, according to this answer you can use the discrete Fourier transform to speed up the computation of the powers of a polynomial.

Community
  • 1
A.P.
  • 9,728
  • Great, thanks! Regarding the answer you linked: The way I understand it, the equation $g_k = ...$ should calculate the coefficients I need (in the case of $k = 3$ and $m = 2$: 1, 3, 3, 1. Is this correct? – justonemorething Jul 05 '15 at 18:48
  • Yes, that is essentially correct. To be clear, while in the author of that answer is showing how to compute the coefficients of $g = f^2$ for some polynomial $f$, here you need to compute the coefficients of $g = f^k$ where $f = 1 + z + \dotsb + z^{m-1}$. I trust that you are able to generalise that method accordingly. – A.P. Jul 05 '15 at 19:19
  • I'd say that $$g = f^a = \sum\limits_{k=0}^{i=an} g_k x^k$$ and furthermore $$g_k = \sum\limits_{i+j+...a=k,0 \le i,j,...a \le n} f_i f_j ... f_a$$. Is that right? Also, since I haven't seen this notation before: I don't know how to interpret a sum sign with no upper bound and multiple "indices". If it's not too much to ask: Could you possibly give a short example of how to calculate the new coefficients with, say, $f = 1 + z$ to the second power? – justonemorething Jul 05 '15 at 23:39
  • I added an explicit formula for $(1+z)^k$. By the way, $(1+z)^2 = 1 + 2z + z^2$. Anyway, when you see a summation (or product) symbol without superscripts, you should interpret it as "sum over the indices satisfying these conditions". In the case you mention, the sum is over all the indices $i,j,\dotsc,a$ between $1$ and $n$ whose sum is $k$. Anyway, for small values of $k$ you may probably still find the multinomial theorem useful. – A.P. Jul 06 '15 at 21:58
  • By the way, if you just wish to pre-compute a handful of coefficients you may find it easier to just use Wolfram Alpha. – A.P. Jul 06 '15 at 22:01
  • Ah, I see. Thanks again! Now I get it. Regarding the pre-computed coefficients: I might do that, but it's nice to know how to do this myself. – justonemorething Jul 09 '15 at 10:29
0

Actually, to factor filters tend to speed up calculations, so you don't have to multiply together to determine the coefficients. I think it would be possible to make faster by applying the uniform filter several times (if implemented fast), especially as the moving average filter has very little overhead. You can store sum of m values, then for the next value in the signal remove the leftmost and add the rightmost - one addition and one subtraction, instead of having to multiply and ackumulate $m$ numbers with a filter.

mathreadler
  • 25,824
  • That's a good idea, I'll try that. – justonemorething Jul 09 '15 at 10:32
  • Even faster is to calculate an integral signal as a cumulative sum. Then the mean value filters to iterate become [1,zeroes,-1], just 2 nonzero taps. And their convolutions N times become N+2 taps partial differential filters of order N-1 with the taps being binomial expansion of (z-1)^N. – mathreadler Jul 09 '15 at 11:18