1 /*************************************************************************** 2 * * 3 * LinuxSampler - modular, streaming capable sampler * 4 * * 5 * Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck * 6 * Copyright (C) 2005 - 2019 Christian Schoenebeck * 7 * * 8 * This program is free software; you can redistribute it and/or modify * 9 * it under the terms of the GNU General Public License as published by * 10 * the Free Software Foundation; either version 2 of the License, or * 11 * (at your option) any later version. * 12 * * 13 * This program is distributed in the hope that it will be useful, * 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 16 * GNU General Public License for more details. * 17 * * 18 * You should have received a copy of the GNU General Public License * 19 * along with this program; if not, write to the Free Software * 20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, * 21 * MA 02111-1307 USA * 22 ***************************************************************************/ 23 24 #ifndef __RT_MATH_H__ 25 #define __RT_MATH_H__ 26 27 #include <math.h> 28 #include <stdint.h> 29 #include "global_private.h" 30 31 /// Needed for calculating frequency ratio used to pitch a sample 32 #define TWELVEHUNDREDTH_ROOT_OF_TWO 1.000577789506555 33 34 enum implementation_t { 35 CPP, 36 ASM_X86_MMX_SSE 37 }; 38 39 /** @brief Real Time Math Base Class 40 * 41 * Math functions for real time operation. This base class contains all 42 * non-template methods. 43 */ 44 class RTMathBase { 45 public: 46 /** 47 * High resolution time stamp. 48 */ 49 typedef uint32_t time_stamp_t; 50 51 typedef uint64_t usecs_t; 52 53 /** 54 * We read the processor's cycle count register as a reference 55 * for the real time. These are of course only abstract values 56 * with arbitrary time entity, but that's not a problem as long 57 * as we calculate relatively. 58 * 59 * @see unsafeMicroSeconds() 60 */ 61 static time_stamp_t CreateTimeStamp(); 62 63 /** 64 * Calculates the frequency ratio for a pitch value given in cents 65 * (assuming equal tempered scale of course, divided into 12 66 * semitones per octave and 100 cents per semitone). 67 * 68 * Note: CONFIG_MAX_PITCH (defined in config.h) has to be defined to an 69 * appropriate value, otherwise the behavior of this function is 70 * undefined, but most probably if CONFIG_MAX_PITCH is too small, the 71 * application will crash due to segmentation fault here. 72 * 73 * @param cents - pitch value in cents (+1200 cents means +1 octave) 74 * @returns frequency ratio (e.g. +2.0 for +1 octave) 75 */ CentsToFreqRatio(double Cents)76 inline static double CentsToFreqRatio(double Cents) { 77 int index_int = (int) (Cents); // integer index 78 float index_fract = Cents - index_int; // fractional part of index 79 return pCentsToFreqTable[index_int] + index_fract * (pCentsToFreqTable[index_int+1] - pCentsToFreqTable[index_int]); 80 } 81 82 /** 83 * Slower version of CentsToFreqRatio, for big values. 84 * 85 * @param cents - pitch value in cents (+1200 cents means +1 octave) 86 * @returns frequency ratio (e.g. +2.0 for +1 octave) 87 */ CentsToFreqRatioUnlimited(double Cents)88 static double CentsToFreqRatioUnlimited(double Cents) { 89 int octaves = int(Cents / 1200); 90 double x = CentsToFreqRatio(Cents - octaves * 1200); 91 return octaves < 0 ? x / (1 << -octaves) : x * (1 << octaves); 92 } 93 94 /** 95 * Inverse function to CentsToFreqRatio(). This function is a bit 96 * slow, so it should not be called too frequently. 97 */ FreqRatioToCents(double FreqRatio)98 static double FreqRatioToCents(double FreqRatio) { 99 return log(FreqRatio) / log(TWELVEHUNDREDTH_ROOT_OF_TWO); 100 } 101 102 /** 103 * Calculates the linear ratio value representation (linear scale) 104 * of the @a decibel value provided (exponential scale). 105 * 106 * The context of audio acoustic sound pressure levels is assumed, and 107 * hence the field version of the dB unit is used here (which uses a 108 * linear factor of 20). This function is a bit slow, so it should 109 * not be called too frequently. 110 * 111 * @param decibel - sound pressure level in dB 112 * @returns linear ratio of the supplied dB value 113 * @see LinRatioToDecibel() as inverse function 114 */ DecibelToLinRatio(float decibel)115 static float DecibelToLinRatio(float decibel) { 116 return powf(10.f, decibel / 20.f); 117 } 118 119 /** 120 * Calculates the decibel value (exponential scale) of the @a linear 121 * ratio value representation (linear scale) provided. 122 * 123 * The context of audio acoustic sound pressure levels is assumed, and 124 * hence the field version of the dB unit is used here (which uses a 125 * linear factor of 20). This function is a bit slow, so it should 126 * not be called too frequently. 127 * 128 * @param linear - sound pressure level as linear ratio value (linear scale) 129 * @returns dB value representation 130 * @see DecibelToLinRatio() as inverse function 131 */ LinRatioToDecibel(float linear)132 static float LinRatioToDecibel(float linear) { 133 return 20.f * log10f(linear); 134 } 135 136 /** 137 * Calculates the relatively summed average of a set of values. 138 * 139 * @param current - the current avaerage value of all previously summed values 140 * @param sample - new value to be applied as summed average to the existing values 141 * @param n - amount of sample values applied so far 142 * @returns new average value of all summed values (including the new @a sample) 143 */ 144 template<typename T_int> RelativeSummedAvg(float current,float sample,T_int n)145 inline static float RelativeSummedAvg(float current, float sample, T_int n) { 146 return current + (sample - current) / float(n); 147 } 148 149 /** 150 * Clock source to use for getting the current time. 151 */ 152 enum clock_source_t { 153 real_clock, ///< Use this to measure time that passed in reality (no matter if process got suspended). 154 process_clock, ///< Use this to measure only the CPU execution time of the current process (if the process got suspended, the clock is paused as well). 155 thread_clock, ///< Use this to measure only the CPU execution time of the current thread (if the process got suspended or another thread is executed, the clock is paused as well). 156 }; 157 158 /** 159 * Returns a time stamp of the current time in microseconds (in 160 * probably real-time @b unsafe way). There is no guarantee about 161 * what the returned amount of microseconds relates to (i.e. 162 * microseconds since epoch, microseconds since system uptime, ...). 163 * So you should only use it to calculate time differences between 164 * values taken with this method. 165 * 166 * @b CAUTION: This method may not @b NOT be real-time safe! On some 167 * systems it could be RT safe, but there is no guarantee whatsoever! 168 * So this method should only be used for debugging, benchmarking and 169 * other developing purposes ! 170 * 171 * For creating time stamps in real-time context, use 172 * CreateTimeStamp() instead. 173 * 174 * @param source - the actual clock to use for getting the current 175 * time, note that the various clock sources may not 176 * be implemented on all systems 177 * @returns time stamp in microseconds 178 * 179 * @see CreateTimeStamp() 180 */ 181 static usecs_t unsafeMicroSeconds(clock_source_t source); 182 183 /** 184 * 'Equalness' comparison between two 32 bit floating point numbers 185 * (that is of single precision data type). 186 * 187 * This methods deals with the expected tolerance of the floating point 188 * representation. 189 * 190 * @return true if the two numbers can be expected to be equal 191 */ 192 static bool fEqual32(float a, float b); 193 194 /** 195 * 'Equalness' comparison between two 64 bit floating point numbers 196 * (that is of double precision data type). 197 * 198 * This methods deals with the expected tolerance of the floating point 199 * representation. 200 * 201 * @return true if the two numbers can be expected to be equal 202 */ 203 static bool fEqual64(double a, double b); 204 205 private: 206 static float* pCentsToFreqTable; 207 208 static float* InitCentsToFreqTable(); 209 }; 210 211 /** @brief Real Time Math 212 * 213 * This is a template which provides customized methods for the desired low 214 * level implementation. The ASM_X86_MMX_SSE implementation of each method 215 * for example doesn't use 387 FPU instruction. This is needed for MMX 216 * algorithms which do not allow mixed MMX and 387 instructions. 217 */ 218 template<implementation_t IMPL = CPP> 219 class __RTMath : public RTMathBase { 220 public: 221 // conversion using truncate Int(const float a)222 inline static int Int(const float a) { 223 switch (IMPL) { 224 #if CONFIG_ASM && ARCH_X86 225 case ASM_X86_MMX_SSE: { 226 int ret; 227 asm ( 228 "cvttss2si %1, %0 # convert to int\n\t" 229 : "=r" (ret) 230 : "m" (a) 231 ); 232 return ret; 233 } 234 #endif // CONFIG_ASM && ARCH_X86 235 default: { 236 return (int) a; 237 } 238 } 239 } 240 241 //for doubles and everything else except floats Int(const T_a a)242 template<class T_a> inline static int Int(const T_a a) { 243 return (int) a; 244 } 245 Float(const int a)246 inline static float Float(const int a) { 247 switch (IMPL) { 248 #if CONFIG_ASM && ARCH_X86 249 case ASM_X86_MMX_SSE: { 250 float ret; 251 asm ( 252 "cvtsi2ss %1, %%xmm0 # convert to float\n\t" 253 "movss %%xmm0,%0 # output\n\t" 254 : "=m" (ret) 255 : "r" (a) 256 ); 257 return ret; 258 } 259 #endif // CONFIG_ASM && ARCH_X86 260 default: { 261 return (float) a; 262 } 263 } 264 } 265 266 #if 0 267 //for everything except ints 268 template<class T_a> inline static float Float(T_a a) { 269 return (float) a; 270 } 271 #endif 272 Sum(const float & a,const float & b)273 inline static float Sum(const float& a, const float& b) { 274 switch (IMPL) { 275 #if CONFIG_ASM && ARCH_X86 276 case ASM_X86_MMX_SSE: { 277 float ret; 278 asm ( 279 "movss %1, %%xmm0 # load a\n\t" 280 "addss %2, %%xmm0 # a + b\n\t" 281 "movss %%xmm0, %0 # output\n\t" 282 : "=m" (ret) 283 : "m" (a), "m" (b) 284 ); 285 return ret; 286 } 287 #endif // CONFIG_ASM && ARCH_X86 288 default: { 289 return (a + b); 290 } 291 } 292 } 293 Sum(const T_a a,const T_b b)294 template<class T_a, class T_b> inline static T_a Sum(const T_a a, const T_b b) { 295 return (a + b); 296 } 297 Sub(const float & a,const float & b)298 inline static float Sub(const float& a, const float& b) { 299 switch (IMPL) { 300 #if CONFIG_ASM && ARCH_X86 301 case ASM_X86_MMX_SSE: { 302 float ret; 303 asm ( 304 "movss %1, %%xmm0 # load a\n\t" 305 "subss %2, %%xmm0 # a - b\n\t" 306 "movss %%xmm0, %0 # output\n\t" 307 : "=m" (ret) 308 : "m" (a), "m" (b) 309 ); 310 return ret; 311 } 312 #endif // CONFIG_ASM && ARCH_X86 313 default: { 314 return (a - b); 315 } 316 } 317 } 318 Sub(const T_a a,const T_b b)319 template<class T_a, class T_b> inline static T_a Sub(const T_a a, const T_b b) { 320 return (a - b); 321 } 322 Mul(const float a,const float b)323 inline static float Mul(const float a, const float b) { 324 switch (IMPL) { 325 #if CONFIG_ASM && ARCH_X86 326 case ASM_X86_MMX_SSE: { 327 float ret; 328 asm ( 329 "movss %1, %%xmm0 # load a\n\t" 330 "mulss %2, %%xmm0 # a * b\n\t" 331 "movss %%xmm0, %0 # output\n\t" 332 : "=m" (ret) 333 : "m" (a), "m" (b) 334 ); 335 return ret; 336 } 337 #endif // CONFIG_ASM && ARCH_X86 338 default: { 339 return (a * b); 340 } 341 } 342 } 343 Mul(const T_a a,const T_b b)344 template<class T_a, class T_b> inline static T_a Mul(const T_a a, const T_b b) { 345 return (a * b); 346 } 347 Div(const float a,const float b)348 inline static float Div(const float a, const float b) { 349 switch (IMPL) { 350 #if CONFIG_ASM && ARCH_X86 351 case ASM_X86_MMX_SSE: { 352 float ret; 353 asm ( 354 "movss %1, %%xmm0 # load a\n\t" 355 "divss %2, %%xmm0 # a / b\n\t" 356 "movss %%xmm0, %0 # output\n\t" 357 : "=m" (ret) 358 : "m" (a), "m" (b) 359 ); 360 return ret; 361 } 362 #endif // CONFIG_ASM && ARCH_X86 363 default: { 364 return (a / b); 365 } 366 } 367 } 368 Div(const T_a a,const T_b b)369 template<class T_a, class T_b> inline static T_a Div(const T_a a, const T_b b) { 370 return (a / b); 371 } 372 Min(const float a,const float b)373 inline static float Min(const float a, const float b) { 374 switch (IMPL) { 375 #if CONFIG_ASM && ARCH_X86 376 case ASM_X86_MMX_SSE: { 377 float ret; 378 asm ( 379 "movss %1, %%xmm0 # load a\n\t" 380 "minss %2, %%xmm0 # Minimum(a, b)\n\t" 381 "movss %%xmm0, %0 # output\n\t" 382 : "=m" (ret) 383 : "m" (a), "m" (b) 384 ); 385 return ret; 386 } 387 #endif // CONFIG_ASM && ARCH_X86 388 default: { 389 return std::min(a, b); 390 } 391 } 392 } 393 Min(const T_a a,const T_b b)394 template<class T_a, class T_b> inline static T_a Min(const T_a a, const T_b b) { 395 return (b < a) ? b : a; 396 } 397 Max(const float a,const float b)398 inline static float Max(const float a, const float b) { 399 switch (IMPL) { 400 #if CONFIG_ASM && ARCH_X86 401 case ASM_X86_MMX_SSE: { 402 float ret; 403 asm ( 404 "movss %1, %%xmm0 # load a\n\t" 405 "maxss %2, %%xmm0 # Maximum(a, b)\n\t" 406 "movss %%xmm0, %0 # output\n\t" 407 : "=m" (ret) 408 : "m" (a), "m" (b) 409 ); 410 return ret; 411 } 412 #endif // CONFIG_ASM && ARCH_X86 413 default: { 414 return std::max(a, b); 415 } 416 } 417 } 418 Max(const T_a a,const T_b b)419 template<class T_a, class T_b> inline static T_a Max(const T_a a, const T_b b) { 420 return (b > a) ? b : a; 421 } 422 Fmodf(const float & a,const float & b)423 inline static float Fmodf(const float &a, const float &b) { 424 switch (IMPL) { 425 #if CONFIG_ASM && ARCH_X86 426 case ASM_X86_MMX_SSE: { 427 float ret; 428 asm ( 429 "movss %1, %%xmm0 # load a\n\t" 430 "movss %2, %%xmm1 # load b\n\t" 431 "movss %%xmm0,%%xmm2\n\t" 432 "divss %%xmm1, %%xmm2 # xmm2 = a / b\n\t" 433 "cvttss2si %%xmm2, %%ecx #convert to int\n\t" 434 "cvtsi2ss %%ecx, %%xmm2 #convert back to float\n\t" 435 "mulss %%xmm1, %%xmm2 # xmm2 = b * int(a/b)\n\t" 436 "subss %%xmm2, %%xmm0 #sub a\n\t" 437 "movss %%xmm0, %0 # output\n\t" 438 : "=m" (ret) 439 : "m" (a), "m" (b) 440 : "%ecx" 441 ); 442 return ret; 443 } 444 #endif // CONFIG_ASM && ARCH_X86 445 default: { 446 return fmodf(a, b); 447 } 448 } 449 } 450 }; 451 452 /// convenience typedef for using the default implementation (which is CPP) 453 typedef __RTMath<> RTMath; 454 455 #endif // __RT_MATH_H__ 456