An Efficient Algorithm for Magnitude Operation

• Conversion of signals in time domain into spatial domain (for different coordinate systems).
• Conversion of signals from one coordinate to another in spatial domain.

• Conversion of spatial domain to rasters. These operations always involve high-speed magnitude computations.

Magnitude operation is generally implemented as lookup tables in real-time applications, with all possible combinations of input arguments as address locations. However, these lookup tables consume huge amounts of memory. While high-speed methods are available, the results they produce are only approximations for magnitude. In "A Baker's Dozen Magnitude Approximations and Their Detection Statistics," (IEEE Transactions on Aerospace and Electronic Systems, January 1976), A.E. Filip notes 13 approximations for magnitude in which the least accurate has a peak error of 27.3 percent and the most accurate has a peak error of 0.97 percent. In this article, I present an alternative approach to magnitude operation that is both fast and efficient in terms of memory. This approach, which is based on the Successive Approximation principle, predicts each bit of the result from MSB to LSB, and corrects it based on feedback. In addition, I've implemented this algorithm in C++ source code; see Listing One.

Algorithm Description

The algorithm itself is fairly straightforward. Note that in this article, I represent magnitude as M(a,b). If Z is a complex number given by Z=a+i×b where i is imaginary unity, and Z' is its conjugate given by Z=a-i×b, then Z×Z'=M×M. Furthermore, the magnitude of two numbers (a,b) is defined as SquareRoot(a×a+b×b).

For two integers (a,b), M(a,b) have the boundaries in Figure 1, where a is the larger number, see Figure 1(a), and b is the smaller number; see Figure 1(b). You define r such that r=M-a; as in Figure 1(c). This implies Figure 1(d). Therefore, if r is computed successively, such that b2=r2 +2×a×r, then (r+a) will give M(a,b). To compute r successively, in each iteration, 1 bit is added to it. Assume this increment is d. This yields Figure 1(e). Figure 2 provides a proof of the theorem in Figure 1(b).

The steps to compute M(a,b) are:

1. Compute b² using shift and add operations. This takes N iterations, where N is the number of bits in b.

2. Set the initial value of increment d to 2ⁿ, where n=(N-1)^th bit.

3. Set r=0.

4. Compute (r+d)²+2×a×(r+d), by computing (d²+2×r×d+2×a×d) and adding to previous value of (r²+2×a×r). Since d is always a power of 2, this can be computed using shift and add operations.

5. If (r+d)²+2×a×(r+d)<= b², update r to r+d; or else do not update.

6. Reduce the increment d to next bit -- 2^(n-1).

7. Repeat steps 4 to 6 until d becomes 1, (N-1 iterations).

8. Return M=r+a.

You can see that the number of iterations depends on the smallest number b, and the result is accurate to 1 LSB. Again, Listing One is C++ code that implements this. In this example, validation of input arguments are omitted.

DDJ

Listing One

unsigned long imag(unsigned int a, unsigned int b)
{
 unsigned int  l, r, d;
 unsigned long d2, b2, R, D;
 R = 0;
 l = 0;
 b2 = 0;
 r  = b;
 do
 {
  d = 1L << l;
  if(b & d)
  {
   b2 += (b << l);
  }
  l++;    // find the number of bits in b
 }while(r >>= 1);
 l--;
 while(l --)
 {
  d   = 1L << l;   // setting the increment
  d2  = d << l;
  D   = d2 + ((a + r) << (l + 1));// D = d**2 + 2*r*d + 2*a*d;
  if((R + D) <= b2)
  {
   r += d;   // updating value of r.
   R += D;
  }
 }
 return (r + a);
}

Back to Article