1 /**
2 * Copyright: Copyright Auburn Sounds 2015-2016
3 * License:   $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
4 * Authors:   Guillaume Piolat
5 */
6 module dplug.core.math;
7 
8 import std.math;
9 
10 version(LDC)
11 {
12     import ldc.intrinsics;
13 }
14 
15 immutable real TAU = PI * 2;
16 
17 /// Map linearly x from the range [a, b] to the range [c, d]
18 T linmap(T)(T value, T a, T b, T c, T d) pure nothrow @nogc
19 {
20     return c + (d - c) * (value - a) / (b - a);
21 }
22 
23 /// map [0..1] to [min..max] logarithmically
24 /// min and max must be all > 0, t in [0..1]
25 T logmap(T)(T t, T min, T max) pure nothrow @nogc
26 {
27     assert(min < max);
28     return min * exp(t * log(max / min));
29 }
30 
31 /// Hermite interpolation.
32 T hermite(T)(T frac_pos, T xm1, T x0, T x1, T x2) pure nothrow @nogc
33 {
34     T c = (x1 - xm1) * 0.5f;
35     T v = x0 - x1;
36     T w = c + v;
37     T a = w + v + (x2 - x0) * 0.5f;
38     T b_neg = w + a;
39     return ((((a * frac_pos) - b_neg) * frac_pos + c) * frac_pos + x0);
40 }
41 
42 /// Convert from dB to float.
43 T deciBelToFloat(T)(T dB) pure nothrow @nogc
44 {
45     static immutable T ln10_20 = cast(T)LN10 / 20;
46     return exp(dB * ln10_20);
47 }
48 
49 /// Convert from float to dB
50 T floatToDeciBel(T)(T x) pure nothrow @nogc
51 {
52     static immutable T f20_ln10 = 20 / cast(T)LN10;
53     return log(x) * f20_ln10;
54 }
55 
56 /// Is this integer odd?
57 bool isOdd(T)(T i) pure nothrow @nogc
58 {
59     return (i & 1) != 0;
60 }
61 
62 /// Is this integer even?
63 bool isEven(T)(T i) pure nothrow @nogc
64 {
65     return (i & 1) == 0;
66 }
67 
68 /// Returns: true of i is a power of 2.
69 bool isPowerOfTwo(int i) pure nothrow @nogc
70 {
71     assert(i >= 0);
72     return (i != 0) && ((i & (i - 1)) == 0);
73 }
74 
75 /// Computes next power of 2.
76 int nextPowerOf2(int i) pure nothrow @nogc
77 {
78     int v = i - 1;
79     v |= v >> 1;
80     v |= v >> 2;
81     v |= v >> 4;
82     v |= v >> 8;
83     v |= v >> 16;
84     v++;
85     return v;
86 }
87 
88 /// Computes next power of 2.
89 long nextPowerOf2(long i) pure nothrow @nogc
90 {
91     long v = i - 1;
92     v |= v >> 1;
93     v |= v >> 2;
94     v |= v >> 4;
95     v |= v >> 8;
96     v |= v >> 16;
97     v |= v >> 32;
98     v++;
99     return v;
100 }
101 
102 /// Returns: x so that (1 << x) >= i
103 int iFloorLog2(int i) pure nothrow @nogc
104 {
105     assert(i >= 1);
106     int result = 0;
107     while (i > 1)
108     {
109         i = i / 2;
110         result = result + 1;
111     }
112     return result;
113 }
114 
115 /// Mapping from MIDI notes to frequency
116 deprecated("use convertMIDINoteToFrequency instead") alias MIDIToFrequency = convertMIDINoteToFrequency;
117 T convertMIDINoteToFrequency(T)(T note) pure nothrow @nogc
118 {
119     return 440.0f * pow(2.0, (note - 69.0f) / 12.0f);
120 }
121 
122 /// Mapping from frequency to MIDI notes
123 deprecated("use convertFrequencyToMIDINote instead") alias frequencyToMIDI = convertFrequencyToMIDINote;
124 T convertFrequencyToMIDINote(T)(T frequency) pure nothrow @nogc
125 {
126     return 69.0f + 12.0f * log2(frequency / 440.0f);
127 }
128 
129 /// Fletcher and Munson equal-loudness curve
130 /// Reference: Xavier Serra thesis (1989).
131 T equalLoudnessCurve(T)(T frequency) pure nothrow @nogc
132 {
133     T x = cast(T)0.05 + 4000 / frequency;
134     return x * ( cast(T)10 ^^ x);
135 }
136 
137 /// Cardinal sine
138 T sinc(T)(T x) pure nothrow @nogc
139 {
140     if (cast(T)(1) + x * x == cast(T)(1))
141         return 1;
142     else
143         return sin(cast(T)PI * x) / (cast(T)PI * x);
144 }
145 
146 /// Params:
147 ///    timeConstantInSeconds time after which the amplitude is only 37% of the original.
148 /// Returns: Multiplier for this time constant and sampling rate.
149 double expDecayFactor(double timeConstantInSeconds, double samplerate) pure nothrow @nogc
150 {
151     // 1 - exp(-time * sampleRate) would yield innacuracies
152     return -expm1(-1.0 / (timeConstantInSeconds * samplerate));
153 }
154 
155 /// Give back a phase between -PI and PI
156 T normalizePhase(T)(T phase) nothrow @nogc
157 {
158     enum bool Assembly = D_InlineAsm_Any && !(is(Unqual!T == real));
159 
160     static if (Assembly)
161     {
162         T k_TAU = PI * 2;
163         T result = phase;
164         asm nothrow @nogc
165         {
166             fld k_TAU;    // TAU
167             fld result;    // phase | TAU
168             fprem1;       // normalized(phase) | TAU
169             fstp result;   // TAU
170             fstp ST(0);   //
171         }
172         return result;
173     }
174     else
175     {
176         T res = fmod(phase, cast(T)TAU);
177         if (res > PI)
178             res -= TAU;
179         if (res < -PI)
180             res += TAU;
181         return res;
182     }
183 }
184 
185 unittest
186 {
187     assert(approxEqual(normalizePhase!real(TAU), 0));
188 
189     assert(approxEqual(normalizePhase!float(0.1f), 0.1f));
190     assert(approxEqual(normalizePhase!float(TAU + 0.1f), 0.1f));
191 
192     assert(approxEqual(normalizePhase!double(-0.1f), -0.1f));
193     assert(approxEqual(normalizePhase!double(-TAU - 0.1f), -0.1f));
194 
195     bool approxEqual(T)(T a, T b) nothrow @nogc
196     {
197         return (a - b) < 1e-7;
198     }
199 }
200 
201 /// Quick and dirty sawtooth for testing purposes.
202 T rawSawtooth(T)(T x) pure nothrow @nogc
203 {
204     return normalizePhase(x) / (cast(T)PI);
205 }
206 
207 /// Quick and dirty triangle for testing purposes.
208 T rawTriangle(T)(T x) pure nothrow @nogc
209 {
210     return 1 - normalizePhase(x) / cast(T)PI_2;
211 }
212 
213 /// Quick and dirty square for testing purposes.
214 T rawSquare(T)(T x) pure nothrow @nogc
215 {
216     return normalizePhase(x) > 0 ? 1 : -1;
217 }
218 
219 T computeRMS(T)(T[] samples) pure nothrow @nogc
220 {
221     double sum = 0;
222     foreach(sample; samples)
223         sum += sample * sample;
224     return sqrt(sum / cast(int)samples.length);
225 }
226 
227 unittest
228 {
229     double[] d = [4, 5, 6];
230     computeRMS(d);
231 }
232 
233 
234 version(D_InlineAsm_X86)
235     private enum D_InlineAsm_Any = true;
236 else version(D_InlineAsm_X86_64)
237     private enum D_InlineAsm_Any = true;
238 else
239     private enum D_InlineAsm_Any = false;
240 
241 // These functions trade correctness for speed
242 // The contract is that they don't check for infinity or NaN
243 // and assume small finite numbers instead.
244 // Don't rely on them being correct for your situation: test them.
245 
246 ///
247 T fast_pow(T)(T val, T power)
248 {
249     version(LDC)
250         return llvm_pow(val, power);
251     else
252         return pow(val, power);
253 }
254 
255 ///
256 T fast_exp(T)(T val, T)
257 {
258     version(LDC)
259         return llvm_exp(val);
260     else
261         return exp(val);
262 }
263 
264 ///
265 T fast_log(T)(T val)
266 {
267     version(LDC)
268         return llvm_log(val);
269     else
270         return log(val);
271 }
272 
273 ///
274 T fast_floor(T)(T val)
275 {
276     version(LDC)
277         return llvm_floor(val);
278     else
279         return log(val);
280 }
281 
282 ///
283 T fast_ceil(T)(T val)
284 {
285     version(LDC)
286         return llvm_ceil(val);
287     else
288         return ceil(val);
289 }
290 
291 ///
292 T fast_trunc(T)(T val)
293 {
294     version(LDC)
295         return llvm_trunc(val);
296     else
297         return trunc(val);
298 }
299 
300 ///
301 T fast_round(T)(T val)
302 {
303     version(LDC)
304         return llvm_round(val);
305     else
306         return round(val);
307 }
308 
309 ///
310 T fast_exp2(T)(T val)
311 {
312     version(LDC)
313         return llvm_exp2(val);
314     else
315         return exp2(val);
316 }
317 
318 ///
319 T fast_log10(T)(T val)
320 {
321     version(LDC)
322         return llvm_log10(val);
323     else
324         return log10(val);
325 }
326 
327 ///
328 T fast_log2(T)(T val)
329 {
330     version(LDC)
331         return llvm_log2(val);
332     else
333         return log2(val);
334 }
335 
336 /// Linear interpolation, akin to GLSL's mix.
337 S lerp(S, T)(S a, S b, T t) pure nothrow @nogc 
338     if (is(typeof(t * b + (1 - t) * a) : S))
339 {
340     return t * b + (1 - t) * a;
341 }
342 
343 /// Same as GLSL smoothstep function.
344 /// See: http://en.wikipedia.org/wiki/Smoothstep
345 T smoothStep(T)(T a, T b, T t) pure nothrow @nogc 
346 {
347     if (t <= a)
348         return 0;
349     else if (t >= b)
350         return 1;
351     else
352     {
353         T x = (t - a) / (b - a);
354         return x * x * (3 - 2 * x);
355     }
356 }
357 
358 /// SSE approximation of reciprocal square root.
359 T inverseSqrt(T)(T x) pure nothrow @nogc if (is(T : float) || is(T: double))
360 {
361     version(AsmX86)
362     {
363         static if (is(T == float))
364         {
365             float result;
366 
367             asm pure nothrow @nogc 
368             {
369                 movss XMM0, x; 
370                 rsqrtss XMM0, XMM0; 
371                 movss result, XMM0; 
372             }
373             return result;
374         }
375         else
376             return 1 / sqrt(x);
377     }
378     else
379         return 1 / sqrt(x);
380 }
381 
382 unittest
383 {
384     assert(abs( inverseSqrt!float(1) - 1) < 1e-3 );
385     assert(abs( inverseSqrt!double(1) - 1) < 1e-3 );
386 }
387 
388 /// Computes a normalized frequency form a frequency.
389 float convertFrequencyToNormalizedFrequency(float frequencyHz, float samplingRate) pure nothrow @nogc
390 {
391     return frequencyHz / samplingRate;
392 }
393 
394 /// Computes a frequency.
395 float convertNormalizedFreqyencyToFrequency(float frequencyCyclesPerSample, float samplingRate) pure nothrow @nogc
396 {
397     return frequencyCyclesPerSample * samplingRate;
398 }