From: Andre <al_fischer-zoth@t-online.de>
Subject:  DSP Trick
Date: 24 Mar 2000 00:00:00 GMT
Newsgroups: comp.dsp
THIS WORK IS PLACED IN THE PUBLIC DOMAIN
Name: Implementation and Enhancement of 'Magnitude Estimator'
Category: Algorithmic
Application: estimation of magnitude of a complex number
Advantages: fast, small, quite accurate
Introduction: please refer to the 'magnitude estimator'
The Trick:
This is an implementation of the 'magnitude estimator' on the ADSP2181.
It replaces sqrt(re^2+im^2).
There are two versions, the first one parting the 45 degree-region in two
regions, the second one dividing it in 3 regions. That means the octagon
representation of the unity-circle is replaced by a 16-gon or a 24-gon,
respectively.
The ADSP-code is quite easy to understand, so it should be easy to implement
it on other processors as well. The coefficients are in 1.15 notation, i.e.
32768 represents 1.0000.
Andre Lodwig, 3/2000
-----------------------
1st version: using 16-gon instead of the octagon:
! magnitude estimation, (re,im) in (ax0,ax1)
! result in mr1
! cycles: 18 (max)
! max error: 1.2% (-38dB)
magni:
 ar = abs ax0;
 ay0 = ar,ar = abs ax1;
 none = ar - ay0;
 if gt jump magni001;
  mx1 = ar;
    ar = ay0;
  jump magni002;
magni001:
  mx1 = ay0;     ! make mx0 the bigger one...
magni002:
 ay0 = mx1;
 my0 =  13573;      ! tan(pi/8) = 0.4142
 mr = ar * my0 (uu);
 none = mr1 - ay0;
 if ge jump magni003;
  my0 = 27491;    ! 0.41*a $lt;= b
  my1 = 18298;
  jump magni005;
magni003:       ! 0.41*a > b
  my0 = 32394;
  my1 = 6369;
magni005:
 mr = ar * my0 (ss);
 mr = mr + mx1 * my1 (rnd);
  rts;
}
2nd version: using a 24gon as a circle-approximation
! magnitude estimation, (re,im) in (ax0,ax1)
! result in mr1
! cycles: 22 (max)
! max error: 0.5% (-46dB)
magni:
 ar = abs ax0;
 ay0 = ar,ar = abs ax1;
 none = ar - ay0;
 if gt jump magni001;
  mx1 = ar;
    ar = ay0;
  jump magni002;
magni001:
  mx1 = ay0;     ! make mx0 the bigger one...
magni002:
 ay0 = mx1;
 my0 =   8779;      ! tan(pi/12) = 0.26
 mr = ar * my0 (rnd);
 none = mr1 - ay0;
 if ge jump magni004;
 my0 =  18920;      ! tan(pi/6) = 0.577
 mr = ar * my0 (uu);
 none = mr1 - ay0;
 if ge jump magni003;
  my0 = 26093;    ! 0.57*a < b
  my1 = 20007;
  jump magni005;
magni003:       ! 0.577*a > b >= 0.26*a
  my0 = 30373;
  my1 = 12583;
  jump magni005;
magni004:       ! 0.26*a >= b
  my0 = 32599;
  my1 =  4298;
magni005:
 mr = ar * my0 (uu);
 mr = mr + mx1 * my1 (rnd);
  rts;