Guest User

4-th order fixed-point sine approx

a guest
Apr 24th, 2017
185
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // 4th-order sine approximation, in 64-bit Q17 -> Q16
  2.  
  3. const double MY_PI = 3.141592653589793238463;
  4.  
  5. typedef int64_t Angl;
  6. int64_t engineer_math_sin(Angl input)
  7. {
  8. // A sine approximation via a fourth-order cosine approx.
  9. // @param x angle (with 2^17 units/circle)
  10. // @return Sine value (Q16)
  11.  
  12. int64_t c, x, x2, y;
  13. static const int64_t
  14. qB = 64,
  15. qN = 15,
  16. qA = 16,
  17. qP = 24,
  18. B = (2<<qP) - int64_t((1<<(qP-2)) * MY_PI), // (2 - pi/4)<<24
  19. C = (1<<qP) - int64_t((1<<(qP-2)) * MY_PI); // (1 - pi/4)<<24
  20.  
  21. x = input;
  22.  
  23. c = x << (qB-2 - qN); // Semi-circle info into carry. 60, 62, 62
  24. x -= 1 << qN; // cosine -> sine calc
  25.  
  26. x = x << (qB-1 - qN); // Mask with PI
  27. x = x >> (qB-1 - qN); // Note: SIGNED shift! (to qN)
  28. x2 = x * x >> (2*qN - qP); // x=x^2 To Q24
  29.  
  30. y = B - (x2 * C >> qP); // B - x^2*C
  31. y = (1 << qA) - (x2*y >> (2*qP-qA)); // A - x^2*(B-x^2*C)
  32.  
  33. return c >= 0 ? y : -y;
  34. }
  35.  
  36. void test_sine()
  37. {
  38. auto sine = engineer_math_sin;
  39.  
  40. struct { Angl a; int64_t check, result; } tests[] = {
  41. {-4<<14, 0<<16 }, // sin(-pi*1) = 0
  42. {-3<<14, -(1<<16)/sqrt(2) }, // sin(-pi*3/4) = -1/sqrt(2)
  43. {-2<<14, -(1<<16) }, // sin(-pi*2/4) = -1
  44. {-1<<14, -(1<<16)/sqrt(2) }, // sin(-pi*1/4) = -1/sqrt(2)
  45.  
  46. { 0<<14, (0<<16) }, // sin(0) = 0
  47. { 1<<14, (1<<16)/sqrt(2) }, // sin(pi*1/4) = 1/sqrt(2)
  48. { 2<<14, 1<<16 }, // sin(pi*2/4) = 1
  49. { 3<<14, (1<<16)/sqrt(2) }, // sin(pi*3/4) = 1/sqrt(2)
  50. { 4<<14, 0<<16 }, // sin(pi*1) = 0
  51. { 5<<14, -(1<<16)/sqrt(2) }, // sin(pi*5/4) = -1/sqrt(2)
  52. { 6<<14, -(1<<16) }, // sin(pi*5/4) = -1/sqrt(2)
  53. { 7<<14, -(1<<16)/sqrt(2) }, // sin(pi*5/4) = -1/sqrt(2)
  54. { 8<<14, (0<<16) }, // sin(pi*5/4) = -1/sqrt(2)
  55. };
  56.  
  57. for (auto& test : tests)
  58. {
  59. test.result = sine(test.a);
  60. }
  61.  
  62. int a = 0; // break here
  63. }
RAW Paste Data