X-Git-Url: https://main.carlh.net/gitweb/?a=blobdiff_plain;f=libs%2Fardour%2Fdsp_filter.cc;h=0eef61a0f4b5a032fe2f346662e724f522ad1d3f;hb=a09fec02139770c90bdadffb56b356826c314dfd;hp=3026666b79efc7e9d77bc2a8a47f55df3cae2331;hpb=7d7f63363b13c06cd02e0f810523d897921e660a;p=ardour.git diff --git a/libs/ardour/dsp_filter.cc b/libs/ardour/dsp_filter.cc index 3026666b79..0eef61a0f4 100644 --- a/libs/ardour/dsp_filter.cc +++ b/libs/ardour/dsp_filter.cc @@ -17,10 +17,20 @@ * */ +#include #include -#include +#include +#include "ardour/dB.h" +#include "ardour/buffer.h" #include "ardour/dsp_filter.h" +#ifdef COMPILER_MSVC +#include +#define isfinite_local(val) (bool)_finite((double)val) +#else +#define isfinite_local std::isfinite +#endif + #ifndef M_PI #define M_PI 3.14159265358979323846 #endif @@ -41,6 +51,57 @@ ARDOUR::DSP::mmult (float *data, float *mult, const uint32_t n_samples) { } } +float +ARDOUR::DSP::log_meter (float power) { + // compare to gtk2_ardour/logmeter.h + static const float lower_db = -192.f; + static const float upper_db = 0.f; + static const float non_linearity = 8.0; + return (power < lower_db ? 0.0 : powf ((power - lower_db) / (upper_db - lower_db), non_linearity)); +} + +float +ARDOUR::DSP::log_meter_coeff (float coeff) { + if (coeff <= 0) return 0; + return log_meter (fast_coefficient_to_dB (coeff)); +} + +void +ARDOUR::DSP::peaks (const float *data, float &min, float &max, uint32_t n_samples) { + for (uint32_t i = 0; i < n_samples; ++i) { + if (data[i] < min) min = data[i]; + if (data[i] > max) max = data[i]; + } +} + +void +ARDOUR::DSP::process_map (BufferSet* bufs, const ChanMapping& in, const ChanMapping& out, pframes_t nframes, framecnt_t offset, const DataType& dt) +{ + const ChanMapping::Mappings& im (in.mappings()); + const ChanMapping::Mappings& om (out.mappings()); + + for (ChanMapping::Mappings::const_iterator tm = im.begin(); tm != im.end(); ++tm) { + if (tm->first != dt) { continue; } + for (ChanMapping::TypeMapping::const_iterator i = tm->second.begin(); i != tm->second.end(); ++i) { + bool valid; + const uint32_t idx = out.get (dt, i->second, &valid); + if (valid && idx != i->first) { + bufs->get (dt, idx).read_from (bufs->get (dt, i->first), nframes, offset, offset); + } + } + } + for (ChanMapping::Mappings::const_iterator tm = im.begin(); tm != im.end(); ++tm) { + if (tm->first != dt) { continue; } + for (ChanMapping::TypeMapping::const_iterator i = tm->second.begin(); i != tm->second.end(); ++i) { + bool valid; + in.get_src (dt, i->first, &valid); + if (!valid) { + bufs->get (dt, i->second).silence (nframes, offset); + } + } + } + +} LowPass::LowPass (double samplerate, float freq) : _rate (samplerate) @@ -66,6 +127,7 @@ LowPass::proc (float *data, const uint32_t n_samples) z = data[i]; } _z = z; + if (!isfinite_local (_z)) { _z = 0; } } void @@ -83,7 +145,7 @@ LowPass::ctrl (float *data, const float val, const uint32_t n_samples) /////////////////////////////////////////////////////////////////////////////// -BiQuad::BiQuad (double samplerate) +Biquad::Biquad (double samplerate) : _rate (samplerate) , _z1 (0.0) , _z2 (0.0) @@ -95,7 +157,7 @@ BiQuad::BiQuad (double samplerate) { } -BiQuad::BiQuad (const BiQuad &other) +Biquad::Biquad (const Biquad &other) : _rate (other._rate) , _z1 (0.0) , _z2 (0.0) @@ -108,7 +170,7 @@ BiQuad::BiQuad (const BiQuad &other) } void -BiQuad::run (float *data, const uint32_t n_samples) +Biquad::run (float *data, const uint32_t n_samples) { for (uint32_t i = 0; i < n_samples; ++i) { const float xn = data[i]; @@ -117,19 +179,36 @@ BiQuad::run (float *data, const uint32_t n_samples) _z2 = _b2 * xn - _a2 * z; data[i] = z; } + + if (!isfinite_local (_z1)) { _z1 = 0; } + if (!isfinite_local (_z2)) { _z2 = 0; } } void -BiQuad::compute (Type type, double freq, double Q, double gain) +Biquad::configure (double a1, double a2, double b0, double b1, double b2) { + _a1 = a1; + _a2 = a2; + _b0 = b0; + _b1 = b1; + _b2 = b2; +} + +void +Biquad::compute (Type type, double freq, double Q, double gain) +{ + if (Q <= .001) { Q = 0.001; } + if (freq <= 1.) { freq = 1.; } + if (freq >= _rate) { freq = _rate; } + /* Compute biquad filter settings. * Based on 'Cookbook formulae for audio EQ biquad filter coefficents' * by Robert Bristow-Johnson */ - const double A = pow (10.0, (gain / 40.0)); + const double A = pow (10.0, (gain / 40.0)); const double W0 = (2.0 * M_PI * freq) / _rate; - const double sinW0 = sin (W0); - const double cosW0 = cos (W0); + const double sinW0 = sin (W0); + const double cosW0 = cos (W0); const double alpha = sinW0 / (2.0 * Q); const double beta = sqrt (A) / Q; @@ -227,3 +306,121 @@ BiQuad::compute (Type type, double freq, double Q, double gain) _a1 /= _a0; _a2 /= _a0; } + +float +Biquad::dB_at_freq (float freq) const +{ + const double W0 = (2.0 * M_PI * freq) / _rate; + const float c1 = cosf (W0); + const float s1 = sinf (W0); + + const float A = _b0 + _b2; + const float B = _b0 - _b2; + const float C = 1.0 + _a2; + const float D = 1.0 - _a2; + + const float a = A * c1 + _b1; + const float b = B * s1; + const float c = C * c1 + _a1; + const float d = D * s1; + +#define SQUARE(x) ( (x) * (x) ) + float rv = 20.f * log10f (sqrtf ((SQUARE(a) + SQUARE(b)) * (SQUARE(c) + SQUARE(d))) / (SQUARE(c) + SQUARE(d))); + if (!isfinite_local (rv)) { rv = 0; } + return std::min (120.f, std::max(-120.f, rv)); +} + + +Glib::Threads::Mutex FFTSpectrum::fft_planner_lock; + +FFTSpectrum::FFTSpectrum (uint32_t window_size, double rate) + : hann_window (0) +{ + init (window_size, rate); +} + +FFTSpectrum::~FFTSpectrum () +{ + { + Glib::Threads::Mutex::Lock lk (fft_planner_lock); + fftwf_destroy_plan (_fftplan); + } + fftwf_free (_fft_data_in); + fftwf_free (_fft_data_out); + free (_fft_power); + free (hann_window); +} + +void +FFTSpectrum::init (uint32_t window_size, double rate) +{ + Glib::Threads::Mutex::Lock lk (fft_planner_lock); + + _fft_window_size = window_size; + _fft_data_size = window_size / 2; + _fft_freq_per_bin = rate / _fft_data_size / 2.f; + + _fft_data_in = (float *) fftwf_malloc (sizeof(float) * _fft_window_size); + _fft_data_out = (float *) fftwf_malloc (sizeof(float) * _fft_window_size); + _fft_power = (float *) malloc (sizeof(float) * _fft_data_size); + + reset (); + + _fftplan = fftwf_plan_r2r_1d (_fft_window_size, _fft_data_in, _fft_data_out, FFTW_R2HC, FFTW_MEASURE); + + hann_window = (float *) malloc(sizeof(float) * window_size); + double sum = 0.0; + + for (uint32_t i = 0; i < window_size; ++i) { + hann_window[i] = 0.5f - (0.5f * (float) cos (2.0f * M_PI * (float)i / (float)(window_size))); + sum += hann_window[i]; + } + const double isum = 2.0 / sum; + for (uint32_t i = 0; i < window_size; ++i) { + hann_window[i] *= isum; + } +} + +void +FFTSpectrum::reset () +{ + for (uint32_t i = 0; i < _fft_data_size; ++i) { + _fft_power[i] = 0; + } + for (uint32_t i = 0; i < _fft_window_size; ++i) { + _fft_data_out[i] = 0; + } +} + +void +FFTSpectrum::set_data_hann (float const * const data, uint32_t n_samples, uint32_t offset) +{ + assert(n_samples + offset <= _fft_window_size); + for (uint32_t i = 0; i < n_samples; ++i) { + _fft_data_in[i + offset] = data[i] * hann_window[i + offset]; + } +} + +void +FFTSpectrum::execute () +{ + fftwf_execute (_fftplan); + + _fft_power[0] = _fft_data_out[0] * _fft_data_out[0]; + +#define FRe (_fft_data_out[i]) +#define FIm (_fft_data_out[_fft_window_size - i]) + for (uint32_t i = 1; i < _fft_data_size - 1; ++i) { + _fft_power[i] = (FRe * FRe) + (FIm * FIm); + //_fft_phase[i] = atan2f (FIm, FRe); + } +#undef FRe +#undef FIm +} + +float +FFTSpectrum::power_at_bin (const uint32_t b, const float norm) const { + assert (b < _fft_data_size); + const float a = _fft_power[b] * norm; + return a > 1e-12 ? 10.0 * fast_log10 (a) : -INFINITY; +}