2 * Copyright (C) 2016 Robin Gareus <robin@gareus.org>
4 * This program is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU General Public License
6 * as published by the Free Software Foundation; either version 2
7 * of the License, or (at your option) any later version.
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 #include "ardour/dB.h"
24 #include "ardour/dsp_filter.h"
28 #define isfinite_local(val) (bool)_finite((double)val)
30 #define isfinite_local std::isfinite
34 #define M_PI 3.14159265358979323846
37 using namespace ARDOUR::DSP;
40 ARDOUR::DSP::memset (float *data, const float val, const uint32_t n_samples) {
41 for (uint32_t i = 0; i < n_samples; ++i) {
47 ARDOUR::DSP::mmult (float *data, float *mult, const uint32_t n_samples) {
48 for (uint32_t i = 0; i < n_samples; ++i) {
54 ARDOUR::DSP::log_meter (float power) {
55 // compare to gtk2_ardour/logmeter.h
56 static const float lower_db = -192.f;
57 static const float upper_db = 0.f;
58 static const float non_linearity = 8.0;
59 return (power < lower_db ? 0.0 : powf ((power - lower_db) / (upper_db - lower_db), non_linearity));
63 ARDOUR::DSP::log_meter_coeff (float coeff) {
64 if (coeff <= 0) return 0;
65 return log_meter (fast_coefficient_to_dB (coeff));
69 ARDOUR::DSP::peaks (float *data, float &min, float &max, uint32_t n_samples) {
70 for (uint32_t i = 0; i < n_samples; ++i) {
71 if (data[i] < min) min = data[i];
72 if (data[i] > max) max = data[i];
76 LowPass::LowPass (double samplerate, float freq)
84 LowPass::set_cutoff (float freq)
86 _a = 1.f - expf (-2.f * M_PI * freq / _rate);
90 LowPass::proc (float *data, const uint32_t n_samples)
95 for (uint32_t i = 0; i < n_samples; ++i) {
96 data[i] += a * (data[i] - z);
100 if (!isfinite_local (_z)) { _z = 0; }
104 LowPass::ctrl (float *data, const float val, const uint32_t n_samples)
106 // localize variables
109 for (uint32_t i = 0; i < n_samples; ++i) {
110 data[i] += a * (val - z);
116 ///////////////////////////////////////////////////////////////////////////////
118 Biquad::Biquad (double samplerate)
130 Biquad::Biquad (const Biquad &other)
131 : _rate (other._rate)
143 Biquad::run (float *data, const uint32_t n_samples)
145 for (uint32_t i = 0; i < n_samples; ++i) {
146 const float xn = data[i];
147 const float z = _b0 * xn + _z1;
148 _z1 = _b1 * xn - _a1 * z + _z2;
149 _z2 = _b2 * xn - _a2 * z;
153 if (!isfinite_local (_z1)) { _z1 = 0; }
154 if (!isfinite_local (_z2)) { _z2 = 0; }
158 Biquad::compute (Type type, double freq, double Q, double gain)
160 if (Q <= .001) { Q = 0.001; }
161 if (freq <= 1.) { freq = 1.; }
162 if (freq >= _rate) { freq = _rate; }
164 /* Compute biquad filter settings.
165 * Based on 'Cookbook formulae for audio EQ biquad filter coefficents'
166 * by Robert Bristow-Johnson
168 const double A = pow (10.0, (gain / 40.0));
169 const double W0 = (2.0 * M_PI * freq) / _rate;
170 const double sinW0 = sin (W0);
171 const double cosW0 = cos (W0);
172 const double alpha = sinW0 / (2.0 * Q);
173 const double beta = sqrt (A) / Q;
179 _b0 = (1.0 - cosW0) / 2.0;
181 _b2 = (1.0 - cosW0) / 2.0;
188 _b0 = (1.0 + cosW0) / 2.0;
189 _b1 = -(1.0 + cosW0);
190 _b2 = (1.0 + cosW0) / 2.0;
196 case BandPassSkirt: /* Constant skirt gain, peak gain = Q */
205 case BandPass0dB: /* Constant 0 dB peak gain */
233 _b0 = 1.0 + (alpha * A);
235 _b2 = 1.0 - (alpha * A);
236 _a0 = 1.0 + (alpha / A);
238 _a2 = 1.0 - (alpha / A);
242 _b0 = A * ((A + 1) - ((A - 1) * cosW0) + (beta * sinW0));
243 _b1 = (2.0 * A) * ((A - 1) - ((A + 1) * cosW0));
244 _b2 = A * ((A + 1) - ((A - 1) * cosW0) - (beta * sinW0));
245 _a0 = (A + 1) + ((A - 1) * cosW0) + (beta * sinW0);
246 _a1 = -2.0 * ((A - 1) + ((A + 1) * cosW0));
247 _a2 = (A + 1) + ((A - 1) * cosW0) - (beta * sinW0);
251 _b0 = A * ((A + 1) + ((A - 1) * cosW0) + (beta * sinW0));
252 _b1 = -(2.0 * A) * ((A - 1) + ((A + 1) * cosW0));
253 _b2 = A * ((A + 1) + ((A - 1) * cosW0) - (beta * sinW0));
254 _a0 = (A + 1) - ((A - 1) * cosW0) + (beta * sinW0);
255 _a1 = 2.0 * ((A - 1) - ((A + 1) * cosW0));
256 _a2 = (A + 1) - ((A - 1) * cosW0) - (beta * sinW0);
259 abort(); /*NOTREACHED*/
271 Biquad::dB_at_freq (float freq) const
273 const double W0 = (2.0 * M_PI * freq) / _rate;
274 const float c1 = cosf (W0);
275 const float s1 = sinf (W0);
277 const float A = _b0 + _b2;
278 const float B = _b0 - _b2;
279 const float C = 1.0 + _a2;
280 const float D = 1.0 - _a2;
282 const float a = A * c1 + _b1;
283 const float b = B * s1;
284 const float c = C * c1 + _a1;
285 const float d = D * s1;
287 #define SQUARE(x) ( (x) * (x) )
288 float rv = 20.f * log10f (sqrtf ((SQUARE(a) + SQUARE(b)) * (SQUARE(c) + SQUARE(d))) / (SQUARE(c) + SQUARE(d)));
289 if (!isfinite_local (rv)) { rv = 0; }
290 return std::min (120.f, std::max(-120.f, rv));