You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

161 lines
4.4KB

  1. // Copyright 2012 Olivier Gillet.
  2. //
  3. // Author: Olivier Gillet (ol.gillet@gmail.com)
  4. //
  5. // Permission is hereby granted, free of charge, to any person obtaining a copy
  6. // of this software and associated documentation files (the "Software"), to deal
  7. // in the Software without restriction, including without limitation the rights
  8. // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  9. // copies of the Software, and to permit persons to whom the Software is
  10. // furnished to do so, subject to the following conditions:
  11. //
  12. // The above copyright notice and this permission notice shall be included in
  13. // all copies or substantial portions of the Software.
  14. //
  15. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  16. // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  18. // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  19. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  20. // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  21. // THE SOFTWARE.
  22. //
  23. // See http://creativecommons.org/licenses/MIT/ for more information.
  24. //
  25. // -----------------------------------------------------------------------------
  26. //
  27. // DSP utility routines.
  28. #ifndef STMLIB_UTILS_DSP_DSP_H_
  29. #define STMLIB_UTILS_DSP_DSP_H_
  30. #include "stmlib/stmlib.h"
  31. #include <cmath>
  32. #include <math.h>
  33. namespace stmlib {
  34. #define MAKE_INTEGRAL_FRACTIONAL(x) \
  35. int32_t x ## _integral = static_cast<int32_t>(x); \
  36. float x ## _fractional = x - static_cast<float>(x ## _integral);
  37. inline float Interpolate(const float* table, float index, float size) {
  38. index *= size;
  39. MAKE_INTEGRAL_FRACTIONAL(index)
  40. float a = table[index_integral];
  41. float b = table[index_integral + 1];
  42. return a + (b - a) * index_fractional;
  43. }
  44. inline float InterpolateHermite(const float* table, float index, float size) {
  45. index *= size;
  46. MAKE_INTEGRAL_FRACTIONAL(index)
  47. const float xm1 = table[index_integral - 1];
  48. const float x0 = table[index_integral + 0];
  49. const float x1 = table[index_integral + 1];
  50. const float x2 = table[index_integral + 2];
  51. const float c = (x1 - xm1) * 0.5f;
  52. const float v = x0 - x1;
  53. const float w = c + v;
  54. const float a = w + v + (x2 - x0) * 0.5f;
  55. const float b_neg = w + a;
  56. const float f = index_fractional;
  57. return (((a * f) - b_neg) * f + c) * f + x0;
  58. }
  59. inline float InterpolateWrap(const float* table, float index, float size) {
  60. index -= static_cast<float>(static_cast<int32_t>(index));
  61. index *= size;
  62. MAKE_INTEGRAL_FRACTIONAL(index)
  63. float a = table[index_integral];
  64. float b = table[index_integral + 1];
  65. return a + (b - a) * index_fractional;
  66. }
  67. #define ONE_POLE(out, in, coefficient) out += (coefficient) * ((in) - out);
  68. #define SLOPE(out, in, positive, negative) { \
  69. float error = (in) - out; \
  70. out += (error > 0 ? positive : negative) * error; \
  71. }
  72. #define SLEW(out, in, delta) { \
  73. float error = (in) - out; \
  74. float d = (delta); \
  75. if (error > d) { \
  76. error = d; \
  77. } else if (error < -d) { \
  78. error = -d; \
  79. } \
  80. out += error; \
  81. }
  82. inline float Crossfade(float a, float b, float fade) {
  83. return a + (b - a) * fade;
  84. }
  85. inline float SoftLimit(float x) {
  86. return x * (27.0f + x * x) / (27.0f + 9.0f * x * x);
  87. }
  88. inline float SoftClip(float x) {
  89. if (x < -3.0f) {
  90. return -1.0f;
  91. } else if (x > 3.0f) {
  92. return 1.0f;
  93. } else {
  94. return SoftLimit(x);
  95. }
  96. }
  97. #ifdef TEST
  98. inline int32_t Clip16(int32_t x) {
  99. if (x < -32768) {
  100. return -32768;
  101. } else if (x > 32767) {
  102. return 32767;
  103. } else {
  104. return x;
  105. }
  106. }
  107. inline uint16_t ClipU16(int32_t x) {
  108. if (x < 0) {
  109. return 0;
  110. } else if (x > 65535) {
  111. return 65535;
  112. } else {
  113. return x;
  114. }
  115. }
  116. #else
  117. inline int32_t Clip16(int32_t x) {
  118. int32_t result;
  119. __asm ("ssat %0, %1, %2" : "=r" (result) : "I" (16), "r" (x) );
  120. return result;
  121. }
  122. inline uint32_t ClipU16(int32_t x) {
  123. uint32_t result;
  124. __asm ("usat %0, %1, %2" : "=r" (result) : "I" (16), "r" (x) );
  125. return result;
  126. }
  127. #endif
  128. #ifdef TEST
  129. inline float Sqrt(float x) {
  130. return sqrtf(x);
  131. }
  132. #else
  133. inline float Sqrt(float x) {
  134. float result;
  135. __asm ("vsqrt.f32 %0, %1" : "=w" (result) : "w" (x) );
  136. return result;
  137. }
  138. #endif
  139. inline int16_t SoftConvert(float x) {
  140. return Clip16(static_cast<int32_t>(SoftLimit(x * 0.5f) * 32768.0f));
  141. }
  142. } // namespace stmlib
  143. #endif // STMLIB_UTILS_DSP_DSP_H_