X-Git-Url: https://main.carlh.net/gitweb/?a=blobdiff_plain;f=libs%2Fardour%2Fdsp_filter.cc;h=3508bf4e290796ea3ed1a8e316246ed12bcad4fc;hb=84272b4e27e537bf2c38c9cd25675c61addea40a;hp=1cf87cf3f9018a50651d391ba3987c2747d2721c;hpb=53c188beb3de47f903482ac3c7812d756e94508f;p=ardour.git

diff --git a/libs/ardour/dsp_filter.cc b/libs/ardour/dsp_filter.cc
index 1cf87cf3f9..3508bf4e29 100644
--- a/libs/ardour/dsp_filter.cc
+++ b/libs/ardour/dsp_filter.cc
@@ -21,6 +21,7 @@
 #include <stdlib.h>
 #include <cmath>
 #include "ardour/dB.h"
+#include "ardour/buffer.h"
 #include "ardour/dsp_filter.h"
 
 #ifdef COMPILER_MSVC
@@ -52,7 +53,7 @@ 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
+	// compare to libs/ardour/log_meter.h
 	static const float lower_db = -192.f;
 	static const float upper_db = 0.f;
 	static const float non_linearity = 8.0;
@@ -66,13 +67,41 @@ ARDOUR::DSP::log_meter_coeff (float coeff) {
 }
 
 void
-ARDOUR::DSP::peaks (float *data, float &min, float &max, uint32_t n_samples) {
+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, samplecnt_t offset, const DataType& dt)
+{
+	const ChanMapping::Mappings& im (in.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)
 	, _z (0)
@@ -115,7 +144,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)
@@ -127,7 +156,7 @@ BiQuad::BiQuad (double samplerate)
 {
 }
 
-BiQuad::BiQuad (const BiQuad &other)
+Biquad::Biquad (const Biquad &other)
 	: _rate (other._rate)
 	, _z1 (0.0)
 	, _z2 (0.0)
@@ -140,7 +169,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];
@@ -155,11 +184,21 @@ BiQuad::run (float *data, const uint32_t n_samples)
 }
 
 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; }
+	if (Q <= .001)  { Q = 0.001; }
+	if (freq <= 1.) { freq = 1.; }
+	if (freq >= 0.4998 * _rate) { freq = 0.4998 * _rate; }
 
 	/* Compute biquad filter settings.
 	 * Based on 'Cookbook formulae for audio EQ biquad filter coefficents'
@@ -268,7 +307,7 @@ BiQuad::compute (Type type, double freq, double Q, double gain)
 }
 
 float
-BiQuad::dB_at_freq (float freq) const
+Biquad::dB_at_freq (float freq) const
 {
 	const double W0 = (2.0 * M_PI * freq) / _rate;
 	const float c1 = cosf (W0);
@@ -289,3 +328,181 @@ BiQuad::dB_at_freq (float freq) const
 	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)
+{
+	assert (window_size > 0);
+	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;
+}
+
+Generator::Generator ()
+	: _type (UniformWhiteNoise)
+	, _rseed (1)
+{
+	set_type (UniformWhiteNoise);
+}
+
+void
+Generator::set_type (Generator::Type t) {
+	_type = t;
+	_b0 = _b1 = _b2 = _b3 = _b4 = _b5 = _b6 = 0;
+	_pass = false;
+	_rn = 0;
+}
+
+void
+Generator::run (float *data, const uint32_t n_samples)
+{
+	switch (_type) {
+		default:
+		case UniformWhiteNoise:
+			for (uint32_t i = 0; i < n_samples; ++i) {
+				data[i] = randf();
+			}
+			break;
+		case GaussianWhiteNoise:
+			for (uint32_t i = 0 ; i < n_samples; ++i) {
+				data[i] = 0.7079f * grandf();
+			}
+			break;
+		case PinkNoise:
+			for (uint32_t i = 0 ; i < n_samples; ++i) {
+				const float white = .39572f * randf ();
+				_b0 = .99886f * _b0 + white * .0555179f;
+				_b1 = .99332f * _b1 + white * .0750759f;
+				_b2 = .96900f * _b2 + white * .1538520f;
+				_b3 = .86650f * _b3 + white * .3104856f;
+				_b4 = .55000f * _b4 + white * .5329522f;
+				_b5 = -.7616f * _b5 - white * .0168980f;
+				data[i] = _b0 + _b1 + _b2 + _b3 + _b4 + _b5 + _b6 + white * 0.5362f;
+				_b6 = white * 0.115926f;
+			}
+			break;
+	}
+}
+
+inline uint32_t
+Generator::randi ()
+{
+	// 31bit Park-Miller-Carta Pseudo-Random Number Generator
+	uint32_t hi, lo;
+	lo = 16807 * (_rseed & 0xffff);
+	hi = 16807 * (_rseed >> 16);
+	lo += (hi & 0x7fff) << 16;
+	lo += hi >> 15;
+	lo = (lo & 0x7fffffff) + (lo >> 31);
+	return (_rseed = lo);
+}
+
+inline float
+Generator::grandf ()
+{
+	float x1, x2, r;
+
+	if (_pass) {
+		_pass = false;
+		return _rn;
+	}
+
+	do {
+		x1 = randf ();
+		x2 = randf ();
+		r = x1 * x1 + x2 * x2;
+	} while ((r >= 1.0f) || (r < 1e-22f));
+
+	r = sqrtf (-2.f * logf (r) / r);
+
+	_pass = true;
+	_rn = r * x2;
+	return r * x1;
+}