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.
22 #include "ardour/dB.h"
23 #include "ardour/dsp_filter.h"
26 #define M_PI 3.14159265358979323846
29 using namespace ARDOUR::DSP;
32 ARDOUR::DSP::memset (float *data, const float val, const uint32_t n_samples) {
33 for (uint32_t i = 0; i < n_samples; ++i) {
39 ARDOUR::DSP::mmult (float *data, float *mult, const uint32_t n_samples) {
40 for (uint32_t i = 0; i < n_samples; ++i) {
46 ARDOUR::DSP::log_meter (float power) {
47 // compare to gtk2_ardour/logmeter.h
48 static const float lower_db = -192.f;
49 static const float upper_db = 0.f;
50 static const float non_linearity = 8.0;
51 return (power < lower_db ? 0.0 : powf ((power - lower_db) / (upper_db - lower_db), non_linearity));
55 ARDOUR::DSP::log_meter_coeff (float coeff) {
56 if (coeff <= 0) return 0;
57 return log_meter (fast_coefficient_to_dB (coeff));
61 ARDOUR::DSP::peaks (float *data, float &min, float &max, uint32_t n_samples) {
62 for (uint32_t i = 0; i < n_samples; ++i) {
63 if (data[i] < min) min = data[i];
64 if (data[i] > max) max = data[i];
68 LowPass::LowPass (double samplerate, float freq)
76 LowPass::set_cutoff (float freq)
78 _a = 1.f - expf (-2.f * M_PI * freq / _rate);
82 LowPass::proc (float *data, const uint32_t n_samples)
87 for (uint32_t i = 0; i < n_samples; ++i) {
88 data[i] += a * (data[i] - z);
95 LowPass::ctrl (float *data, const float val, const uint32_t n_samples)
100 for (uint32_t i = 0; i < n_samples; ++i) {
101 data[i] += a * (val - z);
107 ///////////////////////////////////////////////////////////////////////////////
109 BiQuad::BiQuad (double samplerate)
121 BiQuad::BiQuad (const BiQuad &other)
122 : _rate (other._rate)
134 BiQuad::run (float *data, const uint32_t n_samples)
136 for (uint32_t i = 0; i < n_samples; ++i) {
137 const float xn = data[i];
138 const float z = _b0 * xn + _z1;
139 _z1 = _b1 * xn - _a1 * z + _z2;
140 _z2 = _b2 * xn - _a2 * z;
146 BiQuad::compute (Type type, double freq, double Q, double gain)
148 /* Compute biquad filter settings.
149 * Based on 'Cookbook formulae for audio EQ biquad filter coefficents'
150 * by Robert Bristow-Johnson
152 const double A = pow (10.0, (gain / 40.0));
153 const double W0 = (2.0 * M_PI * freq) / _rate;
154 const double sinW0 = sin (W0);
155 const double cosW0 = cos (W0);
156 const double alpha = sinW0 / (2.0 * Q);
157 const double beta = sqrt (A) / Q;
163 _b0 = (1.0 - cosW0) / 2.0;
165 _b2 = (1.0 - cosW0) / 2.0;
172 _b0 = (1.0 + cosW0) / 2.0;
173 _b1 = -(1.0 + cosW0);
174 _b2 = (1.0 + cosW0) / 2.0;
180 case BandPassSkirt: /* Constant skirt gain, peak gain = Q */
189 case BandPass0dB: /* Constant 0 dB peak gain */
217 _b0 = 1.0 + (alpha * A);
219 _b2 = 1.0 - (alpha * A);
220 _a0 = 1.0 + (alpha / A);
222 _a2 = 1.0 - (alpha / A);
226 _b0 = A * ((A + 1) - ((A - 1) * cosW0) + (beta * sinW0));
227 _b1 = (2.0 * A) * ((A - 1) - ((A + 1) * cosW0));
228 _b2 = A * ((A + 1) - ((A - 1) * cosW0) - (beta * sinW0));
229 _a0 = (A + 1) + ((A - 1) * cosW0) + (beta * sinW0);
230 _a1 = -2.0 * ((A - 1) + ((A + 1) * cosW0));
231 _a2 = (A + 1) + ((A - 1) * cosW0) - (beta * sinW0);
235 _b0 = A * ((A + 1) + ((A - 1) * cosW0) + (beta * sinW0));
236 _b1 = -(2.0 * A) * ((A - 1) + ((A + 1) * cosW0));
237 _b2 = A * ((A + 1) + ((A - 1) * cosW0) - (beta * sinW0));
238 _a0 = (A + 1) - ((A - 1) * cosW0) + (beta * sinW0);
239 _a1 = 2.0 * ((A - 1) - ((A + 1) * cosW0));
240 _a2 = (A + 1) - ((A - 1) * cosW0) - (beta * sinW0);
243 abort(); /*NOTREACHED*/