关于嵌入式性能调优

  本来对嵌入式性能调优方面的一些文章还有些不以为然，但经过自己亲身体验之后，我不得不承认那些都是无比正确的。
  对于性能问题，我一向的观点是，首先优化算法，能够查表的情况下就查表（如三角函数），能展开写的就不循环甚至递归。然后是尽量少的拷贝和分配内存，数据能共享的就共享，能一次分配的内存就一次分配。最后才是代码细节的优化。代码细节的优化大概有以下几点：
  1、能预先算好的就预先算好。比如旋转变换的sin(angle)和cos(angle)就可以预先算好，而不必每个点都重新算一遍。
  2、干掉所有浮点运算。当然，这是不可能的，我们只能做到在最花CPU时间的地方见不到浮点运算，这往往要结合第一点来实现，我们可以预先算好需要的参数，假如它是浮点数，那么我们可以将其乘以2^n，这个n取决于需要的精度，然后取整。在每个需要乘以这个参数的地方加上>>n。
  3、干掉所有除法。我们往往可以将其转化为浮点乘法，然后应用第二点。
  4、在可能的情况下，用移位代替乘除法。
  5、赋值、比较等时候，数据类型要相符，再不济也要加上强制类型转换。
  6、多用宏函数或内联函数。
  下面是一些高效的算法：
一、平方根

1. /* fsqrt.c
2.  *
3.  * A fast square root program adapted from the code of
4.  * Paul Lalonde and Robert Dawson in Graphics Gems I.
5.  * The format of IEEE double precision floating point numbers is:
6.  *
7.  * SEEEEEEEEEEEMMMM MMMMMMMMMMMMMMMM MMMMMMMMMMMMMMMM MMMMMMMMMMMMMMMM
8.  *
9.  * S = Sign bit for whole number
10.  * E = Exponent bit (exponent in excess 1023 form)
11.  * M = Mantissa bit
12.  */
14. #include <stdio.h>
15. #include <math.h>
17. /* MOST_SIG_OFFSET gives the (int *) offset from the address of the double
18.  * to the part of the number containing the sign and exponent.
19.  * You will need to find the relevant offset for your architecture.
20.  */
22. #define MOST_SIG_OFFSET 1
24. /* SQRT_TAB_SIZE - the size of the lookup table - must be a power of four.
25.  */
27. #define SQRT_TAB_SIZE 16384
29. /* MANT_SHIFTS is the number of shifts to move mantissa into position.
30.  * If you quadruple the table size subtract two from this constant,
31.  * if you quarter the table size then add two.
32.  * Valid values are: (16384, 7) (4096, 9) (1024, 11) (256, 13)
33.  */
35. #define MANT_SHIFTS   7
37. #define EXP_BIAS   1023       /* Exponents are always positive     */
38. #define EXP_SHIFTS 20         /* Shifs exponent to least sig. bits */
39. #define EXP_LSB    0x00100000 /* 1 << EXP_SHIFTS                   */
40. #define MANT_MASK  0x000FFFFF /* Mask to extract mantissa          */
42. int        sqrt_tab[SQRT_TAB_SIZE];
44. void
45. init_sqrt_tab()
46. {
47.         int           i;
48.         double        f;
49.         unsigned int  *fi = (unsigned int *) &f + MOST_SIG_OFFSET;
50.         
51.         for (i = 0; i < SQRT_TAB_SIZE/2; i++)
52.         {
53.                 f = 0; /* Clears least sig part */
54.                 *fi = (i << MANT_SHIFTS) | (EXP_BIAS << EXP_SHIFTS);
55.                 f = sqrt(f);
56.                 sqrt_tab[i] = *fi & MANT_MASK;
58.                 f = 0; /* Clears least sig part */
59.                 *fi = (i << MANT_SHIFTS) | ((EXP_BIAS + 1) << EXP_SHIFTS);
60.                 f = sqrt(f);
61.                 sqrt_tab[i + SQRT_TAB_SIZE/2] = *fi & MANT_MASK;
62.         }
63. }
65. double
66. fsqrt(f)
67. double f;
68. {
69.         unsigned int e;
70.         unsigned int   *fi = (unsigned int *) &f + MOST_SIG_OFFSET;
72.         if (f == 0.0) return(0.0);
73.         e = (*fi >> EXP_SHIFTS) - EXP_BIAS;
74.         *fi &= MANT_MASK;
75.         if (e & 1)
76.                 *fi |= EXP_LSB;
77.         e >>= 1;
78.         *fi = (sqrt_tab[*fi >> MANT_SHIFTS]) |
79.               ((e + EXP_BIAS) << EXP_SHIFTS);
80.         return(f);
81. }
83. void
84. dump_sqrt_tab()
85. {
86.         int        i, nl = 0;
88.         printf("unsigned int sqrt_tab[] = {/n");
89.         for (i = 0; i < SQRT_TAB_SIZE-1; i++)
90.         {
91.                 printf("0x%x,", sqrt_tab[i]);
92.                 nl++;
93.                 if (nl > 8) { nl = 0; putchar('/n'); }
94.         }
95.         printf("0x%x/n", sqrt_tab[SQRT_TAB_SIZE-1]);
96.         printf("};/n");
97. }

二、距离计算

