-1

I'm applying an Exponential Moving Average as filter for smooth params within my audio application:

a0 = 0.01
z += a0 * (input - z);

Here's the code and the firsts 50 steps:

#include <iostream>

int main ()
{
    double a0 = 0.1;
    double input = 0.8;
    double z = 0.0;

    std::cout << "init z: " << z << std::endl << std::endl;

    for(int i=0; i < 50; i++) {
        z += a0 * (input - z);
        std::cout << z << std::endl;
    }

    std::cout << std::endl << "final z: " << z << std::endl;
}

I need to check if the prev smoothed value is the same of the current one, which mean that the filter has "finished" its smothing process, and the value will always be the same.

But z will always differs an epsilon from input, so I can't check input == z it wil always be false. Here's an example, with infinite loop.

What will be the epsilon between z and input? So if its within that range, I can check and avoid further operations.

markzzz
  • 47,390
  • 120
  • 299
  • 507
  • 2
    A multiplicative epsilon between `z` and `input` might be easier to characterise, and be more computationally stable. – Bathsheba Oct 31 '17 at 14:35
  • @chux: added the example. – markzzz Oct 31 '17 at 14:35
  • @Bathsheba: z and input can be/become whatever value during the live of the application. – markzzz Oct 31 '17 at 14:36
  • That's an unusual way to describe the wanted result. Obviously, if you apply infinitely many steps, `z` will be `input`. So your value after smoothing has finished is just `input`. No need to do the iterative smoothing. What are you trying to achieve? Would it be more reasonable to apply a fixed number of iterations? – Nico Schertler Oct 31 '17 at 14:38
  • If `input` represents an actual audio stream, it's unclear what exactly it means for a filter to have "finished". On the other hand, if `input` is a fixed quantity, then we know that `z` will always converge on the value of `input`. – Oliver Charlesworth Oct 31 '17 at 14:38
  • @NicoSchertler: you are in double. At some point z will block and won't increment anymore. If it will be different than input, I won't never do a concrete comparison between z and input. – markzzz Oct 31 '17 at 14:39
  • @OliverCharlesworth: yes that's an IIR filter, but `double` in some way are "finite". So there's a stall at some point. – markzzz Oct 31 '17 at 14:41
  • I was referring to the mathematical basis, not the floating point implementation. If you are just interested in that, you could easily compare the two values before and after the modification. But as I said, you will just get the input with some floating point noise. – Nico Schertler Oct 31 '17 at 14:41
  • @NicoSchertler: this will introduce a third variable, the `prevZ`, which I'd like to avoid. – markzzz Oct 31 '17 at 14:42
  • @markzzz - Indeed. I guess I was really just trying to understand what overall problem you're trying to solve here (as it feels somewhat [XY-ish](http://xyproblem.info/) ;) That epsilon will be approximately proportional to the value of `input`, so Bathsheba's suggest is a good one. – Oliver Charlesworth Oct 31 '17 at 14:46
  • @OliverCharlesworth: which suggestion? – markzzz Oct 31 '17 at 14:59
  • 1
    philosophically speaking :), the decision to break the loop has two costs: the distance of the result from the 'true' value *and* the number of iterations; until you clarify in your mind how the total cost is formed ( and hence when the maximum cost is reached and the loop broken ) the question is and always will be ill-formed. – Massimiliano Janes Oct 31 '17 at 15:39
  • 2
    You are wrong to reject comparing the new `z` value to the previous `z` value. It is possible to write code that exactly tests whether the addition of `a0 * (input - z)` will or will not alter the value of `z` by using properties of floating-point arithmetic, but that code is much more complicated than simply copying `z` and comparing new and old values. In the absence of direct comparison or an exact test, then the alternative is simply to test whether the increment is “small” compared to `z`. Then “small” is a matter of preference or is application-specific, and there is no definite answer. – Eric Postpischil Oct 31 '17 at 15:45
  • @EricPostpischil "It is possible to write code that exactly tests whether the addition of a0 * (input - z) will or will not alter the value of z by using properties of floating-point arithmetic" can you show to me an example of this? – markzzz Oct 31 '17 at 15:52
  • just a side note: mathematically `z` can be as closed as a constant `input` as we want if we wait long enough but will never reach `input`. That's the difference between max and sup. Only in floating point can it reach `input` (sometimes). So the epsilon here is not only a floating point problem, as @MassimilianoJanes already noted. – aka.nice Oct 31 '17 at 16:46
  • 1
    @markzzz: I can but will not. It is a considerable nuisance. You can use frexp and ldexp to calculate ½ ULP of `z` and then compare `a0 * (input - z)` to it. Then, if it is exactly ½ ULP, you have to determine the parity of the low bit of the significand of `z`. And you have to deal with various special cases, such as subnormals. It would be an excessive amount of code. Comparing to the previous value of `z` is a superior solution. – Eric Postpischil Oct 31 '17 at 17:54
  • @markzzz Why are you trying to avoid the use of `prevZ`? It's clearly the most straightforward solution, and probably the most efficient one as well. – Sneftel Nov 01 '17 at 15:07
  • @Sneftel: I don't want to compare to `prev` because this means that I need, for every iteration, apply that smooth algo to the input value before comparing with `prev` one. Which is exactly what I need to avoid, saving on CPU for massive processing. – markzzz Nov 10 '17 at 11:27
  • @EricPostpischil: read above, as for Sneftel. – markzzz Nov 10 '17 at 11:27
  • 1
    @markzzz: First you said you did not want to avoid introducing a new variable, now you say you want to save on CPU. The operations required to do any of the tests suggested to you are fast, certainly faster than the operations required to compute whether `z` will change when the increment is added by examining the floating-point characteristics involved rather than directly doing the operations. You have been given good solutions by people who know what they are talking about. Use them. – Eric Postpischil Nov 10 '17 at 15:20

