On Nov 30, 2013, at 8:22, Tom Rondeau <tom@trondeau.com> wrote:
> On Sat, Nov 30, 2013 at 10:55 AM, Kevin Reid <kpreid@switchb.org> wrote:
>> /* don't divide by zero! */
>> if((y_abs < 1.5E-5) && (x_abs < 1.5E-5))
>> return 0.0;
>>
>> Since the preceding multiply_conjugate should square the magnitude, this choice of epsilon accounts for the threshold I found empirically.
>>
>> Given this, I'm going to file an issue requesting that this epsilon be reduced. (The only alternative that comes to mind is to place an AGC block before the demodulator, which seems unlikely to be better from a numerical perspective).
>
> Agreed. This threshold should be reduced. Can you experiment with that
> threshold to see if there's any lower limit that might be problematic
> (my default threshold, for some reason, is 1e-20; would want to make
> sure that doesn't cause something else in the math to blow up).
I have modified a copy of fast_atan2f to have the test "(y_abs == 0.0f) && (x_abs == 0.0f)" and plotted the results of fast_atan2f(k*sin(angle), k*cos(angle)) against both angle and atan2f(k*sin(angle), k*cos(angle)) at various exponents for k.
At small magnitudes fast_atan2f is a _better_ approximation to regular atan2f than it is at magnitudes close to 1, and the error compared to the original angle (which should partly be attributed to representing sin and cos values) has the same patterning at all scales. At no point does fast_atan2f emit complete nonsense; it is consistent all the way down to the point at which both y and x are rounded to 0.0.
My conclusion is that fast_atan2f is as good as it ever is all the way down to the minimum representable values (~ 1e-45), and that one might as well revert the code to an equality test with 0.0 (or perhaps "y_abs <= 0.0" if a linter needs to be satisfied).
That said, I'm no expert on numerics and I might well have missed some particular undesirable behavior that should be avoided. Therefore, I will attach the tools I wrote to test the behavior. "main.cc" generates a table of error values (the parameters are hardcoded), and "run.sh" will build it and plot the values using the gnuplot script.
(Note that there is a special case in main.cc to emit a constant nonzero value when x and y are exactly zero after being scaled down; that value is intended as an out-of-band marker and is not an actual error.)
No comments:
Post a Comment