DSP Trick: Magnitude Estimator

From: Grant Griffin
Subject: DSP Trick: Magnitude Estimator
Date: 02 Oct 1999 00:00:00 GMT
THIS WORK IS PLACED IN THE PUBLIC DOMAIN

Name: Magnitude Estimator

Category: Algorithmic

Application: Signal strength measurement, AM demodulation, etc.

Advantages: This estimation algorithm is very quick compared to calculating magnitude to full precision using a square-root operation.

Introduction:

Given a complex number represented as "I + jQ", the magnitude is: sqrt(I**2 + Q**2). Unfortunately, the square-root operation in this formula is inherently slow and complicated to calculate, either in software or in hardware.

For applications which do not require a full-precision magnitude value, the use of a magnitude _estimation_ can save calculations. This well-known and widely used algorithm is a true gem of DSP because it provides considerable savings in calculation, at the cost of only a minimal loss of accuracy.

The Trick:

Given a complex number whose real part is "I" and whose imaginary part is "Q", the algorithm estimates the magnitude as:

Mag ~=Alpha * max(|I|, |Q|) + Beta * min(|I|, |Q|)

"Alpha" and "Beta" are two constants whose values can be chosen to trade among RMS error, peak error, and implementation complexity.

How It Works:

The absolute value operations folds the complex number into the range of 0-90 degrees, and the min, max operations further fold the complex number into the range of 0-45 degrees. Within this limited range, a linear combination of I and Q are a good approximation of magnitude.

Values for Alpha and Beta:

A program to demonstrate and test this algorithm is given at bottom. It includes a table of the most useful values for Alpha and Beta. The program prints out the following:

=====================================================================
             Alpha * Max + Beta * Min Magnitude Estimator

Name                  Alpha           Beta       Avg Err   RMS   Peak
                                                 (linear)  (dB)  (dB)
---------------------------------------------------------------------
Min RMS Err      0.947543636291 0.392485425092   0.000547 -32.6 -25.6
Min Peak Err     0.960433870103 0.397824734759  -0.013049 -31.4 -28.1
Min RMS w/ Avg=0 0.948059448969 0.392699081699   0.000003 -32.6 -25.7
1, Min RMS Err   1.000000000000 0.323260990000  -0.020865 -28.7 -23.8
1, Min Peak Err  1.000000000000 0.335982538000  -0.025609 -28.3 -25.1
1, 1/2           1.000000000000 0.500000000000  -0.086775 -20.7 -18.6
1, 1/4           1.000000000000 0.250000000000   0.006456 -27.6 -18.7
Frerking         1.000000000000 0.400000000000  -0.049482 -25.1 -22.3
1, 11/32         1.000000000000 0.343750000000  -0.028505 -28.0 -24.8
1, 3/8           1.000000000000 0.375000000000  -0.040159 -26.4 -23.4
15/16, 15/32     0.937500000000 0.468750000000  -0.018851 -29.2 -24.1
15/16, 1/2       0.937500000000 0.500000000000  -0.030505 -26.9 -24.1
31/32, 11/32     0.968750000000 0.343750000000  -0.000371 -31.6 -22.9
31/32, 3/8       0.968750000000 0.375000000000  -0.012024 -31.4 -26.1
61/64, 3/8       0.953125000000 0.375000000000   0.002043 -32.5 -24.3
61/64, 13/32     0.953125000000 0.406250000000  -0.009611 -31.8 -26.6
=====================================================================

Variations:

If you need a more accurate estimation than this algorithm provides, you can use some variation of it. For example, varying values of Alpha and Beta can be taken from a small lookup table, driven by the relative size of min and max values. Another possibility is to use this estimate as the "seed" of an iterative magnitude estimator.

Alternatives:

As an alternative, consider using the CORDIC algorithm, especially in hardware applications.

References:

This algorithm is described in Understanding Digital Signal Processing by Richard G. Lyons [Lyo97], in Digital Signal Processing in Communication Systems by Marvin E. Frerking [Fre94], and elsewhere.

Thanks to Clay S. Turner for providing some of the coefficients.

The Code:

