libs/audiographer/src/general/analyser.cc

   1 /*
   2  * Copyright (C) 2016 Robin Gareus <robin@gareus.org>
   3  *
   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.
   8  *
   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.
  13  *
  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.
  17  */
  18
  19 #include "audiographer/general/analyser.h"
  20 #include "pbd/fastlog.h"
  21 #include "pbd/compose.h"
  22 #include "ardour/debug.h"
  23
  24 using namespace AudioGrapher;
  25
  26 const float Analyser::fft_range_db (120); // dB
  27
  28 Analyser::Analyser (float sample_rate, unsigned int channels, framecnt_t bufsize, framecnt_t n_samples)
  29         : _ebur128_plugin (0)
  30         , _dbtp_plugin (0)
  31         , _sample_rate (sample_rate)
  32         , _channels (channels)
  33         , _bufsize (bufsize / channels)
  34         , _n_samples (n_samples)
  35         , _pos (0)
  36 {
  37         DEBUG_TRACE (PBD::DEBUG::ExportAnalysis, string_compose ("Analyser r:%1 c:%2 f:%3 d:%4\n", sample_rate, channels, bufsize, n_samples));
  38
  39         assert (bufsize % channels == 0);
  40         assert (bufsize > 1);
  41
  42         if (channels > 0 && channels <= 2) {
  43                 using namespace Vamp::HostExt;
  44                 PluginLoader* loader (PluginLoader::getInstance ());
  45                 _ebur128_plugin = loader->loadPlugin ("libardourvampplugins:ebur128", sample_rate, PluginLoader::ADAPT_ALL_SAFE);
  46                 assert (_ebur128_plugin);
  47                 _ebur128_plugin->reset ();
  48                 if (!_ebur128_plugin->initialise (channels, _bufsize, _bufsize)) {
  49                         delete _ebur128_plugin;
  50                         _ebur128_plugin = 0;
  51                 }
  52         }
  53
  54         _dbtp_plugin = (Vamp::Plugin**) malloc (sizeof(Vamp::Plugin*) * channels);
  55         for (unsigned int c = 0; c < _channels; ++c) {
  56                 using namespace Vamp::HostExt;
  57                 PluginLoader* loader (PluginLoader::getInstance ());
  58                 _dbtp_plugin[c] = loader->loadPlugin ("libardourvampplugins:dBTP", sample_rate, PluginLoader::ADAPT_ALL_SAFE);
  59                 assert (_dbtp_plugin[c]);
  60                 _dbtp_plugin[c]->reset ();
  61                 if (!_dbtp_plugin[c]->initialise (1, _bufsize, _bufsize)) {
  62                         delete _dbtp_plugin[c];
  63                         _dbtp_plugin[c] = 0;
  64                 }
  65         }
  66
  67         _bufs[0] = (float*) malloc (sizeof (float) * _bufsize);
  68         _bufs[1] = (float*) malloc (sizeof (float) * _bufsize);
  69
  70         const size_t peaks = sizeof (_result.peaks) / sizeof (ARDOUR::PeakData::PeakDatum) / 4;
  71         _spp = ceil ((_n_samples + 2.f) / (float) peaks);
  72
  73         const size_t swh = sizeof (_result.spectrum) / sizeof (float);
  74         const size_t height = sizeof (_result.spectrum[0]) / sizeof (float);
  75         const size_t width = swh / height;
  76         _fpp = ceil ((_n_samples + 2.f) / (float) width);
  77
  78         _fft_data_size   = _bufsize / 2;
  79         _fft_freq_per_bin = sample_rate / _fft_data_size / 2.f;
  80
  81         _fft_data_in  = (float *) fftwf_malloc (sizeof (float) * _bufsize);
  82         _fft_data_out = (float *) fftwf_malloc (sizeof (float) * _bufsize);
  83         _fft_power    = (float *) malloc (sizeof (float) * _fft_data_size);
  84
  85         for (uint32_t i = 0; i < _fft_data_size; ++i) {
  86                 _fft_power[i] = 0;
  87         }
  88         for (uint32_t i = 0; i < _bufsize; ++i) {
  89                 _fft_data_out[i] = 0;
  90         }
  91
  92         const float nyquist = (sample_rate * .5);
  93 #if 0 // linear
  94 #define YPOS(FREQ) rint (height * (1.0 - FREQ / nyquist))
  95 #else
  96 #define YPOS(FREQ) rint (height * (1 - logf (1.f + .1f * _fft_data_size * FREQ / nyquist) / logf (1.f + .1f * _fft_data_size)))
  97 #endif
  98
  99         _result.freq[0] = YPOS (50);
 100         _result.freq[1] = YPOS (100);
 101         _result.freq[2] = YPOS (500);
 102         _result.freq[3] = YPOS (1000);
 103         _result.freq[4] = YPOS (5000);
 104         _result.freq[5] = YPOS (10000);
 105
 106         _fft_plan = fftwf_plan_r2r_1d (_bufsize, _fft_data_in, _fft_data_out, FFTW_R2HC, FFTW_MEASURE);
 107
 108         _hann_window = (float *) malloc (sizeof (float) * _bufsize);
 109         double sum = 0.0;
 110
 111         for (uint32_t i = 0; i < _bufsize; ++i) {
 112                 _hann_window[i] = 0.5f - (0.5f * (float) cos (2.0f * M_PI * (float)i / (float)(_bufsize)));
 113                 sum += _hann_window[i];
 114         }
 115         const double isum = 2.0 / sum;
 116         for (uint32_t i = 0; i < _bufsize; ++i) {
 117                 _hann_window[i] *= isum;
 118         }
 119
 120         if (channels == 2) {
 121                 _result.n_channels = 2;
 122         } else {
 123                 _result.n_channels = 1;
 124         }
 125 }
 126
 127 Analyser::~Analyser ()
 128 {
 129         delete _ebur128_plugin;
 130         for (unsigned int c = 0; c < _channels; ++c) {
 131                 delete _dbtp_plugin[c];
 132         }
 133         free (_dbtp_plugin);
 134         free (_bufs[0]);
 135         free (_bufs[1]);
 136         fftwf_destroy_plan (_fft_plan);
 137         fftwf_free (_fft_data_in);
 138         fftwf_free (_fft_data_out);
 139         free (_fft_power);
 140         free (_hann_window);
 141 }
 142
 143 void
 144 Analyser::process (ProcessContext<float> const & ctx)
 145 {
 146         const framecnt_t n_samples = ctx.frames () / ctx.channels ();
 147         assert (ctx.channels () == _channels);
 148         assert (ctx.frames () % ctx.channels () == 0);
 149         assert (n_samples <= _bufsize);
 150         //printf ("PROC %p @%ld F: %ld, S: %ld C:%d\n", this, _pos, ctx.frames (), n_samples, ctx.channels ());
 151
 152         // allow 1 sample slack for resampling
 153         if (_pos + n_samples > _n_samples + 1) {
 154                 _pos += n_samples;
 155                 ListedSource<float>::output (ctx);
 156                 return;
 157         }
 158
 159         float const * d = ctx.data ();
 160         framecnt_t s;
 161         const unsigned cmask = _result.n_channels - 1; // [0, 1]
 162         for (s = 0; s < n_samples; ++s) {
 163                 _fft_data_in[s] = 0;
 164                 const framecnt_t pbin = (_pos + s) / _spp;
 165                 for (unsigned int c = 0; c < _channels; ++c) {
 166                         const float v = *d;
 167                         if (fabsf(v) > _result.peak) { _result.peak = fabsf(v); }
 168                         _bufs[c][s] = v;
 169                         const unsigned int cc = c & cmask;
 170                         if (_result.peaks[cc][pbin].min > v) { _result.peaks[cc][pbin].min = *d; }
 171                         if (_result.peaks[cc][pbin].max < v) { _result.peaks[cc][pbin].max = *d; }
 172                         _fft_data_in[s] += v * _hann_window[s] / (float) _channels;
 173                         ++d;
 174                 }
 175         }
 176
 177         for (; s < _bufsize; ++s) {
 178                 _fft_data_in[s] = 0;
 179                 for (unsigned int c = 0; c < _channels; ++c) {
 180                         _bufs[c][s] = 0.f;
 181                 }
 182         }
 183
 184         if (_ebur128_plugin) {
 185                 _ebur128_plugin->process (_bufs, Vamp::RealTime::fromSeconds ((double) _pos / _sample_rate));
 186         }
 187
 188         float const * const data = ctx.data ();
 189         for (unsigned int c = 0; c < _channels; ++c) {
 190                 if (!_dbtp_plugin[c]) { continue; }
 191                 for (s = 0; s < n_samples; ++s) {
 192                         _bufs[0][s] = data[s * _channels + c];
 193                 }
 194                 _dbtp_plugin[c]->process (_bufs, Vamp::RealTime::fromSeconds ((double) _pos / _sample_rate));
 195         }
 196
 197         fftwf_execute (_fft_plan);
 198
 199         _fft_power[0] = _fft_data_out[0] * _fft_data_out[0];
 200 #define FRe (_fft_data_out[i])
 201 #define FIm (_fft_data_out[_bufsize - i])
 202         for (uint32_t i = 1; i < _fft_data_size - 1; ++i) {
 203                 _fft_power[i] = (FRe * FRe) + (FIm * FIm);
 204         }
 205 #undef FRe
 206 #undef FIm
 207
 208         const size_t height = sizeof (_result.spectrum[0]) / sizeof (float);
 209         const framecnt_t x0 = _pos / _fpp;
 210         framecnt_t x1 = (_pos + n_samples) / _fpp;
 211         if (x0 == x1) x1 = x0 + 1;
 212
 213         for (uint32_t i = 0; i < _fft_data_size - 1; ++i) {
 214                 const float level = fft_power_at_bin (i, i);
 215                 if (level < -fft_range_db) continue;
 216                 const float pk = level > 0.0 ? 1.0 : (fft_range_db + level) / fft_range_db;
 217 #if 0 // linear
 218                 const uint32_t y0 = floor (i * (float) height / _fft_data_size);
 219                 uint32_t y1 = ceil ((i + 1.0) * (float) height / _fft_data_size);
 220 #else // logscale
 221                 const uint32_t y0 = floor (height * logf (1.f + .1f * i) / logf (1.f + .1f * _fft_data_size));
 222                 uint32_t y1 = ceilf (height * logf (1.f + .1f * (i + 1.f)) / logf (1.f + .1f * _fft_data_size));
 223 #endif
 224                 assert (y0 < height);
 225                 assert (y1 > 0 && y1 <= height);
 226                 if (y0 == y1) y1 = y0 + 1;
 227                 for (int x = x0; x < x1; ++x) {
 228                         for (uint32_t y = y0; y < y1 && y < height; ++y) {
 229                                 uint32_t yy = height - 1 - y;
 230                                 if (_result.spectrum[x][yy] < pk) { _result.spectrum[x][yy] = pk; }
 231                         }
 232                 }
 233         }
 234
 235         _pos += n_samples;
 236
 237         /* pass audio audio through */
 238         ListedSource<float>::output (ctx);
 239 }
 240
 241 ARDOUR::ExportAnalysisPtr
 242 Analyser::result ()
 243 {
 244         DEBUG_TRACE (PBD::DEBUG::ExportAnalysis, string_compose ("Processed %1 / %2 samples\n", _pos, _n_samples));
 245         if (_pos == 0 || ::llabs (_pos -_n_samples) > 1) {
 246                 printf("NO ANAL\n");
 247                 return ARDOUR::ExportAnalysisPtr ();
 248         }
 249
 250         if (_ebur128_plugin) {
 251                 Vamp::Plugin::FeatureSet features = _ebur128_plugin->getRemainingFeatures ();
 252                 if (!features.empty () && features.size () == 3) {
 253                         _result.loudness = features[0][0].values[0];
 254                         _result.loudness_range = features[1][0].values[0];
 255                         assert (features[2][0].values.size () == 540);
 256                         for (int i = 0; i < 540; ++i) {
 257                                 _result.loudness_hist[i] = features[2][0].values[i];
 258                                 if (_result.loudness_hist[i] > _result.loudness_hist_max) {
 259                                         _result.loudness_hist_max = _result.loudness_hist[i]; }
 260                         }
 261                         _result.have_loudness = true;
 262                 }
 263         }
 264
 265         const unsigned cmask = _result.n_channels - 1; // [0, 1]
 266         for (unsigned int c = 0; c < _channels; ++c) {
 267                 if (!_dbtp_plugin[c]) { continue; }
 268                 Vamp::Plugin::FeatureSet features = _dbtp_plugin[c]->getRemainingFeatures ();
 269                 if (!features.empty () && features.size () == 2) {
 270                         _result.have_dbtp = true;
 271                         float p = features[0][0].values[0];
 272                         if (p > _result.truepeak) { _result.truepeak = p; }
 273
 274                         for (std::vector<float>::const_iterator i = features[1][0].values.begin();
 275                                         i != features[1][0].values.end(); ++i) {
 276                                 const framecnt_t pk = (*i) / _spp;
 277                                 const unsigned int cc = c & cmask;
 278                                 _result.truepeakpos[cc].insert (pk);
 279                         }
 280                 }
 281         }
 282
 283         return ARDOUR::ExportAnalysisPtr (new ARDOUR::ExportAnalysis (_result));
 284 }
 285
 286 float
 287 Analyser::fft_power_at_bin (const uint32_t b, const float norm) const
 288 {
 289         const float a = _fft_power[b] * norm;
 290         return a > 1e-12 ? 10.0 * fast_log10 (a) : -INFINITY;
 291 }