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