1 /* replaygain_synthesis - Routines for applying ReplayGain to a signal
2 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009 Josh Coalson
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2.1 of the License, or (at your option) any later version.
9 * This library 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 GNU
12 * Lesser General Public License for more details.
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with this library; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 * This is an aggregation of pieces of code from John Edwards' WaveGain
20 * program. Mostly cosmetic changes were made; otherwise, the dithering
21 * code is almost untouched and the gain processing was converted from
22 * processing a whole file to processing chunks of samples.
24 * The original copyright notices for WaveGain's dither.c and wavegain.c
28 * (c) 2002 John Edwards
29 * mostly lifted from work by Frank Klemm
30 * random functions for dithering.
33 * Copyright (C) 2002 John Edwards
34 * Additional code by Magnus Holmgren and Gian-Carlo Pascutto
41 #include <string.h> /* for memset() */
43 #include "private/fast_float_math_hack.h"
44 #include "replaygain_synthesis.h"
45 #include "FLAC/assert.h"
47 /* adjust for compilers that can't understand using LL suffix for int64_t literals */
49 #define FLAC__I64L(x) x
51 #define FLAC__I64L(x) x##LL
56 * the following is based on parts of dither.c
61 * This is a simple random number generator with good quality for audio purposes.
62 * It consists of two polycounters with opposite rotation direction and different
63 * periods. The periods are coprime, so the total period is the product of both.
65 * -------------------------------------------------------------------------------------------------
66 * +-> |31:30:29:28:27:26:25:24:23:22:21:20:19:18:17:16:15:14:13:12:11:10: 9: 8: 7: 6: 5: 4: 3: 2: 1: 0|
67 * | -------------------------------------------------------------------------------------------------
69 * | +--+--+--+-XOR-+--------+
71 * +--------------------------------------------------------------------------------------+
73 * -------------------------------------------------------------------------------------------------
74 * |31:30:29:28:27:26:25:24:23:22:21:20:19:18:17:16:15:14:13:12:11:10: 9: 8: 7: 6: 5: 4: 3: 2: 1: 0| <-+
75 * ------------------------------------------------------------------------------------------------- |
77 * +--+----XOR----+--+ |
79 * +----------------------------------------------------------------------------------------+
82 * The first has an period of 3*5*17*257*65537, the second of 7*47*73*178481,
83 * which gives a period of 18.410.713.077.675.721.215. The result is the
84 * XORed values of both generators.
87 static unsigned int random_int_(void)
89 static const unsigned char parity_[256] = {
90 0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
91 1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,
92 1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,
93 0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
94 1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,
95 0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
96 0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
97 1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0
99 static unsigned int r1_ = 1;
100 static unsigned int r2_ = 1;
102 unsigned int t1, t2, t3, t4;
104 /* Parity calculation is done via table lookup, this is also available
105 * on CPUs without parity, can be implemented in C and avoid unpredictable
106 * jumps and slow rotate through the carry flag operations.
108 t3 = t1 = r1_; t4 = t2 = r2_;
109 t1 &= 0xF5; t2 >>= 25;
110 t1 = parity_[t1]; t2 &= 0x63;
111 t1 <<= 31; t2 = parity_[t2];
113 return (r1_ = (t3 >> 1) | t1 ) ^ (r2_ = (t4 + t4) | t2 );
116 /* gives a equal distributed random number */
117 /* between -2^31*mult and +2^31*mult */
118 static double random_equi_(double mult)
120 return mult * (int) random_int_();
123 /* gives a triangular distributed random number */
124 /* between -2^32*mult and +2^32*mult */
125 static double random_triangular_(double mult)
127 return mult * ( (double) (int) random_int_() + (double) (int) random_int_() );
131 static const float F44_0 [16 + 32] = {
132 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
133 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
135 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
136 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
138 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
139 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0
143 static const float F44_1 [16 + 32] = { /* SNR(w) = 4.843163 dB, SNR = -3.192134 dB */
144 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
145 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
146 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
147 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,
149 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
150 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
151 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
152 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,
154 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
155 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
156 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
157 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,
161 static const float F44_2 [16 + 32] = { /* SNR(w) = 10.060213 dB, SNR = -12.766730 dB */
162 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
163 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
164 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
165 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,
167 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
168 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
169 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
170 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,
172 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
173 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
174 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
175 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,
179 static const float F44_3 [16 + 32] = { /* SNR(w) = 15.382598 dB, SNR = -29.402334 dB */
180 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
181 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
182 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
183 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099,
185 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
186 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
187 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
188 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099,
190 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
191 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
192 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
193 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099
197 static double scalar16_(const float* x, const float* y)
200 x[ 0]*y[ 0] + x[ 1]*y[ 1] + x[ 2]*y[ 2] + x[ 3]*y[ 3] +
201 x[ 4]*y[ 4] + x[ 5]*y[ 5] + x[ 6]*y[ 6] + x[ 7]*y[ 7] +
202 x[ 8]*y[ 8] + x[ 9]*y[ 9] + x[10]*y[10] + x[11]*y[11] +
203 x[12]*y[12] + x[13]*y[13] + x[14]*y[14] + x[15]*y[15];
207 void FLAC__replaygain_synthesis__init_dither_context(DitherContext *d, int bits, int shapingtype)
209 static unsigned char default_dither [] = { 92, 92, 88, 84, 81, 78, 74, 67, 0, 0 };
210 static const float* F [] = { F44_0, F44_1, F44_2, F44_3 };
214 if (shapingtype < 0) shapingtype = 0;
215 if (shapingtype > 3) shapingtype = 3;
216 d->ShapingType = (NoiseShaping)shapingtype;
217 index = bits - 11 - shapingtype;
218 if (index < 0) index = 0;
219 if (index > 9) index = 9;
221 memset ( d->ErrorHistory , 0, sizeof (d->ErrorHistory ) );
222 memset ( d->DitherHistory, 0, sizeof (d->DitherHistory) );
224 d->FilterCoeff = F [shapingtype];
225 d->Mask = ((FLAC__uint64)-1) << (32 - bits);
226 d->Add = 0.5 * ((1L << (32 - bits)) - 1);
227 d->Dither = 0.01f*default_dither[index] / (((FLAC__int64)1) << bits);
228 d->LastHistoryIndex = 0;
232 * the following is based on parts of wavegain.c
235 static FLAC__int64 dither_output_(DitherContext *d, FLAC__bool do_dithering, int shapingtype, int i, double Sum, int k)
244 #define ROUND64(x) ( doubletmp.d = (x) + d->Add + (FLAC__int64)FLAC__I64L(0x001FFFFD80000000), doubletmp.i - (FLAC__int64)FLAC__I64L(0x433FFFFD80000000) )
247 if(shapingtype == 0) {
248 double tmp = random_equi_(d->Dither);
249 Sum2 = tmp - d->LastRandomNumber [k];
250 d->LastRandomNumber [k] = (int)tmp;
252 val = ROUND64(Sum2) & d->Mask;
255 Sum2 = random_triangular_(d->Dither) - scalar16_(d->DitherHistory[k], d->FilterCoeff + i);
256 Sum += d->DitherHistory [k] [(-1-i)&15] = (float)Sum2;
257 Sum2 = Sum + scalar16_(d->ErrorHistory [k], d->FilterCoeff + i);
258 val = ROUND64(Sum2) & d->Mask;
259 d->ErrorHistory [k] [(-1-i)&15] = (float)(Sum - val);
278 peak is in the range -32768.0 .. 32767.0
280 /* calculate factors for ReplayGain and ClippingPrevention */
281 *track_gain = GetTitleGain() + settings->man_gain;
282 scale = (float) pow(10., *track_gain * 0.05);
283 if(settings->clip_prev) {
284 factor_clip = (float) (32767./( peak + 1));
285 if(scale < factor_clip)
288 factor_clip /= scale;
289 scale *= factor_clip;
291 new_peak = (float) peak * scale;
293 dB = 20. * log10(scale);
294 *track_gain = (float) dB;
296 const double scale = pow(10., (double)gain * 0.05);
300 size_t FLAC__replaygain_synthesis__apply_gain(FLAC__byte *data_out, FLAC__bool little_endian_data_out, FLAC__bool unsigned_data_out, const FLAC__int32 * const input[], unsigned wide_samples, unsigned channels, const unsigned source_bps, const unsigned target_bps, const double scale, const FLAC__bool hard_limit, FLAC__bool do_dithering, DitherContext *dither_context)
302 static const FLAC__int32 conv_factors_[33] = {
303 -1, /* 0 bits-per-sample (not supported) */
304 -1, /* 1 bits-per-sample (not supported) */
305 -1, /* 2 bits-per-sample (not supported) */
306 -1, /* 3 bits-per-sample (not supported) */
307 268435456, /* 4 bits-per-sample */
308 134217728, /* 5 bits-per-sample */
309 67108864, /* 6 bits-per-sample */
310 33554432, /* 7 bits-per-sample */
311 16777216, /* 8 bits-per-sample */
312 8388608, /* 9 bits-per-sample */
313 4194304, /* 10 bits-per-sample */
314 2097152, /* 11 bits-per-sample */
315 1048576, /* 12 bits-per-sample */
316 524288, /* 13 bits-per-sample */
317 262144, /* 14 bits-per-sample */
318 131072, /* 15 bits-per-sample */
319 65536, /* 16 bits-per-sample */
320 32768, /* 17 bits-per-sample */
321 16384, /* 18 bits-per-sample */
322 8192, /* 19 bits-per-sample */
323 4096, /* 20 bits-per-sample */
324 2048, /* 21 bits-per-sample */
325 1024, /* 22 bits-per-sample */
326 512, /* 23 bits-per-sample */
327 256, /* 24 bits-per-sample */
328 128, /* 25 bits-per-sample */
329 64, /* 26 bits-per-sample */
330 32, /* 27 bits-per-sample */
331 16, /* 28 bits-per-sample */
332 8, /* 29 bits-per-sample */
333 4, /* 30 bits-per-sample */
334 2, /* 31 bits-per-sample */
335 1 /* 32 bits-per-sample */
337 static const FLAC__int64 hard_clip_factors_[33] = {
338 0, /* 0 bits-per-sample (not supported) */
339 0, /* 1 bits-per-sample (not supported) */
340 0, /* 2 bits-per-sample (not supported) */
341 0, /* 3 bits-per-sample (not supported) */
342 -8, /* 4 bits-per-sample */
343 -16, /* 5 bits-per-sample */
344 -32, /* 6 bits-per-sample */
345 -64, /* 7 bits-per-sample */
346 -128, /* 8 bits-per-sample */
347 -256, /* 9 bits-per-sample */
348 -512, /* 10 bits-per-sample */
349 -1024, /* 11 bits-per-sample */
350 -2048, /* 12 bits-per-sample */
351 -4096, /* 13 bits-per-sample */
352 -8192, /* 14 bits-per-sample */
353 -16384, /* 15 bits-per-sample */
354 -32768, /* 16 bits-per-sample */
355 -65536, /* 17 bits-per-sample */
356 -131072, /* 18 bits-per-sample */
357 -262144, /* 19 bits-per-sample */
358 -524288, /* 20 bits-per-sample */
359 -1048576, /* 21 bits-per-sample */
360 -2097152, /* 22 bits-per-sample */
361 -4194304, /* 23 bits-per-sample */
362 -8388608, /* 24 bits-per-sample */
363 -16777216, /* 25 bits-per-sample */
364 -33554432, /* 26 bits-per-sample */
365 -67108864, /* 27 bits-per-sample */
366 -134217728, /* 28 bits-per-sample */
367 -268435456, /* 29 bits-per-sample */
368 -536870912, /* 30 bits-per-sample */
369 -1073741824, /* 31 bits-per-sample */
370 (FLAC__int64)(-1073741824) * 2 /* 32 bits-per-sample */
372 const FLAC__int32 conv_factor = conv_factors_[target_bps];
373 const FLAC__int64 hard_clip_factor = hard_clip_factors_[target_bps];
375 * The integer input coming in has a varying range based on the
376 * source_bps. We want to normalize it to [-1.0, 1.0) so instead
377 * of doing two multiplies on each sample, we just multiple
378 * 'scale' by 1/(2^(source_bps-1))
380 const double multi_scale = scale / (double)(1u << (source_bps-1));
382 FLAC__byte * const start = data_out;
384 const FLAC__int32 *input_;
386 const unsigned bytes_per_sample = target_bps / 8;
387 const unsigned last_history_index = dither_context->LastHistoryIndex;
388 NoiseShaping noise_shaping = dither_context->ShapingType;
392 const FLAC__uint32 twiggle = 1u << (target_bps - 1);
394 FLAC__ASSERT(channels > 0 && channels <= FLAC_SHARE__MAX_SUPPORTED_CHANNELS);
395 FLAC__ASSERT(source_bps >= 4);
396 FLAC__ASSERT(target_bps >= 4);
397 FLAC__ASSERT(source_bps <= 32);
398 FLAC__ASSERT(target_bps < 32);
399 FLAC__ASSERT((target_bps & 7) == 0);
401 for(channel = 0; channel < channels; channel++) {
402 const unsigned incr = bytes_per_sample * channels;
403 data_out = start + bytes_per_sample * channel;
404 input_ = input[channel];
405 for(i = 0; i < wide_samples; i++, data_out += incr) {
406 sample = (double)input_[i] * multi_scale;
409 /* hard 6dB limiting */
411 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
412 else if(sample > 0.5)
413 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
415 sample *= 2147483647.f;
417 val64 = dither_output_(dither_context, do_dithering, noise_shaping, (i + last_history_index) % 32, sample, channel) / conv_factor;
419 val32 = (FLAC__int32)val64;
420 if(val64 >= -hard_clip_factor)
421 val32 = (FLAC__int32)(-(hard_clip_factor+1));
422 else if(val64 < hard_clip_factor)
423 val32 = (FLAC__int32)hard_clip_factor;
425 uval32 = (FLAC__uint32)val32;
426 if (unsigned_data_out)
429 if (little_endian_data_out) {
432 data_out[2] = (FLAC__byte)(uval32 >> 16);
435 data_out[1] = (FLAC__byte)(uval32 >> 8);
438 data_out[0] = (FLAC__byte)uval32;
445 data_out[0] = (FLAC__byte)(uval32 >> 16);
446 data_out[1] = (FLAC__byte)(uval32 >> 8);
447 data_out[2] = (FLAC__byte)uval32;
450 data_out[0] = (FLAC__byte)(uval32 >> 8);
451 data_out[1] = (FLAC__byte)uval32;
454 data_out[0] = (FLAC__byte)uval32;
460 dither_context->LastHistoryIndex = (last_history_index + wide_samples) % 32;
462 return wide_samples * channels * (target_bps/8);