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 <stack>
  25 #include <map>
  26 #include <cmath>
  27 #include <cstdlib>
  28 #include <stdint.h>
  29 #include <stdio.h>
  30
  31 namespace dsp {
  32
  33 /// Set a float to zero
  34 inline void zero(float &v) {
  35     v = 0;
  36 };
  37
  38 /// Set a double to zero
  39 inline void zero(double &v) {
  40     v = 0;
  41 };
  42
  43 /// Set 64-bit unsigned integer value to zero
  44 inline void zero(uint64_t &v) { v = 0; };
  45 /// Set 32-bit unsigned integer value to zero
  46 inline void zero(uint32_t &v) { v = 0; };
  47 /// Set 16-bit unsigned integer value to zero
  48 inline void zero(uint16_t &v) { v = 0; };
  49 /// Set 8-bit unsigned integer value to zero
  50 inline void zero(uint8_t &v) { v = 0; };
  51 /// Set 64-bit signed integer value to zero
  52 inline void zero(int64_t &v) { v = 0; };
  53 /// Set 32-bit signed integer value to zero
  54 inline void zero(int32_t &v) { v = 0; };
  55 /// Set 16-bit signed integer value to zero
  56 inline void zero(int16_t &v) { v = 0; };
  57 /// Set 8-bit signed integer value to zero
  58 inline void zero(int8_t &v) { v = 0; };
  59
  60 /// Set array (buffer or anything similar) to vector of zeroes
  61 template<class T>
  62 void zero(T *data, unsigned int size) {
  63     T value;
  64     dsp::zero(value);
  65     for (unsigned int i=0; i<size; i++)
  66         *data++ = value;
  67 }
  68
  69 template<class T = float>struct stereo_sample {
  70     T left;
  71     T right;
  72     /// default constructor - preserves T's semantics (ie. no implicit initialization to 0)
  73     inline stereo_sample() {
  74     }
  75     inline stereo_sample(T _left, T _right) {
  76         left = _left;
  77         right = _right;
  78     }
  79     inline stereo_sample(T _both) {
  80         left = right = _both;
  81     }
  82     template<typename U>
  83     inline stereo_sample(const stereo_sample<U> &value) {
  84         left = value.left;
  85         right = value.right;
  86     }
  87     inline stereo_sample& operator=(const T &value) {
  88         left = right = value;
  89         return *this;
  90     }
  91     template<typename U>
  92     inline stereo_sample& operator=(const stereo_sample<U> &value) {
  93         left = value.left;
  94         right = value.right;
  95         return *this;
  96     }
  97 /*
  98     inline operator T() const {
  99         return (left+right)/2;
 100     }
 101 */
 102     inline stereo_sample& operator*=(const T &multiplier) {
 103         left *= multiplier;
 104         right *= multiplier;
 105         return *this;
 106     }
 107     inline stereo_sample& operator+=(const stereo_sample<T> &value) {
 108         left += value.left;
 109         right += value.right;
 110         return *this;
 111     }
 112     inline stereo_sample& operator-=(const stereo_sample<T> &value) {
 113         left -= value.left;
 114         right -= value.right;
 115         return *this;
 116     }
 117     template<typename U> inline stereo_sample<U> operator*(const U &value) const {
 118         return stereo_sample<U>(left*value, right*value);
 119     }
 120     /*inline stereo_sample<float> operator*(float value) const {
 121         return stereo_sample<float>(left*value, right*value);
 122     }
 123     inline stereo_sample<double> operator*(double value) const {
 124         return stereo_sample<double>(left*value, right*value);
 125     }*/
 126     inline stereo_sample<T> operator+(const stereo_sample<T> &value) {
 127         return stereo_sample(left+value.left, right+value.right);
 128     }
 129     inline stereo_sample<T> operator-(const stereo_sample<T> &value) {
 130         return stereo_sample(left-value.left, right-value.right);
 131     }
 132     inline stereo_sample<T> operator+(const T &value) {
 133         return stereo_sample(left+value, right+value);
 134     }
 135     inline stereo_sample<T> operator-(const T &value) {
 136         return stereo_sample(left-value, right-value);
 137     }
 138     inline stereo_sample<float> operator+(float value) {
 139         return stereo_sample<float>(left+value, right+value);
 140     }
 141     inline stereo_sample<float> operator-(float value) {
 142         return stereo_sample<float>(left-value, right-value);
 143     }
 144     inline stereo_sample<double> operator+(double value) {
 145         return stereo_sample<double>(left+value, right+value);
 146     }
 147     inline stereo_sample<double> operator-(double value) {
 148         return stereo_sample<double>(left-value, right-value);
 149     }
 150 };
 151
 152 /// Multiply constant by stereo_value
 153 template<class T>
 154 inline stereo_sample<T> operator*(const T &value, const stereo_sample<T> &value2) {
 155     return stereo_sample<T>(value2.left*value, value2.right*value);
 156 }
 157
 158 /// Add constant to stereo_value
 159 template<class T>
 160 inline stereo_sample<T> operator+(const T &value, const stereo_sample<T> &value2) {
 161     return stereo_sample<T>(value2.left+value, value2.right+value);
 162 }
 163
 164 /// Subtract stereo_value from constant (yields stereo_value of course)
 165 template<class T>
 166 inline stereo_sample<T> operator-(const T &value, const stereo_sample<T> &value2) {
 167     return stereo_sample<T>(value-value2.left, value-value2.right);
 168 }
 169
 170 /// Shift value right by 'bits' bits (multiply by 2^-bits)
 171 template<typename T>
 172 inline stereo_sample<T> shr(stereo_sample<T> v, int bits = 1) {
 173     v.left = shr(v.left, bits);
 174     v.right = shr(v.right, bits);
 175     return v;
 176 }
 177
 178 /// Set a stereo_sample<T> value to zero
 179 template<typename T>
 180 inline void zero(stereo_sample<T> &v) {
 181     dsp::zero(v.left);
 182     dsp::zero(v.right);
 183 }
 184
 185 /// 'Small value' for integer and other types
 186 template<typename T>
 187 inline T small_value() {
 188     return 0;
 189 }
 190
 191 /// '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).
 192 template<>
 193 inline float small_value<float>() {
 194     return (1.0/16777216.0); // allows for 2^-24, should be enough for 24-bit DACs at least :)
 195 }
 196
 197 /// 'Small value' for doubles (2^-24) - used for primitive underrun prevention. The value is pretty much arbitrary.
 198 template<>
 199 inline double small_value<double>() {
 200     return (1.0/16777216.0);
 201 }
 202
 203 /// Convert a single value to single value = do nothing :) (but it's a generic with specialisation for stereo_sample)
 204 template<typename T>
 205 inline float mono(T v) {
 206     return v;
 207 }
 208
 209 /// Convert a stereo_sample to single value by averaging two channels
 210 template<typename T>
 211 inline T mono(stereo_sample<T> v) {
 212     return shr(v.left+v.right);
 213 }
 214
 215 /// Clip a value to [min, max]
 216 template<typename T>
 217 inline T clip(T value, T min, T max) {
 218     if (value < min) return min;
 219     if (value > max) return max;
 220     return value;
 221 }
 222
 223 /// Clip a double to [-1.0, +1.0]
 224 inline double clip11(double value) {
 225     double a = fabs(value);
 226     if (a<=1) return value;
 227     return (a<0) ? -1.0 : 1.0;
 228 }
 229
 230 /// Clip a float to [-1.0f, +1.0f]
 231 inline float clip11(float value) {
 232     float a = fabsf(value);
 233     if (a<=1) return value;
 234     return (a<0) ? -1.0f : 1.0f;
 235 }
 236
 237 /// Clip a double to [0.0, +1.0]
 238 inline double clip01(double value) {
 239     double a = fabs(value-0.5);
 240     if (a<=0.5) return value;
 241     return (a<0) ? -0.0 : 1.0;
 242 }
 243
 244 /// Clip a float to [0.0f, +1.0f]
 245 inline float clip01(float value) {
 246     float a = fabsf(value-0.5f);
 247     if (a<=0.5f) return value;
 248     return (a<0) ? -0.0f : 1.0f;
 249 }
 250
 251 // Linear interpolation (mix-way between v1 and v2).
 252 template<typename T, typename U>
 253 inline T lerp(T v1, T v2, U mix) {
 254     return v1+(v2-v1)*mix;
 255 }
 256
 257 // Linear interpolation for stereo values (mix-way between v1 and v2).
 258 template<typename T>
 259 inline stereo_sample<T> lerp(stereo_sample<T> &v1, stereo_sample<T> &v2, float mix) {
 260     return stereo_sample<T>(v1.left+(v2.left-v1.left)*mix, v1.right+(v2.right-v1.right)*mix);
 261 }
 262
 263 /**
 264  * decay-only envelope (linear or exponential); deactivates itself when it goes below a set point (epsilon)
 265  */
 266 class decay
 267 {
 268     double value, initial;
 269     unsigned int age, mask;
 270     bool active;
 271 public:
 272     decay() {
 273         active = false;
 274         mask = 127;
 275         initial = value = 0.0;
 276     }
 277     inline bool get_active() {
 278         return active;
 279     }
 280     inline double get() {
 281         return active ? value : 0.0;
 282     }
 283     inline void set(double v) {
 284         initial = value = v;
 285         active = true;
 286         age = 0;
 287     }
 288     /// reinitialise envelope (must be called if shape changes from linear to exponential or vice versa in the middle of envelope)
 289     inline void reinit()
 290     {
 291         initial = value;
 292         age = 1;
 293     }
 294     inline void add(double v) {
 295         if (active)
 296             value += v;
 297         else
 298             value = v;
 299         initial = value;
 300         age = 0;
 301         active = true;
 302     }
 303     static inline double calc_exp_constant(double times, double cycles)
 304     {
 305         if (cycles < 1.0)
 306             cycles = 1.0;
 307         return pow(times, 1.0 / cycles);
 308     }
 309     inline void age_exp(double constant, double epsilon) {
 310         if (active) {
 311             if (!(age & mask))
 312                 value = initial * pow(constant, (double)age);
 313             else
 314                 value *= constant;
 315             if (value < epsilon)
 316                 active = false;
 317             age++;
 318         }
 319     }
 320     inline void age_lin(double constant, double epsilon) {
 321         if (active) {
 322             if (!(age & mask))
 323                 value = initial - constant * age;
 324             else
 325                 value -= constant;
 326             if (value < epsilon)
 327                 active = false;
 328             age++;
 329         }
 330     }
 331     inline void deactivate() {
 332         active = false;
 333         value = 0;
 334     }
 335 };
 336
 337 class scheduler;
 338
 339 class task {
 340 public:
 341     virtual void execute(scheduler *s)=0;
 342     virtual void dispose() { delete this; }
 343     virtual ~task() {}
 344 };
 345
 346 /// this scheduler is based on std::multimap, so it isn't very fast, I guess
 347 /// maybe some day it should be rewritten to use heapsort or something
 348 /// work in progress, don't use!
 349 class scheduler {
 350     std::multimap<unsigned int, task *> timeline;
 351     unsigned int time, next_task;
 352     bool eob;
 353     class end_buf_task: public task {
 354     public:
 355         scheduler *p;
 356         end_buf_task(scheduler *_p) : p(_p) {}
 357         virtual void execute(scheduler *s) { p->eob = true; }
 358         virtual void dispose() { }
 359     } eobt;
 360 public:
 361
 362     scheduler()
 363     : time(0)
 364     , next_task((unsigned)-1)
 365     , eob(true)
 366     , eobt (this)
 367     {
 368         time = 0;
 369         next_task = (unsigned)-1;
 370         eob = false;
 371     }
 372     inline bool is_next_tick() {
 373         if (time < next_task)
 374             return true;
 375         do_tasks();
 376     }
 377     inline void next_tick() {
 378         time++;
 379     }
 380     void set(int pos, task *t) {
 381         timeline.insert(std::pair<unsigned int, task *>(time+pos, t));
 382         next_task = timeline.begin()->first;
 383     }
 384     void do_tasks() {
 385         std::multimap<unsigned int, task *>::iterator i = timeline.begin();
 386         while(i != timeline.end() && i->first == time) {
 387             i->second->execute(this);
 388             i->second->dispose();
 389             timeline.erase(i);
 390         }
 391     }
 392     bool is_eob() {
 393         return eob;
 394     }
 395     void set_buffer_size(int count) {
 396         set(count, &eobt);
 397     }
 398 };
 399
 400 /**
 401  * Force "small enough" float value to zero
 402  */
 403 inline void sanitize(float &value)
 404 {
 405     if (std::abs(value) < small_value<float>())
 406         value = 0.f;
 407 }
 408
 409 /**
 410  * Force "small enough" double value to zero
 411  */
 412 inline void sanitize(double &value)
 413 {
 414     if (std::abs(value) < small_value<double>())
 415         value = 0.f;
 416 }
 417
 418 /**
 419  * Force "small enough" stereo value to zero
 420  */
 421 template<class T>
 422 inline void sanitize(stereo_sample<T> &value)
 423 {
 424     sanitize(value.left);
 425     sanitize(value.right);
 426 }
 427
 428 inline float fract16(unsigned int value)
 429 {
 430     return (value & 0xFFFF) * (1.0 / 65536.0);
 431 }
 432
 433 /**
 434  * typical precalculated sine table
 435  */
 436 template<class T, int N, int Multiplier>
 437 class sine_table
 438 {
 439 public:
 440     static bool initialized;
 441     static T data[N+1];
 442     sine_table() {
 443         if (initialized)
 444             return;
 445         initialized = true;
 446         for (int i=0; i<N+1; i++)
 447             data[i] = (T)(Multiplier*sin(i*2*M_PI*(1.0/N)));
 448     }
 449 };
 450
 451 template<class T, int N, int Multiplier>
 452 bool sine_table<T,N,Multiplier>::initialized = false;
 453
 454 template<class T, int N, int Multiplier>
 455 T sine_table<T,N,Multiplier>::data[N+1];
 456
 457 /// fast float to int conversion using default rounding mode
 458 inline int fastf2i_drm(float f)
 459 {
 460 #ifdef __X86__
 461     volatile int v;
 462     __asm ( "flds %1; fistpl %0" : "=m"(v) : "m"(f));
 463     return v;
 464 #else
 465     return (int)nearbyintf(f);
 466 #endif
 467 }
 468
 469 /// Convert MIDI note to frequency in Hz.
 470 inline float note_to_hz(double note, double detune_cents = 0.0)
 471 {
 472     return 440 * pow(2.0, (note - 69 + detune_cents/100.0) / 12.0);
 473 }
 474
 475 /// Hermite interpolation between two points and slopes in normalized range (written after Wikipedia article)
 476 /// @arg t normalized x coordinate (0-1 over the interval in question)
 477 /// @arg p0 first point
 478 /// @arg p1 second point
 479 /// @arg m0 first slope (multiply by interval width when using over non-1-wide interval)
 480 /// @arg m1 second slope (multiply by interval width when using over non-1-wide interval)
 481 inline float normalized_hermite(float t, float p0, float p1, float m0, float m1)
 482 {
 483     float t2 = t*t;
 484     float t3 = t2*t;
 485     return (2*t3 - 3*t2 + 1) * p0 + (t3 - 2*t2 + t) * m0 + (-2*t3 + 3*t2) * p1 + (t3-t2) * m1;
 486 }
 487
 488 /// Hermite interpolation between two points and slopes
 489 /// @arg x point within interval (x0 <= x <= x1)
 490 /// @arg x0 interval start
 491 /// @arg x1 interval end
 492 /// @arg p0 value at x0
 493 /// @arg p1 value at x1
 494 /// @arg m0 slope (steepness, tangent) at x0
 495 /// @arg m1 slope at x1
 496 inline float hermite_interpolation(float x, float x0, float x1, float p0, float p1, float m0, float m1)
 497 {
 498     float width = x1 - x0;
 499     float t = (x - x0) / width;
 500     m0 *= width;
 501     m1 *= width;
 502     float t2 = t*t;
 503     float t3 = t2*t;
 504
 505     float ct0 = p0;
 506     float ct1 = m0;
 507     float ct2 = -3 * p0 - 2 * m0 + 3 * p1 - m1;
 508     float ct3 = 2 * p0 + m0  - 2 * p1 + m1;
 509
 510     return ct3 * t3 + ct2 * t2 + ct1 * t + ct0;
 511     //return (2*t3 - 3*t2 + 1) * p0 + (t3 - 2*t2 + t) * m0 + (-2*t3 + 3*t2) * p1 + (t3-t2) * m1;
 512 }
 513
 514 };
 515
 516 #endif