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://scopeplot.com/dg/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);
}
}