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 struct cbox_sincos
  36 {
  37     float sine;
  38     float cosine;
  39     float prewarp;
  40 };
  41
  42 static inline float hz2w(float hz, float sr)
  43 {
  44     return M_PI * hz / (2 * sr);
  45 }
  46
  47 static inline float cerp_naive(float v0, float v1, float v2, float v3, float f)
  48 {
  49     float x0 = -1;
  50     float x1 = 0;
  51     float x2 = 1;
  52     float x3 = 2;
  53
  54     float l0 = ((f - x1) * (f - x2) * (f - x3)) / (            (x0 - x1) * (x0 - x2) * (x0 - x3));
  55     float l1 = ((f - x0) * (f - x2) * (f - x3)) / ((x1 - x0)             * (x1 - x2) * (x1 - x3));
  56     float l2 = ((f - x0) * (f - x1) * (f - x3)) / ((x2 - x0) * (x2 - x1)             * (x2 - x3));
  57     float l3 = ((f - x0) * (f - x1) * (f - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2)            );
  58
  59     return v0 * l0 + v1 * l1 + v2 * l2 + v3 * l3;
  60 }
  61
  62 static inline float cerp(float v0, float v1, float v2, float v3, float f)
  63 {
  64     f += 1;
  65
  66     float d0 = (f - 0);
  67     float d1 = (f - 1);
  68     float d2 = (f - 2);
  69     float d3 = (f - 3);
  70
  71     float d03 = (d0 * d3) * (1.0 / 2.0);
  72     float d12 = (d03 + 1) * (1.0 / 3.0);
  73
  74     float l0 = -d12 * d3;
  75     float l1 = d03 * d2;
  76     float l2 = -d03 * d1;
  77     float l3 = d12 * d0;
  78
  79     float y = v0 * l0 + v1 * l1 + v2 * l2 + v3 * l3;
  80     // printf("%f\n", y - cerp_naive(v0, v1, v2, v3, f - 1));
  81     return y;
  82 }
  83
  84 static inline float sanef(float v)
  85 {
  86     if (fabs(v) < (1.0 / (65536.0 * 65536.0)))
  87         return 0;
  88     return v;
  89 }
  90
  91 static inline void sanebf(float *buf)
  92 {
  93     int i;
  94     for (i = 0; i < CBOX_BLOCK_SIZE; i++)
  95         buf[i] = sanef(buf[i]);
  96 }
  97
  98 static inline void copybf(float *to, float *from)
  99 {
 100     memcpy(to, from, sizeof(float) * CBOX_BLOCK_SIZE);
 101 }
 102
 103 static inline float cent2factor(float cent)
 104 {
 105     return powf(2.0, cent * (1.f / 1200.f)); // I think this may be optimised using exp()
 106 }
 107
 108 static inline float dB2gain(float dB)
 109 {
 110     return powf(2.f, dB * (1.f / 6.f));
 111 }
 112
 113 static inline float dB2gain_simple(float dB)
 114 {
 115     if (dB <= -96)
 116         return 0;
 117     return powf(2.f, dB * (1.f / 6.f));
 118 }
 119
 120 static inline float gain2dB_simple(float gain)
 121 {
 122     static const float sixoverlog2 = 8.656170245333781; // 6.0 / logf(2.f);
 123     if (gain < (1.f / 65536.f))
 124         return -96.f;
 125     return sixoverlog2 * logf(gain);
 126 }
 127
 128 static inline float deg2rad(float deg)
 129 {
 130     return deg * (float)(M_PI / 180.f);
 131 }
 132
 133 static inline float rad2deg(float rad)
 134 {
 135     return rad * (float)(180.f / M_PI);
 136 }
 137
 138 // Do a butterfly operation:
 139 // dst1 = src1 + e^iw_1*src2
 140 // dst2 = src1 + e^iw_2*src2 (w = phase * 2pi / ANALYSIS_BUFFER_SIZE)
 141 static inline void butterfly(complex float *dst1, complex float *dst2, complex float src1, complex float src2, complex float eiw1, complex float eiw2)
 142 {
 143     *dst1 = src1 + eiw1 * src2;
 144     *dst2 = src1 + eiw2 * src2;
 145 }
 146
 147 #endif