一种复数求模的实用方法
在DSP中，经常遇到复数求模的问题。对于复数X
X = R+j*I (1)
其幅度magX为:
magX = sqrt(R^2+I^2) (2)
由上式可知，复数求模的问题可看做是一个实数求根的问题。对这个问题有一些成熟算法，如博文所述：
http://blog.csdn.net/deepdsp/article/details/7539823
这里，针对复数求模这种特殊情况，在精度要求不太高的情况下，可用如下近似方法计算：
MagX = Alpha * max(|R|, |I|) + Beta * min(|R|, |I|) (3)
其中Alpha和Beta为系数，根据不同的准则，如均方根误差最小，峰值误差最小等，这两个系数的取值有所不同。具体数值参见下表。最简单的情况，选取Alpha=1，Beta=0.25，也能得到很不错的结果。
=====================================================================
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
=====================================================================
附录：C实现代码
/*(代码来源于网络，版权归原作者所有)*/
#include
#include
/*********************************************************************
*
*
* 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
*
*********************************************************************/
/********************************************************************/
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
",
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("
Alpha * Max + Beta * Min Magnitude
Estimator
");
printf("Name Alpha Beta Avg Err
RMS Peak
");
printf(" (linear)
(dB) (dB)
");
printf("---------------------------------------------------------------------
");
for (ii = 0; ii < NUM_ALPHA_BETA; ii++) {
test_alpha_beta(coeff[ii].name, coeff[ii].alpha, coeff[ii].beta,
1024);
}
}
{var%20f='http://v.t.sina.com.cn/share/share.php?appkey=1515056452',u=z||d.location,p=['&url=',e(u),'&title=',e(t||d.title),'&source=',e(r),'&sourceUrl=',e(l),'&content=',c||'gb2312','&pic=',e(p||'')].join('');function%20a(){if(!window.open([f,p].join(''),'mb',['toolbar=0,status=0,resizable=1,width=440,height=430,left=',(s.width-440)/2,',top=',(s.height-430)/2].join('')))u.href=[f,p].join('');};if(/Firefox/.test(navigator.userAgent))setTimeout(a,0);else%20a();})(screen,document,encodeURIComponent,'','','https://www.xiaopingtou.cn//data/attach/topic/topicKPo7gB.jpg', '推荐 Dorsey_c 的文章《一种复数求模的实用方法》','https://www.xiaopingtou.net/article-51885.html','页面编码gb2312|utf-8默认gb2312'));){kind=link}