#include <math.h>
#include <stdio.h>
/*********************************************************************
*
*
* Name: mag_est.c
*
* Synopsis:
*
*   Demonstrates and tests the "Alpha * Min + Beta * Max" magnitude
*   estimation algorithm.
*
* Description:
*
*   This program demonstrates the "Alpha, Beta" algorithm for 
*   estimating the magnitude of a complex number.  Compared to
*   calculating the magnitude directly using sqrt(I^2 + Q^2), this
*   estimation is very quick.
*
*   Various values of Alpha and Beta can be used to trade among RMS
*   error, peak error, and coefficient complexity.  This program
*   includes a table of the most useful values, and it prints out the
*   resulting RMS and peak errors.
*
* Copyright 1999  Grant R. Griffin
*
*                    The Wide Open License (WOL)
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without
* fee, provided that the above copyright notice and this license
* appear in all source copies. THIS SOFTWARE IS PROVIDED "AS IS"
* WITHOUT EXPRESS OR IMPLIED WARRANTY OF ANY KIND. See 
* http://www.dspguru.com/wol.htm for more information.
*
*********************************************************************/
/********************************************************************/
double alpha_beta_mag(double alpha, double beta, double inphase,
                      double quadrature)
{
   /* magnitude ~= alpha * max(|I|, |Q|) + beta * min(|I|, |Q|) */
   double abs_inphase = fabs(inphase);
   double abs_quadrature = fabs(quadrature);
   if (abs_inphase > abs_quadrature) {
      return alpha * abs_inphase + beta * abs_quadrature;
   } else {
      return alpha * abs_quadrature + beta * abs_inphase;
   }
}
/*********************************************************************/
double decibels(double linear)
{
   #define SMALL 1e-20
   if (linear <= SMALL) {
      linear = SMALL;
   }
   return 20.0 * log10(linear);
}
/*********************************************************************/
void test_alpha_beta(char *name, double alpha, double beta,
                     int num_points)
{
   #define PI 3.141592653589793
   int ii;
   double phase, real, imag, err, avg_err, rms_err;
   double peak_err = 0.0;
   double sum_err = 0.0;
   double sum_err_sqrd = 0.0;
   double delta_phase = (2.0 * PI) / num_points;
   for (ii = 0; ii < num_points; ii++) {
      phase = delta_phase * ii;
      real = cos(phase);
      imag = sin(phase);
      err = sqrt(real * real + imag * imag)
            - alpha_beta_mag(alpha, beta, real, imag);
      sum_err += err;
      sum_err_sqrd += err * err;
      err = fabs(err);
      if (err > peak_err) {
         peak_err = err;
      }
   }
   avg_err = sum_err / num_points;
   rms_err = sqrt(sum_err_sqrd / num_points);
   printf("%-16s %14.12lf %14.12lf  %9.6lf %4.1lf %4.1lf\n",
          name, alpha, beta, avg_err, decibels(rms_err),
          decibels(peak_err));
}
/*********************************************************************/
void main(void)
{
   #define NUM_CHECK_POINTS 100000
   typedef struct tagALPHA_BETA {
      char *name;
      double alpha;
      double beta;
   } ALPHA_BETA;
   #define NUM_ALPHA_BETA 16
   const ALPHA_BETA coeff[NUM_ALPHA_BETA] = {
      { "Min RMS Err",      0.947543636291, 0.3924854250920 },
      { "Min Peak Err",     0.960433870103, 0.3978247347593 },
      { "Min RMS w/ Avg=0", 0.948059448969, 0.3926990816987 }, 
      { "1, Min RMS Err",              1.0,     0.323260990 },
      { "1, Min Peak Err",             1.0,     0.335982538 },
      { "1, 1/2",                      1.0,      1.0 / 2.0  },
      { "1, 1/4",                      1.0,      1.0 / 4.0  },
      { "Frerking",                    1.0,            0.4  },
      { "1, 11/32",                    1.0,     11.0 / 32.0 },
      { "1, 3/8",                      1.0,      3.0 / 8.0  },
      { "15/16, 15/32",        15.0 / 16.0,     15.0 / 32.0 },
      { "15/16, 1/2",          15.0 / 16.0,      1.0 / 2.0  },
      { "31/32, 11/32",        31.0 / 32.0,     11.0 / 32.0 },
      { "31/32, 3/8",          31.0 / 32.0,      3.0 / 8.0  },
      { "61/64, 3/8",          61.0 / 64.0,      3.0 / 8.0  },
      { "61/64, 13/32",        61.0 / 64.0,     13.0 / 32.0 }
   };
   int ii;
   printf("\n             Alpha * Max + Beta * Min Magnitude
Estimator\n\n");
   printf("Name                  Alpha           Beta       Avg Err  
RMS   Peak\n");
   printf("                                                 (linear) 
(dB)  (dB)\n");
  
printf("---------------------------------------------------------------------\n");
   for (ii = 0; ii < NUM_ALPHA_BETA; ii++) {
      test_alpha_beta(coeff[ii].name, coeff[ii].alpha, coeff[ii].beta,
                      1024);
   }
}