1. //
2. // 这个函数计算从(0,0)到(x,y)的距离，相对误差为3.5%
3. //
4. int FastDistance2D(int x, int y)
5. {
6.     x = abs(x);
7.     y = abs(y);
8.     int mn = MIN(x,y);
9.     return(x+y-(mn>>1)-(mn>>2)+(mn>>4));
10. }
11. //
12. // 该函数计算(0,0,0)到(x,y,z)的距离，相对误差为8%
13. //
14. float FastDistance3D(float fx, float fy, float fz)
15. {
16.     int temp;
17.     int x,y,z;
18.     // 确保所有的值为正
19.     x = int(fabs(fx) * 1024);
20.     y = int(fabs(fy) * 1024);
21.     z = int(fabs(fz) * 1024);
22.     // 排序
23.     if (y < x) SWAP(x,y,temp)
24.     if (z < y) SWAP(y,z,temp)
25.     if (y < x) SWAP(x,y,temp)
26.     int dist = (z + 11 * (y >> 5) + (x >> 2) );
27.     return((float)(dist >> 10));
28. }

三、随机数发生

1. /*
2. A 32 bit random number generator. An implementation in C of the algorithm given by
3. Knuth, the art of computer programming, vol. 2, pp. 26-27. We use e=32, so 
4. we have to evaluate y(n) = y(n - 24) + y(n - 55) mod 2^32, which is implicitly
5. done by unsigned arithmetic.
6. */
8. unsigned int Random32(void) {
9.   /*
10.   random numbers from Mathematica 2.0.
11.   SeedRandom = 1;
12.   Table[Random[Integer, {0, 2^32 - 1}]
13.   */
14.   static const unsigned long x[55] = {
15.     1410651636UL, 3012776752UL, 3497475623UL, 2892145026UL, 1571949714UL,
16.     3253082284UL, 3489895018UL, 387949491UL, 2597396737UL, 1981903553UL,
17.     3160251843UL, 129444464UL, 1851443344UL, 4156445905UL, 224604922UL,
18.     1455067070UL, 3953493484UL, 1460937157UL, 2528362617UL, 317430674UL, 
19.     3229354360UL, 117491133UL, 832845075UL, 1961600170UL, 1321557429UL,
20.     747750121UL, 545747446UL, 810476036UL, 503334515UL, 4088144633UL,
21.     2824216555UL, 3738252341UL, 3493754131UL, 3672533954UL, 29494241UL,
22.     1180928407UL, 4213624418UL, 33062851UL, 3221315737UL, 1145213552UL,
23.     2957984897UL, 4078668503UL, 2262661702UL, 65478801UL, 2527208841UL,
24.     1960622036UL, 315685891UL, 1196037864UL, 804614524UL, 1421733266UL,
25.     2017105031UL, 3882325900UL, 810735053UL, 384606609UL, 2393861397UL };
26.   static int init = 1;
27.   static unsigned long y[55];
28.   static int j, k;
29.   unsigned long ul;
30.   
31.   if (init)
32.   {
33.     int i;
34.     
35.     init = 0;
36.     for (i = 0; i < 55; i++) y[i] = x[i];
37.     j = 24 - 1;
38.     k = 55 - 1;
39.   }
40.   
41.   ul = (y[k] += y[j]);
42.   if (--j < 0) j = 55 - 1;
43.   if (--k < 0) k = 55 - 1;
44.   return((unsigned int)ul);
45. }

四、FFT变换

1. /*
2. This computes an in-place complex-to-complex FFT
3. x and y are the real and imaginary arrays of n=2^m points. 
4. o(n)=n*log2(n) 
5. dir =  1 gives forward transform 
6. dir = -1 gives reverse transform  
7. FFT algorithm by Cooley and Tukey, 1965
8. */
11. bool _stdcall FFT(int dir,int m,double *x,double *y)
12. {
13.     long nn,i,i1,j,k,i2,l,l1,l2; 
14.     double c1,c2,tx,ty,t1,t2,u1,u2,z; 
15.     
16.     /* Calculate the number of points */
17.     nn = 1<<m; 
18.     
19.     /* Do the bit reversal */ 
20.     i2 = nn >> 1; 
21.     j = 0; 
22.     for (i=0;i<nn-1;i++) { 
23.         if (i < j) { 
24.             tx = x[i]; 
25.             ty = y[i]; 
26.             x[i] = x[j]; 
27.             y[i] = y[j]; 
28.             x[j] = tx; 
29.             y[j] = ty; 
30.         } 
31.         k = i2; 
32.         while (k <= j) { 
33.             j -= k; 
34.             k >>= 1; 
35.         } 
36.         j += k; 
37.     } 
38.     
39.     /* Compute the FFT */ 
40.     c1 = -1.0; 
41.     c2 = 0.0; 
42.     l2 = 1; 
43.     for (l=0;l<m;l++) { 
44.         l1 = l2; 
45.         l2 <<= 1; 
46.         u1 = 1.0; 
47.         u2 = 0.0; 
48.         for (j=0;j<l1;j++) { 
49.             for (i=j;i<nn;i+=l2) { 
50.                 i1 = i + l1; 
51.                 t1 = u1 * x[i1] - u2 * y[i1]; 
52.                 t2 = u1 * y[i1] + u2 * x[i1]; 
53.                 x[i1] = x[i] - t1; 
54.                 y[i1] = y[i] - t2; 
55.                 x[i] += t1; 
56.                 y[i] += t2; 
57.             } 
58.             z =  u1 * c1 - u2 * c2; 
59.             u2 = u1 * c2 + u2 * c1; 
60.             u1 = z; 
61.         } 
62.         c2 = sqrt((1.0 - c1) / 2.0); 
63.         if (dir == 1) 
64.             c2 = -c2; 
65.         c1 = sqrt((1.0 + c1) / 2.0); 
66.     } 
67.     
68.     /* Scaling for forward transform */ 
69.     if (dir == 1) { 
70.         for (i=0;i<nn;i++) { 
71.             x[i] /= (double)nn; 
72.             y[i] /= (double)nn; 
73.         }
74.     } 
75.     
76.     return true; 
77. }

  最后，不得不提：STL真的很低效。