4 Answers4

2

Consider the ratio "new z" to "old z", less 1:

(z + a(i - z)) / z - 1

(Which obviously simplifies to ia / z - a). If the magnitude of this is less than, say 1e-6, then accept as finished. If z is zero then always continue. Adjust this multiplicative tolerance to something suitable to your requirements.

(Scientifically the tolerance will be related - dare I suggest even proportional - to the standard deviation of your data stream, but I can't offer any more hints without studying the actual data.)

Bathsheba
  • 231,907
  • 34
  • 361
  • 483
2

For an audio application, you need to consider the number of bits of sampling to know when the result will be inaudible. Each bit represents a power of 2. For example, 16 bits will be 216 or 65536, so an appropriate epsilon would be your sample scale divided by 65536. For 20 bits it's 220 or 1048576.

These limits are significantly larger than most other applications would require.

Mark Ransom
  • 299,747
  • 42
  • 398
  • 622
  • Are you saying that if my audio card is 24bit and the scale is 0.0 to 1.0, I need an epsilon of 1.0/2^24? 0,000000059604644775390625? It could make sense for "gain", but once become to filter or other param which involve not only volume, this is not related. But I don't know... – markzzz Oct 31 '17 at 15:40
  • @markzzz if I understand your problem correctly, you're trying to determine when the effect will become inaudible. It becomes inaudible when the magnitude of the effect is lower than the quantization noise of your input signal. It doesn't matter how many manipulations the signal goes through before or after this step. – Mark Ransom Oct 31 '17 at 15:47
  • 1
    This is not phrased well. The number of bits of sampling has nothing to do with what volume or change in volume the human ear can perceive. The magnitude of sound corresponding to a particular sampling value can have any scale, and there is no reason to think the low bit of the sample is the threshold of human perception. – Eric Postpischil Oct 31 '17 at 15:49
  • @EricPostpischil I didn't say anything about human perception, I said once the error is lower than the error of the input you've gone far enough. And perhaps you're correct that I phrased it badly. – Mark Ransom Oct 31 '17 at 15:51
  • @MarkRansom: You said “inaudible,” which is “unable to be heard,” which depends on the perception of the listener. We can assume there are enough bits in the samples to provide useful data, but we do not know whether there are barely enough, so that even a change in the lowest bit might be detectable by the listener, or there are plenty, so that large changes in several low bits will make no difference. Targeting the epsilon to the magnitude of the ULP is essentially saying, “Well, you have this much quantization noise already, adding that much more will not hurt.” But that is not always true. – Eric Postpischil Nov 01 '17 at 05:16
1

Instead checking the epsilon between z and input, you can check it between new value of z and the previous one.

Germán
  • 581
  • 1
  • 4
  • 19
  • This will introduce a third variable, the prevZ, which I'd like to avoid. – markzzz Oct 31 '17 at 14:46
  • It’s possible to calculate and test the increment (`a0 * (input - z)`) and only sum to `z` if it’s below `epsilon`. – Germán Oct 31 '17 at 14:52
-4

In C++ there is

std::numeric_limits<double>::epsilon()

which returns the machine epsilon, that is, the difference between 1.0 and the next value representable by the floating-point type T.

Here the modified code:

#include <iostream>

int main()
{
    int counter = 0;
    double a0 = 0.1;
    double input = 0.8;
    double z = 0.0;

    std::cout << "init z: " << z << std::endl << std::endl;

    while (true) {
        z += a0 * (input - z);
        std::cout << counter++ << " | process: " << z << std::endl;

        double eps = std::numeric_limits<double>::epsilon();
        double diff = abs(z - input);
        if (diff <= 2 * eps) {
            break;
        }
    }

    std::cout << std::endl << "final z: " << z << std::endl;
}
Simone Cifani
  • 784
  • 4
  • 14
  • This is unlikely to be a very useful approach in practice - what's important here is the [ULP](https://en.wikipedia.org/wiki/Unit_in_the_last_place), not the epsilon. – Oliver Charlesworth Oct 31 '17 at 14:48
  • @OliverCharlesworth `std::numeric_limits::epsilon()` returns the machine epsilon, that is, the difference between 1.0 and the next value representable by the floating-point type T. In my opinion it expresses the ULP. If not I am curious to understand the difference. – Simone Cifani Oct 31 '17 at 14:55
  • 1
    The problem is that the ULP is dynamic - the spacing between floating-point numbers depends on the magnitude. So you essentially need to multiply the epsilon by the current value to get the relevant ULP. – Oliver Charlesworth Oct 31 '17 at 14:58
  • 2
    Everytime I see someone blindly using `std::numeric_limits::epsilon()` as some kind of tolerance, I die a little inside. **There is no justification for doing this.** – Bathsheba Oct 31 '17 at 15:00
  • @OliverCharlesworth That's right! Thanks, I dindn't think this aspect. – Simone Cifani Oct 31 '17 at 15:00
  • What if I use std::nextafter instead, starting from the actual z, and check if it will "exceed" input? Will ULP respected in that case? – markzzz Nov 03 '17 at 22:37