Prev: Timing Recovery for Blind Demodulation
Next: Using FFTW 3 DLL w/ CLR
From: Robert Scott on 4 Oct 2006 12:04
I have an approximation for the Pythagorean distance formula (magnitude of
vector [x,y]) using integer arithmetic that I would like to improve. Currently
I am using this:

if(x<0) x=-x;
if(y<0) y=-y;
if(x < y)
{
int t = x;
x = y;
y = t; // ensures that x >= y
}
z = (y < ((13107 * x)>>15)) ? // * (.4)
(x + ((y * 6310)>>15)) : // * (.192582403)
(((x * 27926)>>15) // * (.852245894)
+ ((y * 18414)>>15)); // * (.561967668)
//..(linear approximation to within 2% of the Pythagorean distance formula)..

This is for an ARM processor in an application where I cannot afford the time
for any floating point operations. The integer values of x and y come from an
array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
above formula. The final result (z) will be shifted right one bit before being
stored back into an array of 16-bit integers, since the distance formula can
potentially extend the range of x and y by one bit. In case you are wondering,
this is part of calculating the power spectrum at the end of an FFT.

The improvement I am looking for is in accuracy. I would like to try for a
4-fold improvement in accuracy (.5%) without substantially increasing the
running time of what I have now. Does anyone know of a better approximation
that is almost as fast?



Robert Scott
Ypsilanti, Michigan
From: Peter Nachtwey on 4 Oct 2006 14:44

Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve. Currently
> I am using this:
>
> if(x<0) x=-x;
> if(y<0) y=-y;
> if(x < y)
> {
> int t = x;
> x = y;
> y = t; // ensures that x >= y
> }
> z = (y < ((13107 * x)>>15)) ? // * (.4)
> (x + ((y * 6310)>>15)) : // * (.192582403)
> (((x * 27926)>>15) // * (.852245894)
> + ((y * 18414)>>15)); // * (.561967668)
> //..(linear approximation to within 2% of the Pythagorean distance formula)..
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations. The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula. The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit. In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy. I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now. Does anyone know of a better approximation
> that is almost as fast?
>
>
>
> Robert Scott
> Ypsilanti, Michigan

I still like the tried and true isqrt((x^2+y^2)<<16). You can see we
scale up the sum of x^+y^ by shifing left 16 times so we have a fixed
point number with the interger in the upper 16 bits. When we take the
square root of at 16.16 we get an 8.8 result. This meets your .5%
criteria. There are plenty of fast and simple square root routines.
We use a table to find the must significant 8 bits and then do one
Newton iteration to get the lower 8 bits but then I have a divide.
If you don't there are other fast and simple routines on the internet.
I think our isqrt routine takes about 8 microseconds on a 40 Mhz 196.
The squaring and summing would take two extra microseconds.

Peter Nachtwey

From: dbell on 4 Oct 2006 20:12
Robert,

I do not have the book handy, but my recollection is that methods very
similar to the method you are using, with varying levels of accuracy,
are detailed in 'Digital Signal Processing in Communications Systems'
by Marvin E. Frerking. You might look there and see if you can find
something useful.

Maybe someone here who has access to the book could see if my
recollection is correct and help you out.

Dirk Bell
DSP Consultant


Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve. Currently
> I am using this:
>
> if(x<0) x=-x;
> if(y<0) y=-y;
> if(x < y)
> {
> int t = x;
> x = y;
> y = t; // ensures that x >= y
> }
> z = (y < ((13107 * x)>>15)) ? // * (.4)
> (x + ((y * 6310)>>15)) : // * (.192582403)
> (((x * 27926)>>15) // * (.852245894)
> + ((y * 18414)>>15)); // * (.561967668)
> //..(linear approximation to within 2% of the Pythagorean distance formula)..
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations. The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula. The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit. In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy. I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now. Does anyone know of a better approximation
> that is almost as fast?
>
>
>
> Robert Scott
> Ypsilanti, Michigan

From: Phil Frisbie, Jr. on 4 Oct 2006 21:27
Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve. Currently
> I am using this:
<snip>
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations. The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula. The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit. In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy. I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now. Does anyone know of a better approximation
> that is almost as fast?

I like this algorithm. It is snipped from my OpenLPC fixed point codec.
It uses only simple operations. It is likely more accurate than you
need, but I think you can simply truncate the algorithm sooner for less
precision. At the end PRECISION is the fractional bits you are using in
your fixed point code.

static fixed32 fixsqrt32(fixed32 x)
{

unsigned long r = 0, s, v = (unsigned long)x;

#define STEP(k) s = r + (1 << k * 2); r >>= 1; \
if (s <= v) { v -= s; r |= (1 << k * 2); }

STEP(15);
STEP(14);
STEP(13);
STEP(12);
STEP(11);
STEP(10);
STEP(9);
STEP(8);
STEP(7);
STEP(6);
STEP(5);
STEP(4);
STEP(3);
STEP(2);
STEP(1);
STEP(0);

return (fixed32)(r << (PRECISION / 2));
}


--
Phil Frisbie, Jr.
Hawk Software
http://www.hawksoft.com
From: Jerry Wolf on 4 Oct 2006 22:49
The comp.dsp "DSP Tricks" at
http://www.dspguru.com/comp.dsp/tricks/tricks.htm contains two entries
on magnitude estimation, which is the same computation. In particular,
this one http://www.dspguru.com/comp.dsp/tricks/alg/mag_est.htm should
tell you all that you want to know (altho IIRC, the C program contains
a typo -- the coefficients Alpha and Beta are printed backwards in the
printed table (which one multiples the max and which the min?)) .

cheers,
jerry

Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve. Currently
> I am using this:
>
> if(x<0) x=-x;
> if(y<0) y=-y;
> if(x < y)
> {
> int t = x;
> x = y;
> y = t; // ensures that x >= y
> }
> z = (y < ((13107 * x)>>15)) ? // * (.4)
> (x + ((y * 6310)>>15)) : // * (.192582403)
> (((x * 27926)>>15) // * (.852245894)
> + ((y * 18414)>>15)); // * (.561967668)
> //..(linear approximation to within 2% of the Pythagorean distance formula)..
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations. The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula. The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit. In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy. I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now. Does anyone know of a better approximation
> that is almost as fast?
>
>
>
> Robert Scott
> Ypsilanti, Michigan

 |  Next  |  Last
Pages: 1 2 3
Prev: Timing Recovery for Blind Demodulation
Next: Using FFTW 3 DLL w/ CLR