指数函数实现

这个问题是21icbbs 上的一个网友提出的，我第一反应就是迭代。在很多场合下，比如计算热电流thermal current，或是其他一些使用如下公式的应用：

         y(t) = K * (1 - exp(-t / T)))                                                           (1)

         y(t) = K * exp(-t / T))                                                                   (2)

    where T is the time constant of the system and K is the magnitude.

 

公式（1），(2)不能直接应用在实际计算中，因为 t并没有确切的数值。 对 (1) 微分后得：

　      dy = - K * exp(-t / T)) /T =  ( K - K*(1 – exp(-t/T))) / T

      =  (K – y) / T                                                                               (3)

对 (2) 微分后得：

         dy = - K * exp(-t / T)) /T = - y /T                                               (4)    

 

当 K=0, 则 （3）成为 （4）。假设采样时间为 Ts, 把（3） 转成迭代方程:

        y(k+1) = y(k) + Ts/T * (K – y(k))                                                 (5)

  

只要知道初始值 y(0) 就可计算，与时间 t 无关。如果在 s domain ，转换则更简单：把 s = (1- z)/Ts 或 s = 2/Ts * (1+z)/(1-z) 带入下列公式，可得迭代公式：

       y(s) = K/(s + 1/T) * x(s)                                                                  (6)

 

 

再看看 exp(x) 公式，如果x 是16-bit整数（网友的要求），那么可以构造出一个矩阵，把 exp(x) 变成乘法问题。

 float  Coeff[4][16] = {

              { 1, exp(0x0001), exp(0x0002), …., exp(0x000F)  },

              { 1, exp(0x0010), exp(0x0020), …., exp(0x00F0)  },

              { 1, exp(0x0100), exp(0x0200), …., exp(0x0F00)  },

              { 1, exp(0x1000), exp(0x2000), …., exp(0xF000)  }

  };

 

  float exp(uint16 x)

  {

          return  Coeff[0][x & 0x0F] *

                      Coeff[1][(x>>4) & 0x0F] *    

                      Coeff[2][(x>>8) & 0x0F] *    

                      Coeff[3][(x>>12) & 0x0F];  

   }

 

经过一些网友的提示，我看了一下 CORDIC 算法，上面的方法是把exp(x) 变成乘法，如果能精心挑选x分解项，把 Coeff 变成2^N，那么意味着在 FPGA 以及 DSP 上只要使用加减和左移右移就能达到目的，这就是CORDIC 的精髓。X 可分解为

     x = k0 + k1 + … + kn ,

exp(x) = exp(k0) * exp(k1) * … * exp(kn)

           =  2^N0 * 2^N1 * … * 2 * (1+ ½) * (1+ 1/4) * (1+ 1/8) …

 

注意当xn 趋向0 时，exp(xn) 趋向1, 所以当 0 < x < ln(2) 时，无法单纯的左移右移, 这时可以让 exp(kn) = 1 + 1/(2^N), 则运算变成

     y*(1 + 1/(2^N)) = y + (y >> N)

 

验证程序如下：

 

#include "stdafx.h"

#include "math.h"

 

#define LN2 0.69314718055994530941723212145818f

 

float C1[] =

{

            128*LN2, 64*LN2, 32*LN2, 16*LN2, 8*LN2, 4*LN2, 2*LN2, LN2

};

 

float C2[] =

{

            0.405465108f, 0.223143551f, 0.117783036f, 0.060624622f,

            0.030771659f, 0.015504187f, 0.00778214f, 0.00389864f,

            0.00195122f, 0.000976086f, 0.000488162f, 0.000244111f

};

 

 

float Exp(float x)

{

      float y = 1.0;

      int i;

      for (i=0; i<8; i++) {

            if (x >= C1[i]) {

                  x -= C1[i];

                  y *= 1L <<(1<<(7-i));

            }

      }

      for (i=0; i<12; i++) {

            if (x >= C2[i]) {

                  x -= C2[i];

                  y += y / (1<<(i+1));

            }

      }

      return y;   

}

 

 

int _tmain(int argc, _TCHAR* argv[])

{

            float y1 = Exp(5.0);

            float y2 = (float) exp(5.0);

 

            printf("%4.5f   %4.5f", y1, y2);

            return 0;

}

 

 

用定点计算的程序

 

#include "stdafx.h"

#include "math.h"

 

#define uint32 unsigned int

#define int32 int

#define ONE 65536

 

uint32 C1[] =

{

            726817, 363408, 181704, 90852, 45426

};

 

uint32 C2[] =

{

            26572, 14623, 7719, 3973, 2016, 1016,

            510, 256, 128, 64, 32, 16, 8, 4, 2

};

 

// Convert: uint32 = float * 65536

// ONE = 0x00010000 <---> 1.0 in float

// ln(65536) = 11.090,

// so x must be less than 11.090*ONE = 726817

uint32 Exp(uint32 x)

{

            uint32 y = ONE;

            if (x >= C1[0]) { x -= C1[0]; y <<= 16; }

            if (x >= C1[1]) { x -= C1[1]; y <<= 8; }

            if (x >= C1[2]) { x -= C1[2]; y <<= 4; }

            if (x >= C1[3]) { x -= C1[3]; y <<= 2; }

            if (x >= C1[4]) { x -= C1[4]; y <<= 1; }

 

            for (int i=0; i<16; i++) {

                        if (x >= C2[i]) { x -= C2[i]; y += (y >> (i+1)); }

            }

            return y;

}

 

int _tmain(int argc, _TCHAR* argv[])

{

            float y1 = (float) Exp((uint32) 5.0f * ONE) / ONE;

            float y2 = (float) exp(5.0);

 

            printf("%4.5f   %4.5f", y1, y2);

            return 0;

}

       

 

结果：

 

148.40210     148.41316

148.41339     148.41316