1 // Copyright 2012 Emilie Gillet.
2 //
3 // Author: Emilie Gillet (emilie.o.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
29 #ifndef STMLIB_UTILS_DSP_DSP_H_
30 #define STMLIB_UTILS_DSP_DSP_H_
31
32 #include "stmlib/stmlib.h"
33
34 #include <cmath>
35 #include <math.h>
36
37 namespace stmlib {
38
39 #define MAKE_INTEGRAL_FRACTIONAL(x) \
40 int32_t x ## _integral = static_cast<int32_t>(x); \
41 float x ## _fractional = x - static_cast<float>(x ## _integral);
42
43 inline float Interpolate(const float* table, float index, float size) {
44 index *= size;
45 MAKE_INTEGRAL_FRACTIONAL(index)
CosineOscillator()46 float a = table[index_integral];
~CosineOscillator()47 float b = table[index_integral + 1];
48 return a + (b - a) * index_fractional;
49 }
Init(float frequency)50
51
52 inline float InterpolateHermite(const float* table, float index, float size) {
53 index *= size;
54 MAKE_INTEGRAL_FRACTIONAL(index)
55 const float xm1 = table[index_integral - 1];
56 const float x0 = table[index_integral + 0];
57 const float x1 = table[index_integral + 1];
58 const float x2 = table[index_integral + 2];
59 const float c = (x1 - xm1) * 0.5f;
60 const float v = x0 - x1;
61 const float w = c + v;
62 const float a = w + v + (x2 - x0) * 0.5f;
63 const float b_neg = w + a;
64 const float f = index_fractional;
65 return (((a * f) - b_neg) * f + c) * f + x0;
66 }
67
68 inline float InterpolateWrap(const float* table, float index, float size) {
69 index -= static_cast<float>(static_cast<int32_t>(index));
70 index *= size;
71 MAKE_INTEGRAL_FRACTIONAL(index)
72 float a = table[index_integral];
73 float b = table[index_integral + 1];
74 return a + (b - a) * index_fractional;
75 }
Start()76
77 #define ONE_POLE(out, in, coefficient) out += (coefficient) * ((in) - out);
78 #define SLOPE(out, in, positive, negative) { \
79 float error = (in) - out; \
80 out += (error > 0 ? positive : negative) * error; \
value()81 }
82 #define SLEW(out, in, delta) { \
83 float error = (in) - out; \
84 float d = (delta); \
85 if (error > d) { \
86 error = d; \
87 } else if (error < -d) { \
88 error = -d; \
89 } \
90 out += error; \
91 }
92
93 inline float Crossfade(float a, float b, float fade) {
94 return a + (b - a) * fade;
95 }
96
97 inline float SoftLimit(float x) {
98 return x * (27.0f + x * x) / (27.0f + 9.0f * x * x);
99 }
100
101 inline float SoftClip(float x) {
102 if (x < -3.0f) {
103 return -1.0f;
104 } else if (x > 3.0f) {
105 return 1.0f;
106 } else {
107 return SoftLimit(x);
108 }
109 }
110
111 #ifdef TEST
112 inline int32_t Clip16(int32_t x) {
113 if (x < -32768) {
114 return -32768;
115 } else if (x > 32767) {
116 return 32767;
117 } else {
118 return x;
119 }
120 }
121 inline uint16_t ClipU16(int32_t x) {
122 if (x < 0) {
123 return 0;
124 } else if (x > 65535) {
125 return 65535;
126 } else {
127 return x;
128 }
129 }
130 #else
131 inline int32_t Clip16(int32_t x) {
132 int32_t result;
133 __asm ("ssat %0, %1, %2" : "=r" (result) : "I" (16), "r" (x) );
134 return result;
135 }
136 inline uint32_t ClipU16(int32_t x) {
137 uint32_t result;
138 __asm ("usat %0, %1, %2" : "=r" (result) : "I" (16), "r" (x) );
139 return result;
140 }
141 #endif
142
143 #ifdef TEST
144 inline float Sqrt(float x) {
145 return sqrtf(x);
146 }
147 #else
148 inline float Sqrt(float x) {
149 float result;
150 __asm ("vsqrt.f32 %0, %1" : "=w" (result) : "w" (x) );
151 return result;
152 }
153 #endif
154
155 inline int16_t SoftConvert(float x) {
156 return Clip16(static_cast<int32_t>(SoftLimit(x * 0.5f) * 32768.0f));
157 }
158
159 } // namespace stmlib
160
161 #endif // STMLIB_UTILS_DSP_DSP_H_