add a const - just because.
[ardour.git] / libs / ardour / dsp_filter.cc
index 3026666b79efc7e9d77bc2a8a47f55df3cae2331..0eef61a0f4b5a032fe2f346662e724f522ad1d3f 100644 (file)
  *
  */
 
+#include <algorithm>
 #include <stdlib.h>
-#include <math.h>
+#include <cmath>
+#include "ardour/dB.h"
+#include "ardour/buffer.h"
 #include "ardour/dsp_filter.h"
 
+#ifdef COMPILER_MSVC
+#include <float.h>
+#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;
+}