gst/replaygain/rganalysis.c

   1 /* GStreamer ReplayGain analysis
   2  *
   3  * Copyright (C) 2006 Rene Stadler <mail@renestadler.de>
   4  * Copyright (C) 2001 David Robinson <David@Robinson.org>
   5  *                    Glen Sawyer <glensawyer@hotmail.com>
   6  *
   7  * rganalysis.c: Analyze raw audio data in accordance with ReplayGain
   8  *
   9  * This library is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1 of
  12  * the License, or (at your option) any later version.
  13  *
  14  * This library is distributed in the hope that it will be useful, but
  15  * WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with this library; if not, write to the Free Software
  21  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  22  * 02110-1301 USA
  23  */
  24
  25 /* Based on code with Copyright (C) 2001 David Robinson
  26  * <David@Robinson.org> and Glen Sawyer <glensawyer@hotmail.com>,
  27  * which is distributed under the LGPL as part of the vorbisgain
  28  * program.  The original code also mentions Frank Klemm
  29  * (http://www.uni-jena.de/~pfk/mpp/) for having contributed lots of
  30  * good code.  Specifically, this is based on the file
  31  * "gain_analysis.c" from vorbisgain version 0.34.
  32  */
  33
  34 /* Room for future improvement: Mono data is currently in fact copied
  35  * to two channels which get processed normally.  This means that mono
  36  * input data is processed twice.
  37  */
  38
  39 /* Helpful information for understanding this code: The two IIR
  40  * filters depend on previous input _and_ previous output samples (up
  41  * to the filter's order number of samples).  This explains the whole
  42  * lot of memcpy'ing done in rg_analysis_analyze and why the context
  43  * holds so many buffers.
  44  */
  45
  46 #include <math.h>
  47 #include <string.h>
  48 #include <glib.h>
  49
  50 #include "rganalysis.h"
  51
  52 #define YULE_ORDER         10
  53 #define BUTTER_ORDER        2
  54 /* Percentile which is louder than the proposed level: */
  55 #define RMS_PERCENTILE     95
  56 /* Duration of RMS window in milliseconds: */
  57 #define RMS_WINDOW_MSECS   50
  58 /* Histogram array elements per dB: */
  59 #define STEPS_PER_DB      100
  60 /* Histogram upper bound in dB (normal max. values in the wild are
  61  * assumed to be around 70, 80 dB): */
  62 #define MAX_DB            120
  63 /* Calibration value: */
  64 #define PINK_REF           64.82        /* 298640883795 */
  65
  66 #define MAX_ORDER         MAX (BUTTER_ORDER, YULE_ORDER)
  67 #define MAX_SAMPLE_RATE   48000
  68 /* The + 999 has the effect of ceil()ing: */
  69 #define MAX_SAMPLE_WINDOW (guint) \
  70   ((MAX_SAMPLE_RATE * RMS_WINDOW_MSECS + 999) / 1000)
  71
  72 /* Analysis result accumulator. */
  73
  74 struct _RgAnalysisAcc
  75 {
  76   guint32 histogram[STEPS_PER_DB * MAX_DB];
  77   gdouble peak;
  78 };
  79
  80 typedef struct _RgAnalysisAcc RgAnalysisAcc;
  81
  82 /* Analysis context. */
  83
  84 struct _RgAnalysisCtx
  85 {
  86   /* Filter buffers for left channel. */
  87   gfloat inprebuf_l[MAX_ORDER * 2];
  88   gfloat *inpre_l;
  89   gfloat stepbuf_l[MAX_SAMPLE_WINDOW + MAX_ORDER];
  90   gfloat *step_l;
  91   gfloat outbuf_l[MAX_SAMPLE_WINDOW + MAX_ORDER];
  92   gfloat *out_l;
  93   /* Filter buffers for right channel. */
  94   gfloat inprebuf_r[MAX_ORDER * 2];
  95   gfloat *inpre_r;
  96   gfloat stepbuf_r[MAX_SAMPLE_WINDOW + MAX_ORDER];
  97   gfloat *step_r;
  98   gfloat outbuf_r[MAX_SAMPLE_WINDOW + MAX_ORDER];
  99   gfloat *out_r;
 100
 101   /* Number of samples to reach duration of the RMS window: */
 102   guint window_n_samples;
 103   /* Progress of the running window: */
 104   guint window_n_samples_done;
 105   gdouble window_square_sum;
 106
 107   gint sample_rate;
 108   gint sample_rate_index;
 109
 110   RgAnalysisAcc track;
 111   RgAnalysisAcc album;
 112   void (*post_message) (gpointer analysis,
 113       GstClockTime timestamp, GstClockTime duration, gdouble rglevel);
 114   gpointer analysis;
 115   /* The timestamp of the current incoming buffer. */
 116   GstClockTime buffer_timestamp;
 117   /* Number of samples processed in current buffer, during emit_signal,
 118      this will always be on an RMS window boundary. */
 119   guint buffer_n_samples_done;
 120 };
 121
 122 /* Filter coefficients for the IIR filters that form the equal
 123  * loudness filter.  XFilter[ctx->sample_rate_index] gives the array
 124  * of the X coefficients (A or B) for the configured sample rate. */
 125
 126 #ifdef _MSC_VER
 127 /* Disable double-to-float warning: */
 128 /* A better solution would be to append 'f' to each constant, but that
 129  * makes the code ugly. */
 130 #pragma warning ( disable : 4305 )
 131 #endif
 132
 133 static const gfloat AYule[9][11] = {
 134   {1., -3.84664617118067, 7.81501653005538, -11.34170355132042,
 135         13.05504219327545, -12.28759895145294, 9.48293806319790,
 136         -5.87257861775999, 2.75465861874613, -0.86984376593551,
 137       0.13919314567432},
 138   {1., -3.47845948550071, 6.36317777566148, -8.54751527471874, 9.47693607801280,
 139         -8.81498681370155, 6.85401540936998, -4.39470996079559,
 140       2.19611684890774, -0.75104302451432, 0.13149317958808},
 141   {1., -2.37898834973084, 2.84868151156327, -2.64577170229825, 2.23697657451713,
 142         -1.67148153367602, 1.00595954808547, -0.45953458054983,
 143       0.16378164858596, -0.05032077717131, 0.02347897407020},
 144   {1., -1.61273165137247, 1.07977492259970, -0.25656257754070,
 145         -0.16276719120440, -0.22638893773906, 0.39120800788284,
 146         -0.22138138954925, 0.04500235387352, 0.02005851806501,
 147       0.00302439095741},
 148   {1., -1.49858979367799, 0.87350271418188, 0.12205022308084, -0.80774944671438,
 149         0.47854794562326, -0.12453458140019, -0.04067510197014,
 150       0.08333755284107, -0.04237348025746, 0.02977207319925},
 151   {1., -0.62820619233671, 0.29661783706366, -0.37256372942400, 0.00213767857124,
 152         -0.42029820170918, 0.22199650564824, 0.00613424350682, 0.06747620744683,
 153       0.05784820375801, 0.03222754072173},
 154   {1., -1.04800335126349, 0.29156311971249, -0.26806001042947, 0.00819999645858,
 155         0.45054734505008, -0.33032403314006, 0.06739368333110,
 156       -0.04784254229033, 0.01639907836189, 0.01807364323573},
 157   {1., -0.51035327095184, -0.31863563325245, -0.20256413484477,
 158         0.14728154134330, 0.38952639978999, -0.23313271880868,
 159         -0.05246019024463, -0.02505961724053, 0.02442357316099,
 160       0.01818801111503},
 161   {1., -0.25049871956020, -0.43193942311114, -0.03424681017675,
 162         -0.04678328784242, 0.26408300200955, 0.15113130533216,
 163         -0.17556493366449, -0.18823009262115, 0.05477720428674,
 164       0.04704409688120}
 165 };
 166
 167 static const gfloat BYule[9][11] = {
 168   {0.03857599435200, -0.02160367184185, -0.00123395316851, -0.00009291677959,
 169         -0.01655260341619, 0.02161526843274, -0.02074045215285,
 170       0.00594298065125, 0.00306428023191, 0.00012025322027, 0.00288463683916},
 171   {0.05418656406430, -0.02911007808948, -0.00848709379851, -0.00851165645469,
 172         -0.00834990904936, 0.02245293253339, -0.02596338512915,
 173         0.01624864962975, -0.00240879051584, 0.00674613682247,
 174       -0.00187763777362},
 175   {0.15457299681924, -0.09331049056315, -0.06247880153653, 0.02163541888798,
 176         -0.05588393329856, 0.04781476674921, 0.00222312597743, 0.03174092540049,
 177       -0.01390589421898, 0.00651420667831, -0.00881362733839},
 178   {0.30296907319327, -0.22613988682123, -0.08587323730772, 0.03282930172664,
 179         -0.00915702933434, -0.02364141202522, -0.00584456039913,
 180         0.06276101321749, -0.00000828086748, 0.00205861885564,
 181       -0.02950134983287},
 182   {0.33642304856132, -0.25572241425570, -0.11828570177555, 0.11921148675203,
 183         -0.07834489609479, -0.00469977914380, -0.00589500224440,
 184         0.05724228140351, 0.00832043980773, -0.01635381384540,
 185       -0.01760176568150},
 186   {0.44915256608450, -0.14351757464547, -0.22784394429749, -0.01419140100551,
 187         0.04078262797139, -0.12398163381748, 0.04097565135648, 0.10478503600251,
 188       -0.01863887810927, -0.03193428438915, 0.00541907748707},
 189   {0.56619470757641, -0.75464456939302, 0.16242137742230, 0.16744243493672,
 190         -0.18901604199609, 0.30931782841830, -0.27562961986224,
 191         0.00647310677246, 0.08647503780351, -0.03788984554840,
 192       -0.00588215443421},
 193   {0.58100494960553, -0.53174909058578, -0.14289799034253, 0.17520704835522,
 194         0.02377945217615, 0.15558449135573, -0.25344790059353, 0.01628462406333,
 195       0.06920467763959, -0.03721611395801, -0.00749618797172},
 196   {0.53648789255105, -0.42163034350696, -0.00275953611929, 0.04267842219415,
 197         -0.10214864179676, 0.14590772289388, -0.02459864859345,
 198         -0.11202315195388, -0.04060034127000, 0.04788665548180,
 199       -0.02217936801134}
 200 };
 201
 202 static const gfloat AButter[9][3] = {
 203   {1., -1.97223372919527, 0.97261396931306},
 204   {1., -1.96977855582618, 0.97022847566350},
 205   {1., -1.95835380975398, 0.95920349965459},
 206   {1., -1.95002759149878, 0.95124613669835},
 207   {1., -1.94561023566527, 0.94705070426118},
 208   {1., -1.92783286977036, 0.93034775234268},
 209   {1., -1.91858953033784, 0.92177618768381},
 210   {1., -1.91542108074780, 0.91885558323625},
 211   {1., -1.88903307939452, 0.89487434461664}
 212 };
 213
 214 static const gfloat BButter[9][3] = {
 215   {0.98621192462708, -1.97242384925416, 0.98621192462708},
 216   {0.98500175787242, -1.97000351574484, 0.98500175787242},
 217   {0.97938932735214, -1.95877865470428, 0.97938932735214},
 218   {0.97531843204928, -1.95063686409857, 0.97531843204928},
 219   {0.97316523498161, -1.94633046996323, 0.97316523498161},
 220   {0.96454515552826, -1.92909031105652, 0.96454515552826},
 221   {0.96009142950541, -1.92018285901082, 0.96009142950541},
 222   {0.95856916599601, -1.91713833199203, 0.95856916599601},
 223   {0.94597685600279, -1.89195371200558, 0.94597685600279}
 224 };
 225
 226 #ifdef _MSC_VER
 227 #pragma warning ( default : 4305 )
 228 #endif
 229
 230 /* Filter functions.  These access elements with negative indices of
 231  * the input and output arrays (up to the filter's order). */
 232
 233 /* For much better performance, the function below has been
 234  * implemented by unrolling the inner loop for our two use cases. */
 235
 236 /*
 237  * static inline void
 238  * apply_filter (const gfloat * input, gfloat * output, guint n_samples,
 239  *     const gfloat * a, const gfloat * b, guint order)
 240  * {
 241  *   gfloat y;
 242  *   gint i, k;
 243  *
 244  *   for (i = 0; i < n_samples; i++) {
 245  *     y = input[i] * b[0];
 246  *     for (k = 1; k <= order; k++)
 247  *       y += input[i - k] * b[k] - output[i - k] * a[k];
 248  *     output[i] = y;
 249  *   }
 250  * }
 251  */
 252
 253 static inline void
 254 yule_filter (const gfloat * input, gfloat * output,
 255     const gfloat * a, const gfloat * b)
 256 {
 257   /* 1e-10 is added below to avoid running into denormals when operating on
 258    * near silence. */
 259
 260   output[0] = 1e-10 + input[0] * b[0]
 261       + input[-1] * b[1] - output[-1] * a[1]
 262       + input[-2] * b[2] - output[-2] * a[2]
 263       + input[-3] * b[3] - output[-3] * a[3]
 264       + input[-4] * b[4] - output[-4] * a[4]
 265       + input[-5] * b[5] - output[-5] * a[5]
 266       + input[-6] * b[6] - output[-6] * a[6]
 267       + input[-7] * b[7] - output[-7] * a[7]
 268       + input[-8] * b[8] - output[-8] * a[8]
 269       + input[-9] * b[9] - output[-9] * a[9]
 270       + input[-10] * b[10] - output[-10] * a[10];
 271 }
 272
 273 static inline void
 274 butter_filter (const gfloat * input, gfloat * output,
 275     const gfloat * a, const gfloat * b)
 276 {
 277   output[0] = input[0] * b[0]
 278       + input[-1] * b[1] - output[-1] * a[1]
 279       + input[-2] * b[2] - output[-2] * a[2];
 280 }
 281
 282 /* Because butter_filter and yule_filter are inlined, this function is
 283  * a bit blown-up (code-size wise), but not inlining gives a ca. 40%
 284  * performance penalty. */
 285
 286 static inline void
 287 apply_filters (const RgAnalysisCtx * ctx, const gfloat * input_l,
 288     const gfloat * input_r, guint n_samples)
 289 {
 290   const gfloat *ayule = AYule[ctx->sample_rate_index];
 291   const gfloat *byule = BYule[ctx->sample_rate_index];
 292   const gfloat *abutter = AButter[ctx->sample_rate_index];
 293   const gfloat *bbutter = BButter[ctx->sample_rate_index];
 294   gint pos = ctx->window_n_samples_done;
 295   gint i;
 296
 297   for (i = 0; i < n_samples; i++, pos++) {
 298     yule_filter (input_l + i, ctx->step_l + pos, ayule, byule);
 299     butter_filter (ctx->step_l + pos, ctx->out_l + pos, abutter, bbutter);
 300
 301     yule_filter (input_r + i, ctx->step_r + pos, ayule, byule);
 302     butter_filter (ctx->step_r + pos, ctx->out_r + pos, abutter, bbutter);
 303   }
 304 }
 305
 306 /* Clear filter buffer state and current RMS window. */
 307
 308 static void
 309 reset_filters (RgAnalysisCtx * ctx)
 310 {
 311   gint i;
 312
 313   for (i = 0; i < MAX_ORDER; i++) {
 314
 315     ctx->inprebuf_l[i] = 0.;
 316     ctx->stepbuf_l[i] = 0.;
 317     ctx->outbuf_l[i] = 0.;
 318
 319     ctx->inprebuf_r[i] = 0.;
 320     ctx->stepbuf_r[i] = 0.;
 321     ctx->outbuf_r[i] = 0.;
 322   }
 323
 324   ctx->window_square_sum = 0.;
 325   ctx->window_n_samples_done = 0;
 326 }
 327
 328 /* Accumulator functions. */
 329
 330 /* Add two accumulators in-place.  The sum is defined as the result of
 331  * the vector sum of the histogram array and the maximum value of the
 332  * peak field.  Thus "adding" the accumulators for all tracks yields
 333  * the correct result for obtaining the album gain and peak. */
 334
 335 static void
 336 accumulator_add (RgAnalysisAcc * acc, const RgAnalysisAcc * acc_other)
 337 {
 338   gint i;
 339
 340   for (i = 0; i < G_N_ELEMENTS (acc->histogram); i++)
 341     acc->histogram[i] += acc_other->histogram[i];
 342
 343   acc->peak = MAX (acc->peak, acc_other->peak);
 344 }
 345
 346 /* Reset an accumulator to zero. */
 347
 348 static void
 349 accumulator_clear (RgAnalysisAcc * acc)
 350 {
 351   memset (acc->histogram, 0, sizeof (acc->histogram));
 352   acc->peak = 0.;
 353 }
 354
 355 /* Obtain final analysis result from an accumulator.  Returns TRUE on
 356  * success, FALSE on error (if accumulator is still zero). */
 357
 358 static gboolean
 359 accumulator_result (const RgAnalysisAcc * acc, gdouble * result_gain,
 360     gdouble * result_peak)
 361 {
 362   guint32 sum = 0;
 363   guint32 upper;
 364   guint i;
 365
 366   for (i = 0; i < G_N_ELEMENTS (acc->histogram); i++)
 367     sum += acc->histogram[i];
 368
 369   if (sum == 0)
 370     /* All entries are 0: We got less than 50ms of data. */
 371     return FALSE;
 372
 373   upper = (guint32) ceil (sum * (1. - (gdouble) (RMS_PERCENTILE / 100.)));
 374
 375   for (i = G_N_ELEMENTS (acc->histogram); i--;) {
 376     if (upper <= acc->histogram[i])
 377       break;
 378     upper -= acc->histogram[i];
 379   }
 380
 381   if (result_peak != NULL)
 382     *result_peak = acc->peak;
 383   if (result_gain != NULL)
 384     *result_gain = PINK_REF - (gdouble) i / STEPS_PER_DB;
 385
 386   return TRUE;
 387 }
 388
 389 /* Functions that operate on contexts, for external usage. */
 390
 391 /* Create a new context.  Before it can be used, a sample rate must be
 392  * configured using rg_analysis_set_sample_rate. */
 393
 394 RgAnalysisCtx *
 395 rg_analysis_new (void)
 396 {
 397   RgAnalysisCtx *ctx;
 398
 399   ctx = g_new (RgAnalysisCtx, 1);
 400
 401   ctx->inpre_l = ctx->inprebuf_l + MAX_ORDER;
 402   ctx->step_l = ctx->stepbuf_l + MAX_ORDER;
 403   ctx->out_l = ctx->outbuf_l + MAX_ORDER;
 404
 405   ctx->inpre_r = ctx->inprebuf_r + MAX_ORDER;
 406   ctx->step_r = ctx->stepbuf_r + MAX_ORDER;
 407   ctx->out_r = ctx->outbuf_r + MAX_ORDER;
 408
 409   ctx->sample_rate = 0;
 410
 411   accumulator_clear (&ctx->track);
 412   accumulator_clear (&ctx->album);
 413
 414   return ctx;
 415 }
 416
 417 static void
 418 reset_silence_detection (RgAnalysisCtx * ctx)
 419 {
 420   ctx->buffer_timestamp = GST_CLOCK_TIME_NONE;
 421   ctx->buffer_n_samples_done = 0;
 422 }
 423
 424 /* Adapt to given sample rate.  Does nothing if already the current
 425  * rate (returns TRUE then).  Returns FALSE only if given sample rate
 426  * is not supported.  If the configured rate changes, the last
 427  * unprocessed incomplete 50ms chunk of data is dropped because the
 428  * filters are reset. */
 429
 430 gboolean
 431 rg_analysis_set_sample_rate (RgAnalysisCtx * ctx, gint sample_rate)
 432 {
 433   g_return_val_if_fail (ctx != NULL, FALSE);
 434
 435   if (ctx->sample_rate == sample_rate)
 436     return TRUE;
 437
 438   switch (sample_rate) {
 439     case 48000:
 440       ctx->sample_rate_index = 0;
 441       break;
 442     case 44100:
 443       ctx->sample_rate_index = 1;
 444       break;
 445     case 32000:
 446       ctx->sample_rate_index = 2;
 447       break;
 448     case 24000:
 449       ctx->sample_rate_index = 3;
 450       break;
 451     case 22050:
 452       ctx->sample_rate_index = 4;
 453       break;
 454     case 16000:
 455       ctx->sample_rate_index = 5;
 456       break;
 457     case 12000:
 458       ctx->sample_rate_index = 6;
 459       break;
 460     case 11025:
 461       ctx->sample_rate_index = 7;
 462       break;
 463     case 8000:
 464       ctx->sample_rate_index = 8;
 465       break;
 466     default:
 467       return FALSE;
 468   }
 469
 470   ctx->sample_rate = sample_rate;
 471   /* The + 999 has the effect of ceil()ing: */
 472   ctx->window_n_samples = (guint) ((sample_rate * RMS_WINDOW_MSECS + 999)
 473       / 1000);
 474
 475   reset_filters (ctx);
 476   reset_silence_detection (ctx);
 477
 478   return TRUE;
 479 }
 480
 481 void
 482 rg_analysis_init_silence_detection (RgAnalysisCtx * ctx,
 483     void (*post_message) (gpointer analysis, GstClockTime timestamp,
 484         GstClockTime duration, gdouble rglevel), gpointer analysis)
 485 {
 486   ctx->post_message = post_message;
 487   ctx->analysis = analysis;
 488   reset_silence_detection (ctx);
 489 }
 490
 491 void
 492 rg_analysis_start_buffer (RgAnalysisCtx * ctx, GstClockTime buffer_timestamp)
 493 {
 494   ctx->buffer_timestamp = buffer_timestamp;
 495   ctx->buffer_n_samples_done = 0;
 496 }
 497
 498 void
 499 rg_analysis_destroy (RgAnalysisCtx * ctx)
 500 {
 501   g_free (ctx);
 502 }
 503
 504 /* Entry points for analyzing sample data in common raw data formats.
 505  * The stereo format functions expect interleaved frames.  It is
 506  * possible to pass data in different formats for the same context,
 507  * there are no restrictions.  All functions have the same signature;
 508  * the depth argument for the float functions is not variable and must
 509  * be given the value 32. */
 510
 511 void
 512 rg_analysis_analyze_mono_float (RgAnalysisCtx * ctx, gconstpointer data,
 513     gsize size, guint depth)
 514 {
 515   gfloat conv_samples[512];
 516   const gfloat *samples = (gfloat *) data;
 517   guint n_samples = size / sizeof (gfloat);
 518   gint i;
 519
 520   g_return_if_fail (depth == 32);
 521   g_return_if_fail (size % sizeof (gfloat) == 0);
 522
 523   while (n_samples) {
 524     gint n = MIN (n_samples, G_N_ELEMENTS (conv_samples));
 525
 526     n_samples -= n;
 527     memcpy (conv_samples, samples, n * sizeof (gfloat));
 528     for (i = 0; i < n; i++) {
 529       ctx->track.peak = MAX (ctx->track.peak, fabs (conv_samples[i]));
 530       conv_samples[i] *= 32768.;
 531     }
 532     samples += n;
 533     rg_analysis_analyze (ctx, conv_samples, NULL, n);
 534   }
 535 }
 536
 537 void
 538 rg_analysis_analyze_stereo_float (RgAnalysisCtx * ctx, gconstpointer data,
 539     gsize size, guint depth)
 540 {
 541   gfloat conv_samples_l[256];
 542   gfloat conv_samples_r[256];
 543   const gfloat *samples = (gfloat *) data;
 544   guint n_frames = size / (sizeof (gfloat) * 2);
 545   gint i;
 546
 547   g_return_if_fail (depth == 32);
 548   g_return_if_fail (size % (sizeof (gfloat) * 2) == 0);
 549
 550   while (n_frames) {
 551     gint n = MIN (n_frames, G_N_ELEMENTS (conv_samples_l));
 552
 553     n_frames -= n;
 554     for (i = 0; i < n; i++) {
 555       gfloat old_sample;
 556
 557       old_sample = samples[2 * i];
 558       ctx->track.peak = MAX (ctx->track.peak, fabs (old_sample));
 559       conv_samples_l[i] = old_sample * 32768.;
 560
 561       old_sample = samples[2 * i + 1];
 562       ctx->track.peak = MAX (ctx->track.peak, fabs (old_sample));
 563       conv_samples_r[i] = old_sample * 32768.;
 564     }
 565     samples += 2 * n;
 566     rg_analysis_analyze (ctx, conv_samples_l, conv_samples_r, n);
 567   }
 568 }
 569
 570 void
 571 rg_analysis_analyze_mono_int16 (RgAnalysisCtx * ctx, gconstpointer data,
 572     gsize size, guint depth)
 573 {
 574   gfloat conv_samples[512];
 575   gint32 peak_sample = 0;
 576   const gint16 *samples = (gint16 *) data;
 577   guint n_samples = size / sizeof (gint16);
 578   gint shift = 1 << (sizeof (gint16) * 8 - depth);
 579   gint i;
 580
 581   g_return_if_fail (depth <= (sizeof (gint16) * 8));
 582   g_return_if_fail (size % sizeof (gint16) == 0);
 583
 584   while (n_samples) {
 585     gint n = MIN (n_samples, G_N_ELEMENTS (conv_samples));
 586
 587     n_samples -= n;
 588     for (i = 0; i < n; i++) {
 589       gint16 old_sample = samples[i] * shift;
 590
 591       peak_sample = MAX (peak_sample, ABS ((gint32) old_sample));
 592       conv_samples[i] = (gfloat) old_sample;
 593     }
 594     samples += n;
 595     rg_analysis_analyze (ctx, conv_samples, NULL, n);
 596   }
 597   ctx->track.peak = MAX (ctx->track.peak,
 598       (gdouble) peak_sample / ((gdouble) (1u << 15)));
 599 }
 600
 601 void
 602 rg_analysis_analyze_stereo_int16 (RgAnalysisCtx * ctx, gconstpointer data,
 603     gsize size, guint depth)
 604 {
 605   gfloat conv_samples_l[256];
 606   gfloat conv_samples_r[256];
 607   gint32 peak_sample = 0;
 608   const gint16 *samples = (gint16 *) data;
 609   guint n_frames = size / (sizeof (gint16) * 2);
 610   gint shift = 1 << (sizeof (gint16) * 8 - depth);
 611   gint i;
 612
 613   g_return_if_fail (depth <= (sizeof (gint16) * 8));
 614   g_return_if_fail (size % (sizeof (gint16) * 2) == 0);
 615
 616   while (n_frames) {
 617     gint n = MIN (n_frames, G_N_ELEMENTS (conv_samples_l));
 618
 619     n_frames -= n;
 620     for (i = 0; i < n; i++) {
 621       gint16 old_sample;
 622
 623       old_sample = samples[2 * i] * shift;
 624       peak_sample = MAX (peak_sample, ABS ((gint32) old_sample));
 625       conv_samples_l[i] = (gfloat) old_sample;
 626
 627       old_sample = samples[2 * i + 1] * shift;
 628       peak_sample = MAX (peak_sample, ABS ((gint32) old_sample));
 629       conv_samples_r[i] = (gfloat) old_sample;
 630     }
 631     samples += 2 * n;
 632     rg_analysis_analyze (ctx, conv_samples_l, conv_samples_r, n);
 633   }
 634   ctx->track.peak = MAX (ctx->track.peak,
 635       (gdouble) peak_sample / ((gdouble) (1u << 15)));
 636 }
 637
 638 /* Analyze the given chunk of samples.  The sample data is given in
 639  * floating point format but should be scaled such that the values
 640  * +/-32768.0 correspond to the -0dBFS reference amplitude.
 641  *
 642  * samples_l: Buffer with sample data for the left channel or of the
 643  * mono channel.
 644  *
 645  * samples_r: Buffer with sample data for the right channel or NULL
 646  * for mono.
 647  *
 648  * n_samples: Number of samples passed in each buffer.
 649  */
 650
 651 void
 652 rg_analysis_analyze (RgAnalysisCtx * ctx, const gfloat * samples_l,
 653     const gfloat * samples_r, guint n_samples)
 654 {
 655   const gfloat *input_l, *input_r;
 656   guint n_samples_done;
 657   gint i;
 658
 659   g_return_if_fail (ctx != NULL);
 660   g_return_if_fail (samples_l != NULL);
 661   g_return_if_fail (ctx->sample_rate != 0);
 662
 663   if (n_samples == 0)
 664     return;
 665
 666   if (samples_r == NULL)
 667     /* Mono. */
 668     samples_r = samples_l;
 669
 670   memcpy (ctx->inpre_l, samples_l,
 671       MIN (n_samples, MAX_ORDER) * sizeof (gfloat));
 672   memcpy (ctx->inpre_r, samples_r,
 673       MIN (n_samples, MAX_ORDER) * sizeof (gfloat));
 674
 675   n_samples_done = 0;
 676   while (n_samples_done < n_samples) {
 677     /* Limit number of samples to be processed in this iteration to
 678      * the number needed to complete the next window: */
 679     guint n_samples_current = MIN (n_samples - n_samples_done,
 680         ctx->window_n_samples - ctx->window_n_samples_done);
 681
 682     if (n_samples_done < MAX_ORDER) {
 683       input_l = ctx->inpre_l + n_samples_done;
 684       input_r = ctx->inpre_r + n_samples_done;
 685       n_samples_current = MIN (n_samples_current, MAX_ORDER - n_samples_done);
 686     } else {
 687       input_l = samples_l + n_samples_done;
 688       input_r = samples_r + n_samples_done;
 689     }
 690
 691     apply_filters (ctx, input_l, input_r, n_samples_current);
 692
 693     /* Update the square sum. */
 694     for (i = 0; i < n_samples_current; i++)
 695       ctx->window_square_sum += ctx->out_l[ctx->window_n_samples_done + i]
 696           * ctx->out_l[ctx->window_n_samples_done + i]
 697           + ctx->out_r[ctx->window_n_samples_done + i]
 698           * ctx->out_r[ctx->window_n_samples_done + i];
 699
 700     ctx->window_n_samples_done += n_samples_current;
 701     ctx->buffer_n_samples_done += n_samples_current;
 702
 703     g_return_if_fail (ctx->window_n_samples_done <= ctx->window_n_samples);
 704
 705     if (ctx->window_n_samples_done == ctx->window_n_samples) {
 706       /* Get the Root Mean Square (RMS) for this set of samples. */
 707       gdouble val = STEPS_PER_DB * 10. * log10 (ctx->window_square_sum /
 708           ctx->window_n_samples * 0.5 + 1.e-37);
 709       gint ival = CLAMP ((gint) val, 0,
 710           (gint) G_N_ELEMENTS (ctx->track.histogram) - 1);
 711       /* Compute the per-window gain */
 712       const gdouble gain = PINK_REF - (gdouble) ival / STEPS_PER_DB;
 713       const GstClockTime timestamp = ctx->buffer_timestamp
 714           + gst_util_uint64_scale_int_ceil (GST_SECOND,
 715           ctx->buffer_n_samples_done,
 716           ctx->sample_rate)
 717           - RMS_WINDOW_MSECS * GST_MSECOND;
 718
 719       ctx->post_message (ctx->analysis, timestamp,
 720           RMS_WINDOW_MSECS * GST_MSECOND, -gain);
 721
 722
 723       ctx->track.histogram[ival]++;
 724       ctx->window_square_sum = 0.;
 725       ctx->window_n_samples_done = 0;
 726
 727       /* No need for memmove here, the areas never overlap: Even for
 728        * the smallest sample rate, the number of samples needed for
 729        * the window is greater than MAX_ORDER. */
 730
 731       memcpy (ctx->stepbuf_l, ctx->stepbuf_l + ctx->window_n_samples,
 732           MAX_ORDER * sizeof (gfloat));
 733       memcpy (ctx->outbuf_l, ctx->outbuf_l + ctx->window_n_samples,
 734           MAX_ORDER * sizeof (gfloat));
 735
 736       memcpy (ctx->stepbuf_r, ctx->stepbuf_r + ctx->window_n_samples,
 737           MAX_ORDER * sizeof (gfloat));
 738       memcpy (ctx->outbuf_r, ctx->outbuf_r + ctx->window_n_samples,
 739           MAX_ORDER * sizeof (gfloat));
 740     }
 741
 742     n_samples_done += n_samples_current;
 743   }
 744
 745   if (n_samples >= MAX_ORDER) {
 746
 747     memcpy (ctx->inprebuf_l, samples_l + n_samples - MAX_ORDER,
 748         MAX_ORDER * sizeof (gfloat));
 749
 750     memcpy (ctx->inprebuf_r, samples_r + n_samples - MAX_ORDER,
 751         MAX_ORDER * sizeof (gfloat));
 752
 753   } else {
 754
 755     memmove (ctx->inprebuf_l, ctx->inprebuf_l + n_samples,
 756         (MAX_ORDER - n_samples) * sizeof (gfloat));
 757     memcpy (ctx->inprebuf_l + MAX_ORDER - n_samples, samples_l,
 758         n_samples * sizeof (gfloat));
 759
 760     memmove (ctx->inprebuf_r, ctx->inprebuf_r + n_samples,
 761         (MAX_ORDER - n_samples) * sizeof (gfloat));
 762     memcpy (ctx->inprebuf_r + MAX_ORDER - n_samples, samples_r,
 763         n_samples * sizeof (gfloat));
 764
 765   }
 766 }
 767
 768 /* Obtain track gain and peak.  Returns TRUE on success.  Can fail if
 769  * not enough samples have been processed.  Updates album accumulator.
 770  * Resets track accumulator. */
 771
 772 gboolean
 773 rg_analysis_track_result (RgAnalysisCtx * ctx, gdouble * gain, gdouble * peak)
 774 {
 775   gboolean result;
 776
 777   g_return_val_if_fail (ctx != NULL, FALSE);
 778
 779   accumulator_add (&ctx->album, &ctx->track);
 780   result = accumulator_result (&ctx->track, gain, peak);
 781   accumulator_clear (&ctx->track);
 782
 783   reset_filters (ctx);
 784   reset_silence_detection (ctx);
 785
 786   return result;
 787 }
 788
 789 /* Obtain album gain and peak.  Returns TRUE on success.  Can fail if
 790  * not enough samples have been processed.  Resets album
 791  * accumulator. */
 792
 793 gboolean
 794 rg_analysis_album_result (RgAnalysisCtx * ctx, gdouble * gain, gdouble * peak)
 795 {
 796   gboolean result;
 797
 798   g_return_val_if_fail (ctx != NULL, FALSE);
 799
 800   result = accumulator_result (&ctx->album, gain, peak);
 801   accumulator_clear (&ctx->album);
 802
 803   return result;
 804 }
 805
 806 void
 807 rg_analysis_reset_album (RgAnalysisCtx * ctx)
 808 {
 809   accumulator_clear (&ctx->album);
 810 }
 811
 812 /* Reset internal buffers as well as track and album accumulators.
 813  * Configured sample rate is kept intact. */
 814
 815 void
 816 rg_analysis_reset (RgAnalysisCtx * ctx)
 817 {
 818   g_return_if_fail (ctx != NULL);
 819
 820   reset_filters (ctx);
 821   accumulator_clear (&ctx->track);
 822   accumulator_clear (&ctx->album);
 823   reset_silence_detection (ctx);
 824 }