dspmath.h

   1 /*
   2 Calf Box, an open source musical instrument.
   3 Copyright (C) 2010 Krzysztof Foltman
   4
   5 This program is free software: you can redistribute it and/or modify
   6 it under the terms of the GNU General Public License as published by
   7 the Free Software Foundation, either version 3 of the License, or
   8 (at your option) any later version.
   9
  10 This program is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program.  If not, see <http://www.gnu.org/licenses/>.
  17 */
  18 #ifndef CBOX_DSPMATH_H
  19 #define CBOX_DSPMATH_H
  20
  21 #define CBOX_BLOCK_SIZE 16
  22
  23 #include <complex.h>
  24 #include <stdlib.h>
  25 #include <math.h>
  26 #include <memory.h>
  27
  28 #ifndef M_PI
  29 #include <glib.h>
  30 #define M_PI G_PI
  31 #endif
  32
  33 typedef float cbox_sample_t;
  34
  35 static inline float hz2w(float hz, float sr)
  36 {
  37     return M_PI * hz / (2 * sr);
  38 }
  39
  40 static inline float cerp_naive(float v0, float v1, float v2, float v3, float f)
  41 {
  42     float x0 = -1;
  43     float x1 = 0;
  44     float x2 = 1;
  45     float x3 = 2;
  46
  47     float l0 = ((f - x1) * (f - x2) * (f - x3)) / (            (x0 - x1) * (x0 - x2) * (x0 - x3));
  48     float l1 = ((f - x0) * (f - x2) * (f - x3)) / ((x1 - x0)             * (x1 - x2) * (x1 - x3));
  49     float l2 = ((f - x0) * (f - x1) * (f - x3)) / ((x2 - x0) * (x2 - x1)             * (x2 - x3));
  50     float l3 = ((f - x0) * (f - x1) * (f - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2)            );
  51
  52     return v0 * l0 + v1 * l1 + v2 * l2 + v3 * l3;
  53 }
  54
  55 static inline float cerp(float v0, float v1, float v2, float v3, float f)
  56 {
  57     f += 1;
  58
  59     float d0 = (f - 0);
  60     float d1 = (f - 1);
  61     float d2 = (f - 2);
  62     float d3 = (f - 3);
  63
  64     float d03 = (d0 * d3) * (1.0 / 2.0);
  65     float d12 = (d03 + 1) * (1.0 / 3.0);
  66
  67     float l0 = -d12 * d3;
  68     float l1 = d03 * d2;
  69     float l2 = -d03 * d1;
  70     float l3 = d12 * d0;
  71
  72     float y = v0 * l0 + v1 * l1 + v2 * l2 + v3 * l3;
  73     // printf("%f\n", y - cerp_naive(v0, v1, v2, v3, f - 1));
  74     return y;
  75 }
  76
  77 static inline float sanef(float v)
  78 {
  79     if (fabs(v) < (1.0 / (65536.0 * 65536.0)))
  80         return 0;
  81     return v;
  82 }
  83
  84 static inline void sanebf(float *buf)
  85 {
  86     int i;
  87     for (i = 0; i < CBOX_BLOCK_SIZE; i++)
  88         buf[i] = sanef(buf[i]);
  89 }
  90
  91 static inline void copybf(float *to, float *from)
  92 {
  93     memcpy(to, from, sizeof(float) * CBOX_BLOCK_SIZE);
  94 }
  95
  96 static inline float cent2factor(float cent)
  97 {
  98     return pow(2.0, cent / 1200.0); // I think this may be optimised using exp()
  99 }
 100
 101 static inline float dB2gain(float dB)
 102 {
 103     return pow(2.0, dB / 6.0);
 104 }
 105
 106 static inline float dB2gain_simple(float dB)
 107 {
 108     if (dB <= -96)
 109         return 0;
 110     return pow(2.0, dB / 6.0);
 111 }
 112
 113 static inline float gain2dB_simple(float gain)
 114 {
 115     if (gain < pow(2.0, -96.0 / 6.0))
 116         return -96;
 117     return 6.0 * log(gain) / log(2.0);
 118 }
 119
 120 static inline float deg2rad(float deg)
 121 {
 122     return deg * M_PI / 180;
 123 }
 124
 125 static inline float rad2deg(float rad)
 126 {
 127     return rad * 180 / M_PI;
 128 }
 129
 130 // Do a butterfly operation:
 131 // dst1 = src1 + e^iw_1*src2
 132 // dst2 = src1 + e^iw_2*src2 (w = phase * 2pi / ANALYSIS_BUFFER_SIZE)
 133 static inline void butterfly(complex float *dst1, complex float *dst2, complex float src1, complex float src2, complex float eiw1, complex float eiw2)
 134 {
 135     *dst1 = src1 + eiw1 * src2;
 136     *dst2 = src1 + eiw2 * src2;
 137 }
 138
 139 #endif