How to get "reliable" rounding?

Here are some sample times in seconds from some equipment nominally sampled at just below 500Hz (~2ms per sample), the times are stored with 6 decimal places:

x = [650.979531;650.981479;650.983426
    650.985373;650.987321;650.989267
    650.991215;650.993162;650.995109
    650.997057;650.999003;651.000951
    651.002897;651.004846;651.006792
    651.008738;651.010688;651.012633
    651.014580;651.016528;651.018474];

MATLAB loads it as double precision. Lets look at the diff to see the increments in milliseconds (*1e3):

>> diff(x*1e3)
ans =
   1.948000000091270
   1.946999999927357
   1.947000000043772
   1.947999999974854
   1.945999999996275
   1.947999999974854
   ...

Looks good, every step is ~1.9ms so if I round to the nearest millisecond I naively expect it to all round to 2ms steps. But round(x, 3) [3 decimal places] doesn't get what I naively expect, note I inserted a < to the millisecond value just to highlight the millisecond value:

>> round(x,3)
ans =
   1.0e+02 *
   6.50980<0000000000
   6.50981<0000000000
   6.50982<9999999999
   6.50985<0000000000
   6.50986<9999999999
   ...

0, 1, 2, 5, 6 the rounding is not what I expect (and my analysis code expects times to have 2ms inter-sample interval). I tried to first convert to single precision (given the input has 6 decimals) but that didn't help. I also naively tried "succesive" rounding (6to5, 5to4, 4to3) hoping to accumulate each decimal:

>> round(round(round(x,5),4),3)
ans =
   1.0e+02 *
   6.50980<0000000000
   6.50981<9999999999
   6.50982<9999999999
   6.50985<0000000000
   6.50986<9999999999
   ...

I don't know why I get values like 6.509829999999999 when I expect 6.509830000000000. I also tried to convert to int32 with no luck:

   >> int32(x*1e3)
ans =
  21×1 int32 column vector
   650980
   650981
   650983
   650985
   650987
   650989
   650991
   650993
   650995
   650997
   650999
   651001
   651003
   651005
   651007
   651009
   651011
   651013
   651015
   651017
   651018

Mostly steps of 2ms with occasional steps of 1ms (first and last step). If I convert to single then int32 [int32(single(x)*1e3)] I still get mixed results. I'm sure that is my ignorance of the core of computer science and floating point values, but I just want to find a way if possible to reliably round to the nearest 2ms (given the clear diff values). Any help much appreciated, thanks!

Answer 1

[Rik]

Your problem is not with how floats work, but with how decimal values work.

If you look at the calender, you might notice something weird: February sometimes has an extra day. That's because the sidereal day does not match the solar day. It is close though: 23h56m04s. If you round that, you would get to 1 day. So if you now do what you did to your data, that would mean the leap days are no longer needed. Good news for people who can't remember the rules, but bad news for people who don't want the seasons to have moved 100 days since the early 1600s.

You can round your data to 2ms intervals OR to the closest integer milisecond, but not both. The choice is yours.

Alternatively, you can resample your data to match the timing interval you want.

[Ian]

Thank you, resampling may make the most sense in this case, and I appreciate you sharing your knowledge.

Answer 2

[Walter Roberson]

temp = 6.50983000000000

temp = 6.5098

fprintf('%.999g\n', temp)

6.50983000000000000540012479177676141262054443359375

You can see that 6.50983000000000 exactly does not exist in binary floating point. Values must always be higher or lower than that value.

>> diff(x*1e3)

ans =

6.50983000000000

We can tell from that output that you are using format long g . But here's the thing: format long g uses %.15g . The majority of double precision numbers need 16 decimal digits to exactly replicate, with their being portion that need 17 decimal digits to exactly replicate. So where 6.50983000000000 shows up in output, the value might be less than or greater than the value you were hoping for

format long g
t2 = temp*(1-eps);
t2

t2 = 
                   6.50983

fprintf('%.999g\n', t2);

6.50982999999999822904328539152629673480987548828125

t3 = temp*(1+eps);
t3

t3 = 
                   6.50983

fprintf('%.999g\n', t3)

6.50983000000000178175696419202722609043121337890625

t2 - temp

ans = 
     -1.77635683940025e-15

t3 - temp

ans = 
      1.77635683940025e-15

t3 - t2

ans = 
       3.5527136788005e-15

The value, t2, is slightly less than the target value, but displays the same as the target value does. The value, t3, is slightly more than the target value, but displays the same way that the target value does.

... and that means that where 6.50983000000000 showed up in your output, it might have been a value slightly smaller than you expected, one that does not round the way you expect.

[Ian]

That is a very clear explanation, thank you so much! It is somewhat annoying I can only accept one answer as both your answer and Riks contributed helpful knowledge to myproblem and my lack of understanding...

[Rik]

You can still give this answer an upvote if you want to give explicit recognition.

[Ian]

Done, thanks Rik!