calculate.c 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616
  1. /******************************************************************************
  2. 版权所有:
  3. 文件名称: calculate.c
  4. 文件版本: 01.01
  5. 创建作者: sunxi
  6. 创建日期: 2008-08-25
  7. 功能说明: 算术计算公式。
  8. 其它说明: Q16定点数定义----如果一个32bits数,高16代表整数部分,低16位代表小数部分,
  9. 我们就将它称为Q16定点数(Qxx代表xx位小数).
  10. 修改记录:
  11. */
  12. /*------------------------------- 头文件 --------------------------------------
  13. */
  14. #ifndef WIN32 // VC6.0 中定义的宏
  15. #include "bspconfig.h"
  16. #include "calculate.h"
  17. #define __int64 long long
  18. #endif
  19. #include "rt.h"
  20. /*------------------------------- 宏定义 --------------------------------------
  21. */
  22. #if CFG_BSP_DEBUG
  23. #define _CAL_DEBUG
  24. #endif
  25. /*------------------------------ 全局变量 -------------------------------------
  26. */
  27. const long asin_table[257] =
  28. {
  29. 0x00000000,
  30. 0x0000394B,
  31. 0x00007297,
  32. 0x0000ABE4,
  33. 0x0000E531,
  34. 0x00011E7F,
  35. 0x000157CE,
  36. 0x0001911E,
  37. 0x0001CA70,
  38. 0x000203C4,
  39. 0x00023D1A,
  40. 0x00027672,
  41. 0x0002AFCD,
  42. 0x0002E92A,
  43. 0x0003228A,
  44. 0x00035BED,
  45. 0x00039554,
  46. 0x0003CEBE,
  47. 0x0004082C,
  48. 0x0004419F,
  49. 0x00047B15,
  50. 0x0004B490,
  51. 0x0004EE10,
  52. 0x00052795,
  53. 0x0005611E,
  54. 0x00059AAE,
  55. 0x0005D443,
  56. 0x00060DDE,
  57. 0x0006477F,
  58. 0x00068126,
  59. 0x0006BAD5,
  60. 0x0006F489,
  61. 0x00072E46,
  62. 0x00076809,
  63. 0x0007A1D4,
  64. 0x0007DBA7,
  65. 0x00081581,
  66. 0x00084F65,
  67. 0x00088950,
  68. 0x0008C345,
  69. 0x0008FD42,
  70. 0x00093749,
  71. 0x00097159,
  72. 0x0009AB74,
  73. 0x0009E598,
  74. 0x000A1FC6,
  75. 0x000A59FF,
  76. 0x000A9443,
  77. 0x000ACE92,
  78. 0x000B08EC,
  79. 0x000B4352,
  80. 0x000B7DC4,
  81. 0x000BB842,
  82. 0x000BF2CC,
  83. 0x000C2D63,
  84. 0x000C6807,
  85. 0x000CA2B8,
  86. 0x000CDD76,
  87. 0x000D1843,
  88. 0x000D531D,
  89. 0x000D8E06,
  90. 0x000DC8FD,
  91. 0x000E0403,
  92. 0x000E3F19,
  93. 0x000E7A3E,
  94. 0x000EB572,
  95. 0x000EF0B7,
  96. 0x000F2C0C,
  97. 0x000F6772,
  98. 0x000FA2E9,
  99. 0x000FDE71,
  100. 0x00101A0B,
  101. 0x001055B6,
  102. 0x00109174,
  103. 0x0010CD45,
  104. 0x00110928,
  105. 0x0011451F,
  106. 0x00118129,
  107. 0x0011BD46,
  108. 0x0011F979,
  109. 0x001235BF,
  110. 0x0012721B,
  111. 0x0012AE8C,
  112. 0x0012EB12,
  113. 0x001327AE,
  114. 0x00136461,
  115. 0x0013A12A,
  116. 0x0013DE0A,
  117. 0x00141B02,
  118. 0x00145812,
  119. 0x00149539,
  120. 0x0014D27A,
  121. 0x00150FD3,
  122. 0x00154D45,
  123. 0x00158AD2,
  124. 0x0015C878,
  125. 0x00160639,
  126. 0x00164415,
  127. 0x0016820C,
  128. 0x0016C01F,
  129. 0x0016FE4F,
  130. 0x00173C9B,
  131. 0x00177B04,
  132. 0x0017B98B,
  133. 0x0017F830,
  134. 0x001836F3,
  135. 0x001875D5,
  136. 0x0018B4D7,
  137. 0x0018F3F9,
  138. 0x0019333B,
  139. 0x0019729E,
  140. 0x0019B222,
  141. 0x0019F1C9,
  142. 0x001A3192,
  143. 0x001A717E,
  144. 0x001AB18D,
  145. 0x001AF1C1,
  146. 0x001B3219,
  147. 0x001B7297,
  148. 0x001BB33A,
  149. 0x001BF404,
  150. 0x001C34F4,
  151. 0x001C760C,
  152. 0x001CB74D,
  153. 0x001CF8B6,
  154. 0x001D3A48,
  155. 0x001D7C05,
  156. 0x001DBDED,
  157. 0x001E0000,
  158. 0x001E423F,
  159. 0x001E84AA,
  160. 0x001EC744,
  161. 0x001F0A0B,
  162. 0x001F4D02,
  163. 0x001F9028,
  164. 0x001FD37F,
  165. 0x00201707,
  166. 0x00205AC1,
  167. 0x00209EAD,
  168. 0x0020E2CE,
  169. 0x00212723,
  170. 0x00216BAE,
  171. 0x0021B06F,
  172. 0x0021F566,
  173. 0x00223A97,
  174. 0x00228000,
  175. 0x0022C5A3,
  176. 0x00230B81,
  177. 0x0023519B,
  178. 0x002397F2,
  179. 0x0023DE87,
  180. 0x0024255B,
  181. 0x00246C6F,
  182. 0x0024B3C4,
  183. 0x0024FB5C,
  184. 0x00254338,
  185. 0x00258B58,
  186. 0x0025D3BF,
  187. 0x00261C6D,
  188. 0x00266563,
  189. 0x0026AEA3,
  190. 0x0026F82F,
  191. 0x00274207,
  192. 0x00278C2E,
  193. 0x0027D6A4,
  194. 0x0028216B,
  195. 0x00286C84,
  196. 0x0028B7F2,
  197. 0x002903B6,
  198. 0x00294FD1,
  199. 0x00299C46,
  200. 0x0029E915,
  201. 0x002A3641,
  202. 0x002A83CD,
  203. 0x002AD1B8,
  204. 0x002B2007,
  205. 0x002B6EBA,
  206. 0x002BBDD4,
  207. 0x002C0D58,
  208. 0x002C5D46,
  209. 0x002CADA2,
  210. 0x002CFE6E,
  211. 0x002D4FAD,
  212. 0x002DA161,
  213. 0x002DF38D,
  214. 0x002E4633,
  215. 0x002E9957,
  216. 0x002EECFB,
  217. 0x002F4122,
  218. 0x002F95D0,
  219. 0x002FEB08,
  220. 0x003040CD,
  221. 0x00309723,
  222. 0x0030EE0D,
  223. 0x00314590,
  224. 0x00319DAF,
  225. 0x0031F66E,
  226. 0x00324FD3,
  227. 0x0032A9E0,
  228. 0x0033049C,
  229. 0x0033600B,
  230. 0x0033BC31,
  231. 0x00341915,
  232. 0x003476BC,
  233. 0x0034D52C,
  234. 0x0035346B,
  235. 0x00359480,
  236. 0x0035F571,
  237. 0x00365746,
  238. 0x0036BA05,
  239. 0x00371DB8,
  240. 0x00378265,
  241. 0x0037E817,
  242. 0x00384ED6,
  243. 0x0038B6AD,
  244. 0x00391FA5,
  245. 0x003989CB,
  246. 0x0039F529,
  247. 0x003A61CC,
  248. 0x003ACFC2,
  249. 0x003B3F19,
  250. 0x003BAFE0,
  251. 0x003C2228,
  252. 0x003C9603,
  253. 0x003D0B83,
  254. 0x003D82BD,
  255. 0x003DFBC7,
  256. 0x003E76BA,
  257. 0x003EF3AF,
  258. 0x003F72C2,
  259. 0x003FF414,
  260. 0x004077C6,
  261. 0x0040FDFE,
  262. 0x004186E6,
  263. 0x004212AC,
  264. 0x0042A183,
  265. 0x004333A6,
  266. 0x0043C955,
  267. 0x004462DB,
  268. 0x0045008B,
  269. 0x0045A2C8,
  270. 0x00464A00,
  271. 0x0046F6B9,
  272. 0x0047A98D,
  273. 0x00486337,
  274. 0x00492498,
  275. 0x0049EEC8,
  276. 0x004AC325,
  277. 0x004BA374,
  278. 0x004C9209,
  279. 0x004D921C,
  280. 0x004EA84C,
  281. 0x004FDBB3,
  282. 0x00513846,
  283. 0x0052D556,
  284. 0x0054EF1F,
  285. 0x005A0000,
  286. };
  287. /*------------------------------ 函数声明 -------------------------------------
  288. */
  289. long _cal_asin(long x);
  290. /*------------------------------ 外部函数 -------------------------------------
  291. */
  292. /******************************************************************************
  293. 函数名称: cal_complex_magnitude
  294. 函数版本: 01.01
  295. 创建作者: sunxi
  296. 创建日期: 2008-08-26
  297. 函数说明: 得到复数的幅值(sqrt(real*real + imag*imag))。算法如下:
  298. l=max{|real|,|imag|},s=min{|real|,|imag|}
  299. m = l + 5*s*s/(9*l+3*s);
  300. 此公式的实数算法精度(最大相对误差值)为0.1735%.
  301. 参数说明:
  302. real: 复数的实部
  303. imag: 复数的虚部
  304. 返回值: 复数的幅值.
  305. 修改记录:
  306. */
  307. unsigned long cal_complex_magnitude(long real, long imag)
  308. {
  309. unsigned long m = 0, l, t, s;
  310. if (real == 0 && imag == 0)
  311. {
  312. return 0;
  313. }
  314. if (real < 0)
  315. {
  316. real = -real;
  317. }
  318. if (imag < 0)
  319. {
  320. imag = -imag;
  321. }
  322. if (real < imag)
  323. {
  324. l = (unsigned long)imag;
  325. s = (unsigned long)real;
  326. }
  327. else
  328. {
  329. l = (unsigned long)real;
  330. s = (unsigned long)imag;
  331. }
  332. // 计算结果,并进行四舍五入.
  333. t = 9 * l + 3 * s;
  334. if (t)
  335. {
  336. m = (unsigned long)(l + (5 * (__int64)s * s + (t >> 1)) / t);
  337. }
  338. return m;
  339. }
  340. /******************************************************************************
  341. 函数名称: cal_complex_phase_cos
  342. 函数版本: 01.01
  343. 创建作者: sunxi
  344. 创建日期: 2008-08-26
  345. 函数说明: 得到复数极坐标表示法的相位角的余弦
  346. 参数说明:
  347. real: 复数的实部
  348. imag: 复数的虚部
  349. 返回值: 复数相位角的余弦值(Q16定点数)
  350. 修改记录:
  351. */
  352. long cal_complex_phase_cos(long real, long imag)
  353. {
  354. unsigned long m;
  355. m = cal_complex_magnitude(real, imag);
  356. if (m == 0)
  357. {
  358. return 0;
  359. }
  360. return (long)((((__int64)real << 16) + (m >> 1)) / m);
  361. }
  362. /******************************************************************************
  363. 函数名称: cal_complex_phase
  364. 函数版本: 01.01
  365. 创建作者: sunxi
  366. 创建日期: 2008-08-26
  367. 函数说明: 得到复数的相位角.
  368. 参数说明:
  369. real: 复数的实部
  370. imag: 复数的虚部
  371. 返回值: 复数的相位角(-180~180,Q16定点数)
  372. 修改记录:
  373. */
  374. #define ANGLE_90 (90 << 16)
  375. #define ANGLE_180 (180 << 16)
  376. long cal_complex_phase(long real, long imag)
  377. {
  378. long real_abs, imag_abs;
  379. long angle;
  380. // 得到绝对值
  381. real_abs = real >= 0 ? real : -real;
  382. imag_abs = imag >= 0 ? imag : -imag;
  383. // 对于sin函数来说,0~45度的线性要好于45~90度的线性。
  384. // 所以我们只使用0~45度。
  385. if (real_abs >= imag_abs)
  386. {
  387. //<=45度,用sin去算
  388. angle = _cal_asin(cal_complex_phase_cos(imag_abs, real_abs));
  389. }
  390. else
  391. {
  392. //>45度,用cos去算
  393. angle = _cal_asin(cal_complex_phase_cos(real_abs, imag_abs));
  394. angle = ANGLE_90 - angle;
  395. }
  396. // 象限处理
  397. if (real >= 0)
  398. {
  399. if (imag >= 0)
  400. {
  401. return angle; // 第一象限
  402. }
  403. else
  404. {
  405. return (0 - angle); // 第四象限
  406. }
  407. }
  408. else
  409. {
  410. if (imag >= 0)
  411. {
  412. return ANGLE_180 - angle; // 第二象限
  413. }
  414. else
  415. {
  416. return (-ANGLE_180 + angle); // 第三象限
  417. }
  418. }
  419. }
  420. /******************************************************************************
  421. 函数名称: cal_q16_multi
  422. 函数版本: 01.01
  423. 创建作者: sunxi
  424. 创建日期: 2008-08-26
  425. 函数说明: 两个Q16定点数相乘,得到一个Q16数.计算公式为:
  426. (a*b)>>16
  427. 参数说明:
  428. a: Q16定点数1
  429. b: Q16定点数2
  430. 返回值: 乘积(Q16定点数).注:相乘后的结果应不超过48位,否则出现溢出错误
  431. 修改记录:
  432. */
  433. long cal_q16_multi(long a, long b)
  434. {
  435. __int64 ll;
  436. ll = (__int64)a * b;
  437. return (long)(ll >> 16);
  438. }
  439. /******************************************************************************
  440. 函数名称: cal_multi_divide
  441. 函数版本: 01.01
  442. 创建作者: sunxi
  443. 创建日期: 2008-08-26
  444. 函数说明: 使用64位算法计算a*b/c.有符号算法.注意:结果仍有可能溢出.
  445. 参数说明:
  446. a: 被乘数
  447. b: 乘数
  448. c: 除数
  449. 返回值: 结果
  450. 修改记录:
  451. */
  452. long cal_multi_divide(long a, long b, long c)
  453. {
  454. __int64 ll;
  455. #ifdef __KERNEL__
  456. ll = a * b;
  457. div_s64(ll, c);
  458. #else
  459. ll = (__int64)a * b / c;
  460. #endif
  461. return (long)ll;
  462. }
  463. /******************************************************************************
  464. 函数名称: cal_multi_divide_unsigned
  465. 函数版本: 01.01
  466. 创建作者: sunxi
  467. 创建日期: 2008-08-26
  468. 函数说明: 使用64位算法计算a*b/c.无符号算法.注意:结果仍有可能溢出.
  469. 参数说明:
  470. a: 被乘数
  471. b: 乘数
  472. c: 除数
  473. 返回值: 结果
  474. 修改记录:
  475. */
  476. unsigned long cal_multi_divide_unsigned(unsigned long a, unsigned long b, unsigned long c)
  477. {
  478. unsigned __int64 ll;
  479. ll = (unsigned __int64)a * b / c;
  480. return (unsigned long)ll;
  481. }
  482. /*------------------------------ 内部函数 -------------------------------------
  483. */
  484. /******************************************************************************
  485. 函数名称: _cal_asin
  486. 函数版本: 01.01
  487. 创建作者: sunxi
  488. 创建日期: 2008-08-26
  489. 函数说明: 反sin函数
  490. 参数说明:
  491. x: 反sin函数输入值,Q16定点数,0<=x<=1。
  492. 返回值: 结果y,Q16定点数,0<=y<=90.
  493. 修改记录:
  494. */
  495. long _cal_asin(long x)
  496. {
  497. int high, low;
  498. // x>=1,返回90度。
  499. if (x >= 0x10000)
  500. {
  501. return asin_table[256];
  502. }
  503. high = (x >> 8) & 0xff;
  504. low = x & 0xff;
  505. return ((asin_table[high + 1] - asin_table[high]) * low + (256 >> 1)) / 256 + asin_table[high];
  506. }
  507. /*------------------------------ 测试函数 -------------------------------------
  508. */
  509. #ifndef WIN32 // VC6.0 中定义的宏
  510. #ifdef _CAL_DEBUG
  511. #include "rt.h"
  512. int cal_test(void)
  513. {
  514. long angle, a, b, ans;
  515. // int i1,i2;
  516. rt_printf("cal_test begin!\r\n");
  517. a = 0x77777;
  518. b = 0x66666;
  519. ans = cal_q16_multi(a, b);
  520. rt_printf("cal_q16_multi:%d*%d=%d\r\n", a, b, ans);
  521. a = 234;
  522. b = 567;
  523. angle = cal_complex_phase(a, b);
  524. rt_printf("angle=%d.%04d(real=%d,imag=%d)\r\n", angle >> 16, (angle & 0xffff) * 10000 / 0x10000, a, b);
  525. rt_printf("cal_test end!\r\n\r\n");
  526. return 0;
  527. }
  528. #endif
  529. #else
  530. #include "stdio.h"
  531. #include "math.h"
  532. int cal_test_vc(void)
  533. {
  534. long i, j, angle;
  535. double angle_db1, angle_db2, diff, diff_max;
  536. angle = cal_complex_phase(136, -8010);
  537. diff_max = 0;
  538. for (i = 100; i < 65536; i++)
  539. {
  540. for (j = 100; j < 65536; j++)
  541. {
  542. angle = cal_complex_phase(i, j);
  543. angle_db1 = (double)angle / 65536;
  544. angle_db2 = atan((double)j / i);
  545. angle_db2 = angle_db2 * 180 / 3.1416;
  546. diff = angle_db2 - angle_db1;
  547. if (diff < 0)
  548. {
  549. diff = -diff;
  550. }
  551. if (diff > diff_max)
  552. {
  553. diff_max = diff;
  554. rt_printf("cal_test_vc:i=%d,j=%d,diff_max=%f\n", i, j, diff);
  555. }
  556. }
  557. }
  558. return 0;
  559. }
  560. #endif