1 // Licensed GNU LGPL v3 or later: http://www.gnu.org/licenses/lgpl.html
2
3 #include "smmath.hh"
4 #include <assert.h>
5
6 namespace SpectMorph {
7
8 double
db_to_factor(double dB)9 db_to_factor (double dB)
10 {
11 double factor = dB / 20; /* Bell */
12 return pow (10, factor);
13 }
14
15 double
db_from_factor(double factor,double min_dB)16 db_from_factor (double factor, double min_dB)
17 {
18 if (factor > 0)
19 {
20 double dB = log10 (factor); /* Bell */
21 dB *= 20;
22 return dB;
23 }
24 else
25 return min_dB;
26 }
27
28 int
sm_factor2delta_idb(double factor)29 sm_factor2delta_idb (double factor)
30 {
31 return int (sm_factor2idb (factor)) - (512 * 64);
32 }
33
34 double
sm_idb2factor_slow(uint16_t idb)35 sm_idb2factor_slow (uint16_t idb)
36 {
37 double db = idb / 64.0 - 512;
38 return db_to_factor (db);
39 }
40
41 #define FAC 6000.0
42 #define ADD (3 * FAC)
43
44 uint16_t
sm_freq2ifreq(double freq)45 sm_freq2ifreq (double freq)
46 {
47 return sm_bound (0, sm_round_positive (log (freq) * FAC + ADD), 65535);
48 }
49
50 double
sm_ifreq2freq_slow(uint16_t ifreq)51 sm_ifreq2freq_slow (uint16_t ifreq)
52 {
53 return exp ((ifreq - ADD) / FAC);
54 }
55
56 /* tables for:
57 *
58 * - fast idb -> factor conversion
59 * - fast ifreq -> freq conversion
60 *
61 * exp (high + low) = exp (high) * exp (low)
62 */
63 float MathTables::idb2f_high[256];
64 float MathTables::idb2f_low[256];
65
66 float MathTables::ifreq2f_high[256];
67 float MathTables::ifreq2f_low[256];
68
69 void
sm_math_init()70 sm_math_init()
71 {
72 for (size_t i = 0; i < 256; i++)
73 {
74 MathTables::idb2f_high[i] = sm_idb2factor_slow (i * 256);
75 MathTables::idb2f_low[i] = sm_idb2factor_slow (64 * 512 + i);
76
77 MathTables::ifreq2f_high[i] = sm_ifreq2freq_slow (i * 256);
78 MathTables::ifreq2f_low[i] = sm_ifreq2freq_slow (ADD + i);
79 }
80
81 // ensure proper rounding mode
82 assert (sm_fpu_okround());
83
84 assert (sm_round_positive (42.51) == 43);
85 assert (sm_round_positive (3.14) == 3);
86 assert (sm_round_positive (2.1) == 2);
87 assert (sm_round_positive (0.7) == 1);
88 assert (sm_round_positive (0.2) == 0);
89 }
90
91 int
sm_fpu_okround()92 sm_fpu_okround()
93 {
94 typedef unsigned short int BseFpuState;
95
96 BseFpuState cv;
97 __asm__ ("fnstcw %0"
98 : "=m" (*&cv));
99 return !(cv & 0x0c00);
100 }
101
102 double
sm_lowpass1_factor(double mix_freq,double freq)103 sm_lowpass1_factor (double mix_freq, double freq)
104 {
105 const double delta_t = (1 / mix_freq);
106
107 return 2 * M_PI * delta_t * freq / (2 * M_PI * delta_t * freq + 1);
108 }
109
110 double
sm_xparam(double x,double slope)111 sm_xparam (double x, double slope)
112 {
113 return pow (x, slope);
114 }
115
116 double
sm_xparam_inv(double x,double slope)117 sm_xparam_inv (double x, double slope)
118 {
119 return pow (x, 1 / slope);
120 }
121
122 double
sm_bessel_i0(double x)123 sm_bessel_i0 (double x)
124 {
125 /* http://www.vibrationdata.com/Bessel.htm */
126
127 /* 1 + (x/2)^2/(1!^2)
128 * + (x/2)^4/(2!^2)
129 * + (x/2)^6/(3!^2) ... */
130
131 double delta = 1;
132 double result = 1;
133 const double sqr_x_2 = (x/2)*(x/2);
134
135 for (int i = 1; i < 500; i++)
136 {
137 delta *= sqr_x_2 / (i * i);
138 result += delta;
139
140 if (delta < 1e-14 * result)
141 break;
142 }
143 return result;
144 }
145
146 double
velocity_to_gain(double velocity,double vrange_db)147 velocity_to_gain (double velocity, double vrange_db)
148 {
149 /* vrange_db should be positive or zero */
150 g_return_val_if_fail (vrange_db > -0.01, 0);
151
152 /* convert, so that
153 * - gain (0) is at the -vrange_db
154 * - gain (1) is 0 db
155 * - sqrt(gain(v)) is a straight line
156 *
157 * See Roger B. Dannenberg: The Interpretation of Midi Velocity
158 */
159 const double b = db_to_factor (-vrange_db * 0.5);
160 const double x = (1 - b) * velocity + b;
161
162 return x * x;
163 }
164
165 }
166