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 }