src/calf/primitives.h

   1 /* Calf DSP Library
   2  * DSP primitives.
   3  *
   4  * Copyright (C) 2001-2007 Krzysztof Foltman
   5  *
   6  * This program is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2 of the License, or (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General
  17  * Public License along with this program; if not, write to the
  18  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  19  * Boston, MA 02111-1307, USA.
  20  */
  21 #ifndef __CALF_PRIMITIVES_H
  22 #define __CALF_PRIMITIVES_H
  23
  24 #include <assert.h>
  25 #include <stdint.h>
  26 #include <stdio.h>
  27 #include <string.h>
  28 #include <cmath>
  29 #include <cstdlib>
  30 #include <map>
  31
  32 namespace dsp {
  33
  34 /// Set a float to zero
  35 inline void zero(float &v) {
  36     v = 0.f;
  37 };
  38
  39 /// Set a double to zero
  40 inline void zero(double &v) {
  41     v = 0.0;
  42 };
  43
  44 /// Set 64-bit unsigned integer value to zero
  45 inline void zero(uint64_t &v) { v = 0; };
  46 /// Set 32-bit unsigned integer value to zero
  47 inline void zero(uint32_t &v) { v = 0; };
  48 /// Set 16-bit unsigned integer value to zero
  49 inline void zero(uint16_t &v) { v = 0; };
  50 /// Set 8-bit unsigned integer value to zero
  51 inline void zero(uint8_t &v) { v = 0; };
  52 /// Set 64-bit signed integer value to zero
  53 inline void zero(int64_t &v) { v = 0; };
  54 /// Set 32-bit signed integer value to zero
  55 inline void zero(int32_t &v) { v = 0; };
  56 /// Set 16-bit signed integer value to zero
  57 inline void zero(int16_t &v) { v = 0; };
  58 /// Set 8-bit signed integer value to zero
  59 inline void zero(int8_t &v) { v = 0; };
  60
  61 /// Set array (buffer or anything similar) to vector of zeroes
  62 template<class T>
  63 void zero(T *data, unsigned int size) {
  64     T value;
  65     dsp::zero(value);
  66     for (unsigned int i=0; i<size; i++)
  67         *data++ = value;
  68 }
  69
  70 /// Set array (buffer or anything similar) to vector of values
  71 template<class T>
  72 void fill(T *data, T value, unsigned int size) {
  73     for (unsigned int i=0; i<size; i++)
  74         *data++ = value;
  75 }
  76
  77 template<class T = float>struct stereo_sample {
  78     T left;
  79     T right;
  80     /// default constructor - preserves T's semantics (ie. no implicit initialization to 0)
  81     inline stereo_sample() {
  82     }
  83     inline stereo_sample(T _left, T _right) {
  84         left = _left;
  85         right = _right;
  86     }
  87     inline stereo_sample(T _both) {
  88         left = right = _both;
  89     }
  90     template<typename U>
  91     inline stereo_sample(const stereo_sample<U> &value) {
  92         left = value.left;
  93         right = value.right;
  94     }
  95     inline stereo_sample& operator=(const T &value) {
  96         left = right = value;
  97         return *this;
  98     }
  99     template<typename U>
 100     inline stereo_sample& operator=(const stereo_sample<U> &value) {
 101         left = value.left;
 102         right = value.right;
 103         return *this;
 104     }
 105 /*
 106     inline operator T() const {
 107         return (left+right)/2;
 108     }
 109 */
 110     inline stereo_sample& operator*=(const T &multiplier) {
 111         left *= multiplier;
 112         right *= multiplier;
 113         return *this;
 114     }
 115     inline stereo_sample& operator+=(const stereo_sample<T> &value) {
 116         left += value.left;
 117         right += value.right;
 118         return *this;
 119     }
 120     inline stereo_sample& operator-=(const stereo_sample<T> &value) {
 121         left -= value.left;
 122         right -= value.right;
 123         return *this;
 124     }
 125     template<typename U> inline stereo_sample<U> operator*(const U &value) const {
 126         return stereo_sample<U>(left*value, right*value);
 127     }
 128     /*inline stereo_sample<float> operator*(float value) const {
 129         return stereo_sample<float>(left*value, right*value);
 130     }
 131     inline stereo_sample<double> operator*(double value) const {
 132         return stereo_sample<double>(left*value, right*value);
 133     }*/
 134     inline stereo_sample<T> operator+(const stereo_sample<T> &value) {
 135         return stereo_sample(left+value.left, right+value.right);
 136     }
 137     inline stereo_sample<T> operator-(const stereo_sample<T> &value) {
 138         return stereo_sample(left-value.left, right-value.right);
 139     }
 140     inline stereo_sample<T> operator+(const T &value) {
 141         return stereo_sample(left+value, right+value);
 142     }
 143     inline stereo_sample<T> operator-(const T &value) {
 144         return stereo_sample(left-value, right-value);
 145     }
 146     inline stereo_sample<float> operator+(float value) {
 147         return stereo_sample<float>(left+value, right+value);
 148     }
 149     inline stereo_sample<float> operator-(float value) {
 150         return stereo_sample<float>(left-value, right-value);
 151     }
 152     inline stereo_sample<double> operator+(double value) {
 153         return stereo_sample<double>(left+value, right+value);
 154     }
 155     inline stereo_sample<double> operator-(double value) {
 156         return stereo_sample<double>(left-value, right-value);
 157     }
 158 };
 159
 160 /// Multiply constant by stereo_value
 161 template<class T>
 162 inline stereo_sample<T> operator*(const T &value, const stereo_sample<T> &value2) {
 163     return stereo_sample<T>(value2.left*value, value2.right*value);
 164 }
 165
 166 /// Add constant to stereo_value
 167 template<class T>
 168 inline stereo_sample<T> operator+(const T &value, const stereo_sample<T> &value2) {
 169     return stereo_sample<T>(value2.left+value, value2.right+value);
 170 }
 171
 172 /// Subtract stereo_value from constant (yields stereo_value of course)
 173 template<class T>
 174 inline stereo_sample<T> operator-(const T &value, const stereo_sample<T> &value2) {
 175     return stereo_sample<T>(value-value2.left, value-value2.right);
 176 }
 177
 178 /// Shift value right by 'bits' bits (multiply by 2^-bits)
 179 template<typename T>
 180 inline stereo_sample<T> shr(stereo_sample<T> v, int bits = 1) {
 181     v.left = shr(v.left, bits);
 182     v.right = shr(v.right, bits);
 183     return v;
 184 }
 185
 186 /// Set a stereo_sample<T> value to zero
 187 template<typename T>
 188 inline void zero(stereo_sample<T> &v) {
 189     dsp::zero(v.left);
 190     dsp::zero(v.right);
 191 }
 192
 193 /// 'Small value' for integer and other types
 194 template<typename T>
 195 inline T small_value() {
 196     return 0;
 197 }
 198
 199 /// 'Small value' for floats (2^-24) - used for primitive underrun prevention. The value is pretty much arbitrary (allowing for 24-bit signals normalized to 1.0).
 200 template<>
 201 inline float small_value<float>() {
 202     return (1.0/16777216.0); // allows for 2^-24, should be enough for 24-bit DACs at least :)
 203 }
 204
 205 /// 'Small value' for doubles (2^-24) - used for primitive underrun prevention. The value is pretty much arbitrary.
 206 template<>
 207 inline double small_value<double>() {
 208     return (1.0/16777216.0);
 209 }
 210
 211 /// Convert a single value to single value = do nothing :) (but it's a generic with specialisation for stereo_sample)
 212 template<typename T>
 213 inline float mono(T v) {
 214     return v;
 215 }
 216
 217 /// Convert a stereo_sample to single value by averaging two channels
 218 template<typename T>
 219 inline T mono(stereo_sample<T> v) {
 220     return shr(v.left+v.right);
 221 }
 222
 223 /// Clip a value to [min, max]
 224 template<typename T>
 225 inline T clip(T value, T min, T max) {
 226     if (value < min) return min;
 227     if (value > max) return max;
 228     return value;
 229 }
 230
 231 /// Clip a double to [-1.0, +1.0]
 232 inline double clip11(double value) {
 233     double a = fabs(value);
 234     if (a<=1) return value;
 235     return (value<0) ? -1.0 : 1.0;
 236 }
 237
 238 /// Clip a float to [-1.0f, +1.0f]
 239 inline float clip11(float value) {
 240     float a = fabsf(value);
 241     if (a<=1) return value;
 242     return (value<0) ? -1.0f : 1.0f;
 243 }
 244
 245 /// Clip a double to [0.0, +1.0]
 246 inline double clip01(double value) {
 247     double a = fabs(value-0.5);
 248     if (a<=0.5) return value;
 249     return (a<0) ? -0.0 : 1.0;
 250 }
 251
 252 /// Clip a float to [0.0f, +1.0f]
 253 inline float clip01(float value) {
 254     float a = fabsf(value-0.5f);
 255     if (a<=0.5f) return value;
 256     return (value < 0) ? -0.0f : 1.0f;
 257 }
 258
 259 // Linear interpolation (mix-way between v1 and v2).
 260 template<typename T, typename U>
 261 inline T lerp(T v1, T v2, U mix) {
 262     return v1+(v2-v1)*mix;
 263 }
 264
 265 // Linear interpolation for stereo values (mix-way between v1 and v2).
 266 template<typename T>
 267 inline stereo_sample<T> lerp(stereo_sample<T> &v1, stereo_sample<T> &v2, float mix) {
 268     return stereo_sample<T>(v1.left+(v2.left-v1.left)*mix, v1.right+(v2.right-v1.right)*mix);
 269 }
 270
 271 /**
 272  * decay-only envelope (linear or exponential); deactivates itself when it goes below a set point (epsilon)
 273  */
 274 class decay
 275 {
 276     double value, initial;
 277     unsigned int age, mask;
 278     bool active;
 279 public:
 280     decay() {
 281         active = false;
 282         mask = 127;
 283         initial = value = 0.0;
 284     }
 285     inline bool get_active() {
 286         return active;
 287     }
 288     inline double get() {
 289         return active ? value : 0.0;
 290     }
 291     inline void set(double v) {
 292         initial = value = v;
 293         active = true;
 294         age = 0;
 295     }
 296     /// reinitialise envelope (must be called if shape changes from linear to exponential or vice versa in the middle of envelope)
 297     inline void reinit()
 298     {
 299         initial = value;
 300         age = 1;
 301     }
 302     inline void add(double v) {
 303         if (active)
 304             value += v;
 305         else
 306             value = v;
 307         initial = value;
 308         age = 0;
 309         active = true;
 310     }
 311     static inline double calc_exp_constant(double times, double cycles)
 312     {
 313         if (cycles < 1.0)
 314             cycles = 1.0;
 315         return pow(times, 1.0 / cycles);
 316     }
 317     inline void age_exp(double constant, double epsilon) {
 318         if (active) {
 319             if (!(age & mask))
 320                 value = initial * pow(constant, (double)age);
 321             else
 322                 value *= constant;
 323             if (value < epsilon)
 324                 active = false;
 325             age++;
 326         }
 327     }
 328     inline void age_lin(double constant, double epsilon) {
 329         if (active) {
 330             if (!(age & mask))
 331                 value = initial - constant * age;
 332             else
 333                 value -= constant;
 334             if (value < epsilon)
 335                 active = false;
 336             age++;
 337         }
 338     }
 339     inline void deactivate() {
 340         active = false;
 341         value = 0;
 342     }
 343 };
 344
 345 class scheduler;
 346
 347 class task {
 348 public:
 349     virtual void execute(scheduler *s)=0;
 350     virtual void dispose() { delete this; }
 351     virtual ~task() {}
 352 };
 353
 354 /// this scheduler is based on std::multimap, so it isn't very fast, I guess
 355 /// maybe some day it should be rewritten to use heapsort or something
 356 /// work in progress, don't use!
 357 class scheduler {
 358     std::multimap<unsigned int, task *> timeline;
 359     unsigned int time, next_task;
 360     bool eob;
 361     class end_buf_task: public task {
 362     public:
 363         scheduler *p;
 364         end_buf_task(scheduler *_p) : p(_p) {}
 365         virtual void execute(scheduler *s) { p->eob = true; }
 366         virtual void dispose() { }
 367     } eobt;
 368 public:
 369
 370     scheduler()
 371     : time(0)
 372     , next_task((unsigned)-1)
 373     , eob(true)
 374     , eobt (this)
 375     {
 376         time = 0;
 377         next_task = (unsigned)-1;
 378         eob = false;
 379     }
 380     inline bool is_next_tick() {
 381         if (time < next_task)
 382             return true;
 383         do_tasks();
 384     }
 385     inline void next_tick() {
 386         time++;
 387     }
 388     void set(int pos, task *t) {
 389         timeline.insert(std::pair<unsigned int, task *>(time+pos, t));
 390         next_task = timeline.begin()->first;
 391     }
 392     void do_tasks() {
 393         std::multimap<unsigned int, task *>::iterator i = timeline.begin();
 394         while(i != timeline.end() && i->first == time) {
 395             i->second->execute(this);
 396             i->second->dispose();
 397             timeline.erase(i);
 398         }
 399     }
 400     bool is_eob() {
 401         return eob;
 402     }
 403     void set_buffer_size(int count) {
 404         set(count, &eobt);
 405     }
 406 };
 407
 408 /**
 409  * Force "small enough" float value to zero
 410  */
 411 inline void sanitize(float &value)
 412 {
 413     // real number?
 414     if (std::abs(value) < small_value<float>())
 415         value = 0.f;
 416     // close to 0?
 417     const int val = *reinterpret_cast <const int *> (&value);
 418     if ((val & 0x7F800000) == 0 && (val & 0x007FFFFF) != 0)
 419         value = 0.f;
 420 }
 421 inline float _sanitize(float value)
 422 {
 423     if (std::abs(value) < small_value<float>())
 424         return 0.f;
 425     return value;
 426 }
 427
 428 /**
 429  * Force already-denormal float value to zero
 430  */
 431 inline void sanitize_denormal(float& value)
 432 {
 433     if (!std::isnormal(value))
 434          value = 0.f;
 435 }
 436
 437 /**
 438  * Force already-denormal float value to zero
 439  */
 440 inline void sanitize_denormal(double & value)
 441 {
 442     if (!std::isnormal(value))
 443          value = 0.f;
 444 }
 445
 446 /**
 447  * Force "small enough" double value to zero
 448  */
 449 inline void sanitize(double &value)
 450 {
 451     if (std::abs(value) < small_value<double>())
 452         value = 0.0;
 453 }
 454
 455 inline double _sanitize(double value)
 456 {
 457     if (std::abs(value) < small_value<double>())
 458         return 0.0;
 459
 460     return value;
 461 }
 462 /**
 463  * Force "small enough" stereo value to zero
 464  */
 465 template<class T>
 466 inline void sanitize(stereo_sample<T> &value)
 467 {
 468     sanitize(value.left);
 469     sanitize(value.right);
 470 }
 471
 472 inline float fract16(unsigned int value)
 473 {
 474     return (value & 0xFFFF) * (1.0 / 65536.0);
 475 }
 476
 477 /**
 478  * typical precalculated sine table
 479  */
 480 template<class T, int N, int Multiplier>
 481 class sine_table
 482 {
 483 public:
 484     static bool initialized;
 485     static T data[N+1];
 486     sine_table() {
 487         if (initialized)
 488             return;
 489         initialized = true;
 490         for (int i=0; i<N+1; i++)
 491             data[i] = (T)(Multiplier*sin(i*2*M_PI*(1.0/N)));
 492     }
 493 };
 494
 495 template<class T, int N, int Multiplier>
 496 bool sine_table<T,N,Multiplier>::initialized = false;
 497
 498 template<class T, int N, int Multiplier>
 499 T sine_table<T,N,Multiplier>::data[N+1];
 500
 501 /// fast float to int conversion using default rounding mode
 502 inline int fastf2i_drm(float f)
 503 {
 504 #ifdef __X86__
 505     volatile int v;
 506     __asm ( "flds %1; fistpl %0" : "=m"(v) : "m"(f));
 507     return v;
 508 #else
 509     return (int)nearbyintf(f);
 510 #endif
 511 }
 512
 513 /// Convert MIDI note to frequency in Hz.
 514 inline float note_to_hz(double note, double detune_cents = 0.0)
 515 {
 516     return 440 * pow(2.0, (note - 69 + detune_cents/100.0) / 12.0);
 517 }
 518
 519 /// Hermite interpolation between two points and slopes in normalized range (written after Wikipedia article)
 520 /// @arg t normalized x coordinate (0-1 over the interval in question)
 521 /// @arg p0 first point
 522 /// @arg p1 second point
 523 /// @arg m0 first slope (multiply by interval width when using over non-1-wide interval)
 524 /// @arg m1 second slope (multiply by interval width when using over non-1-wide interval)
 525 inline float normalized_hermite(float t, float p0, float p1, float m0, float m1)
 526 {
 527     float t2 = t*t;
 528     float t3 = t2*t;
 529     return (2*t3 - 3*t2 + 1) * p0 + (t3 - 2*t2 + t) * m0 + (-2*t3 + 3*t2) * p1 + (t3-t2) * m1;
 530 }
 531
 532 /// Hermite interpolation between two points and slopes
 533 /// @arg x point within interval (x0 <= x <= x1)
 534 /// @arg x0 interval start
 535 /// @arg x1 interval end
 536 /// @arg p0 value at x0
 537 /// @arg p1 value at x1
 538 /// @arg m0 slope (steepness, tangent) at x0
 539 /// @arg m1 slope at x1
 540 inline float hermite_interpolation(float x, float x0, float x1, float p0, float p1, float m0, float m1)
 541 {
 542     float width = x1 - x0;
 543     float t = (x - x0) / width;
 544     m0 *= width;
 545     m1 *= width;
 546     float t2 = t*t;
 547     float t3 = t2*t;
 548
 549     float ct0 = p0;
 550     float ct1 = m0;
 551     float ct2 = -3 * p0 - 2 * m0 + 3 * p1 - m1;
 552     float ct3 = 2 * p0 + m0  - 2 * p1 + m1;
 553
 554     return ct3 * t3 + ct2 * t2 + ct1 * t + ct0;
 555     //return (2*t3 - 3*t2 + 1) * p0 + (t3 - 2*t2 + t) * m0 + (-2*t3 + 3*t2) * p1 + (t3-t2) * m1;
 556 }
 557
 558 /// convert amplitude value to dB
 559 inline float amp2dB(float amp)
 560 {
 561     return 20 * log10(amp);
 562 }
 563 /// convert dB to amplitude value
 564 inline float dB2amp(float db)
 565 {
 566     return exp((db / 20.0) * log(10.0));
 567 }
 568
 569 /// print binary of any data type
 570 /// assumes little endian
 571 inline void print_bits(size_t const size, void const * const ptr)
 572 {
 573     unsigned char *b = (unsigned char*) ptr;
 574     unsigned char byte;
 575     for (int i = size - 1; i >=0 ; i--) {
 576         for (int j = 7; j >= 0; j--) {
 577             byte = b[i] & (1<<j);
 578             byte >>= j;
 579             printf("%u", byte);
 580         }
 581     }
 582     puts("");
 583 }
 584
 585 };
 586 #endif