一种复数求模的实用方法

一种复数求模的实用方法

        在DSP中，经常遇到复数求模的问题。对于复数X

                        X = R+j*I   (1)

其幅度magX为:

                     magX = sqrt(R^2+I^2)  (2)

由上式可知，复数求模的问题可看做是一个实数求根的问题。对这个问题有一些成熟算法，如博文所述：http://blog.youkuaiyun.com/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 <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

*

*********************************************************************/

/********************************************************************/

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);

   }

}