+ 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)
+{
+ 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;