Re: [Discuss-gnuradio] Serious bug in tag propagation with non-integer relative rate

[Andy Walls]

Hi Marcus:

On Sun, 2017-12-31 at 15:09 +0000, Müller, Marcus (CEL) wrote:
> Hi Andy, hi Eugene
> 
> Hm, coming back to an idea I had not so long ago:
> 
> tag offset should not be 64bit unsihned integers only, but also have
> a
> 64 bit fractional part.
> 
> That would not immediately solve the computational inaccurracy during
> resampling, but would at least for multi (as in >2) rate systems (eg.
> rational resampling for sample rate reduction, arbitrary resampling
> for
> clock recovery) at least give the option of having consistent
> tagging.

IMO, adding a fractional part does not solve the problem at hand -
correct tag propagation through 1 block - but it probably is required
for precise propagation through multiple rate changing blocks.

I'm keeping it in mind, but I want to avoid changing existing
interfaces right now.  I.e. get_tags_in_[range|window]() runs into a
backward compatability issue for a tag that comes < 0.5 samples before
the beginning of the requested range: should the tag offset -0.49
samples realtive to the start of the window be returned or not?


> Personally, my gut feeling is that for this kind of divisional math,
> floating point is not the number format of choice

Agreed.

>  (for multiple (as in
> > 2) reasons).
> 
> This actually doesn't apply to the arbitrary resampling case as much,
> but the d_relative_rate property of a block shouldn't be a floating
> point,

Internal to the gnuradio-runtime there are a number of places where,
for estimating, using a floating point d_relative_rate is just fine. 
It is expedient and avoids the problem of interger overflow for cases
where one has large integer interpolation and decimation factors.

> but a ratio of two integers, one of them being signed. Maybe I
> should be adding a set_relative_rate(int64_t numerator, uint_64_t
> denominator) to block, and make set_relative_rate(double) a wrapper
> for
> that (and of course change the block_executor to do fractional rather
> than floating point math);

So that's what my change does. :)  I also changed all the blocks, that
could easily easily call the integer version of
set_relataive_rate(interpolation, decimation), to do so.

To avoid integer overflow, I picked an algoritm in the new
set_relative_rate(double) wrapper, to convert the relative rate to a
ratio of two integers that got within 1 ppm of the passed in relative
rate.  Because sane users only use relative rates of 1/25 or 1001/1000
or 255/256 or 2/3 or something like that, and whacky relative rates
using large (>16 bit) integers in the ratio are just a product of
rounding, right?  Well I was wrong. :(

Unfortunately, the fractional resampler is a corner case that actually
uses the rounded (in binary) version of the resampling rate, that the
user specifies, in its actual operation. Computing exact integer ratios
to represent the double can lead to 53 bit numbers (IIRC), which can
lead to integer overflows when multiplying uint64_t's together. :(

So what I'm left with are two broad options to deal with handling a
double relative rate (like the fractional resampler uses):

1. Implement an method like Eugene suggests: do a double multiply, get
a rounded result, do a double division to get back to an initial
offset, subtract off that initial offset, do a double multiply again,
and then a final add.

or

2. Allow and compute large, exact integer ratios of
interpolation/decimation rate and use multiple precision (128 bit)
arithmetic when overflow out of 64 bits is possible.
Maybe using boost::multiprescision with the MIPR library on the
backend.
http://www.boost.org/doc/libs/1_66_0/libs/multiprecision/doc/html/index
.html
http://mpir.org/
(MPIR is supposedly better for Windows compatability than GMP.)


I'm leaning towards option 2.

> I find that I never calculate a rate other
> than by calculating the floating point approximation of a fraction of
> integers, and that there's always been subtle problems when running
> odd
> ratios for prolonged times. 

Technically a finite precision double can always be represented by a
ratio of integers, but those integers can be very (2^53-ish) large. 

Thanks for the feedback.

Regards,
Andy

> Best regards,
> Marcus
> 
> On Sat, 2017-12-30 at 18:24 -0500, Andy Walls wrote:
> > Hi Eugene:
> > 
> > On Wed, 2017-12-27 at 16:18 -0500, Andy Walls wrote:
> > > Hi Eugene
> > > 
> > > > From:   Eugene Grayver
> > > > Date:   Thu, 9 Nov 2017 19:52:35 +0000
> > > > 
> > > > There is a major problem with the way tags are propagated in
> > > > blocks
> > > > with non-integer relative rate. If the ratio is not a power of
> > > > two,
> > > > the numerical accuracy of the floating point will cause the
> > > > output
> > > > tags to diverge from the input tags.  Consider the fractional
> > > > resampler. It accumulates the timing offset in the range of 0
> > > > to
> > > > 1.
> > > > However tag propagation multiplies the ratio by the sample
> > > > number.
> > > > As
> > > > the sample number grows the LSB accuracy of the ratio gets
> > > > scaled
> > > > by
> > > > the ever larger value. For a ratio of 1.001 we saw divergence
> > > > of
> > > > 1000s of samples over a few minutes at 10msps.
> > > 
> > > Could you please test the following branch to see if it fixes the
> > > problem?  Maybe test something simple first, like an FIR filter
> > > decimating by 5 or 3?
> > > 
> > > https://github.com/awalls-cx18/gnuradio.git branch: tag_fix3
> > 
> > Don't bother testing.  See below.
> > 
> > > Or if you have a GRC or python script I can use myself for
> > > testing,
> > > that would be great.
> > > 
> > > >  I rewrote tag propagation for the resampler but did not rework
> > > > the
> > > > generic logic. I think the key point is to use the delta
> > > > between
> > > > read
> > > > and written items to take out the large integer difference and
> > > > then
> > > > apply the scaling to a local delta within the current window.
> > > > 
> > > 
> > > The fix that I have made stores the relative_rate as an integer
> > > numerator and an integer denominator, and it uses integer
> > > arithmetic
> > > to
> > > propagate tags. (Except if enable_update_rate() is True, in which
> > > case,
> > > precision tag placement was abandonded by the block author
> > > anyway.)
> > 
> > So this fix makes the fraction resampler tag propagation actually
> > perform worse.  The reason appears to be that, at least for
> > resamp_ratios very close to 1.0, the fractional resampler isn't
> > really
> > running at the requested rate, but some rounded (in binary) rate.
> > 
> > So I have a test flowgraph with a single fractional respampler in
> > it,
> > with a resample ratio of 1.001 specified.  Here are the parameters
> > reported by my patched GNURadio:
> > 
> >     [    address@hidden     grcs]$ ./tag_prop_test.py 
> >     Fractional Resampler resamp ratio:  1.00100004673
> >     Block relative rate:  0.999000952364
> >     Block relative rate i:  1000
> >     Block relative rate d:  1001
> > 
> > So we can see that the block is really running with a resampling
> > ratio
> > of 1.00100004673.  The relative rate of 0.999000952364 does appear
> > to
> > be the correct reciprocal of 1.00100004673.  My patch's back-
> > figuring
> > of a relative rate of 1000/1001 (= 0.999000999 = 1/1.001) is in
> > fact
> > what the user wanted.  But when propagating tags using those
> > integers,
> > tags start sliding very soon, since 0.999000999 is not the
> > 0.999000952364 that the block is actually operating at.
> > 
> > I'll have to think about how to deal with this sort of situation.
> > 
> > Regards,
> > Andy
> > 
> > _______________________________________________
> > Discuss-gnuradio mailing list
> > address@hidden
> > https://lists.gnu.org/mailman/listinfo/discuss-gnuradio