1 /* replaygain_synthesis - Routines for applying ReplayGain to a signal
2 * Copyright (C) 2002-2009 Josh Coalson
3 * Copyright (C) 2011-2016 Xiph.Org Foundation
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2.1 of the License, or (at your option) any later version.
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with this library; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 * This is an aggregation of pieces of code from John Edwards' WaveGain
21 * program. Mostly cosmetic changes were made; otherwise, the dithering
22 * code is almost untouched and the gain processing was converted from
23 * processing a whole file to processing chunks of samples.
25 * The original copyright notices for WaveGain's dither.c and wavegain.c
29 * (c) 2002 John Edwards
30 * mostly lifted from work by Frank Klemm
31 * random functions for dithering.
34 * Copyright (C) 2002 John Edwards
35 * Additional code by Magnus Holmgren and Gian-Carlo Pascutto
42 #include <string.h> /* for memset() */
44 #include "share/replaygain_synthesis.h"
45 #include "FLAC/assert.h"
47 #define FLAC__I64L(x) x##LL
51 * the following is based on parts of dither.c
56 * This is a simple random number generator with good quality for audio purposes.
57 * It consists of two polycounters with opposite rotation direction and different
58 * periods. The periods are coprime, so the total period is the product of both.
60 * -------------------------------------------------------------------------------------------------
61 * +-> |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|
62 * | -------------------------------------------------------------------------------------------------
64 * | +--+--+--+-XOR-+--------+
66 * +--------------------------------------------------------------------------------------+
68 * -------------------------------------------------------------------------------------------------
69 * |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| <-+
70 * ------------------------------------------------------------------------------------------------- |
72 * +--+----XOR----+--+ |
74 * +----------------------------------------------------------------------------------------+
77 * The first has an period of 3*5*17*257*65537, the second of 7*47*73*178481,
78 * which gives a period of 18.410.713.077.675.721.215. The result is the
79 * XORed values of both generators.
82 static unsigned int random_int_(void)
84 static const unsigned char parity_[256] = {
85 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,
86 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,
87 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,
88 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,
89 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,
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 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,
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
94 static unsigned int r1_ = 1;
95 static unsigned int r2_ = 1;
97 unsigned int t1, t2, t3, t4;
99 /* Parity calculation is done via table lookup, this is also available
100 * on CPUs without parity, can be implemented in C and avoid unpredictable
101 * jumps and slow rotate through the carry flag operations.
103 t3 = t1 = r1_; t4 = t2 = r2_;
104 t1 &= 0xF5; t2 >>= 25;
105 t1 = parity_[t1]; t2 &= 0x63;
106 t1 <<= 31; t2 = parity_[t2];
108 return (r1_ = (t3 >> 1) | t1 ) ^ (r2_ = (t4 + t4) | t2 );
111 /* gives a equal distributed random number */
112 /* between -2^31*mult and +2^31*mult */
113 static double random_equi_(double mult)
115 return mult * (int) random_int_();
118 /* gives a triangular distributed random number */
119 /* between -2^32*mult and +2^32*mult */
120 static double random_triangular_(double mult)
122 return mult * ( (double) (int) random_int_() + (double) (int) random_int_() );
126 static const float F44_0 [16 + 32] = {
127 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
128 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
130 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
131 (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,
134 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0
138 static const float F44_1 [16 + 32] = { /* SNR(w) = 4.843163 dB, SNR = -3.192134 dB */
139 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
140 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
141 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
142 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,
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,
156 static const float F44_2 [16 + 32] = { /* SNR(w) = 10.060213 dB, SNR = -12.766730 dB */
157 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
158 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
159 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
160 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,
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,
174 static const float F44_3 [16 + 32] = { /* SNR(w) = 15.382598 dB, SNR = -29.402334 dB */
175 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
176 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
177 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
178 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099,
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
192 static double scalar16_(const float* x, const float* y)
195 x[ 0]*y[ 0] + x[ 1]*y[ 1] + x[ 2]*y[ 2] + x[ 3]*y[ 3] +
196 x[ 4]*y[ 4] + x[ 5]*y[ 5] + x[ 6]*y[ 6] + x[ 7]*y[ 7] +
197 x[ 8]*y[ 8] + x[ 9]*y[ 9] + x[10]*y[10] + x[11]*y[11] +
198 x[12]*y[12] + x[13]*y[13] + x[14]*y[14] + x[15]*y[15];
202 void FLAC__replaygain_synthesis__init_dither_context(DitherContext *d, int bits, int shapingtype)
204 static unsigned char default_dither [] = { 92, 92, 88, 84, 81, 78, 74, 67, 0, 0 };
205 static const float* F [] = { F44_0, F44_1, F44_2, F44_3 };
209 if (shapingtype < 0) shapingtype = 0;
210 if (shapingtype > 3) shapingtype = 3;
211 d->ShapingType = (NoiseShaping)shapingtype;
212 indx = bits - 11 - shapingtype;
213 if (indx < 0) indx = 0;
214 if (indx > 9) indx = 9;
216 memset ( d->ErrorHistory , 0, sizeof (d->ErrorHistory ) );
217 memset ( d->DitherHistory, 0, sizeof (d->DitherHistory) );
219 d->FilterCoeff = F [shapingtype];
220 d->Mask = ((FLAC__uint64)-1) << (32 - bits);
221 d->Add = 0.5 * ((1L << (32 - bits)) - 1);
222 d->Dither = 0.01f*default_dither[indx] / (((FLAC__int64)1) << bits);
223 d->LastHistoryIndex = 0;
227 * the following is based on parts of wavegain.c
230 static FLAC__int64 dither_output_(DitherContext *d, FLAC__bool do_dithering, int shapingtype, int i, double Sum, int k)
239 #define ROUND64(x) ( doubletmp.d = (x) + d->Add + (FLAC__int64)FLAC__I64L(0x001FFFFD80000000), doubletmp.i - (FLAC__int64)FLAC__I64L(0x433FFFFD80000000) )
242 if(shapingtype == 0) {
243 double tmp = random_equi_(d->Dither);
244 Sum2 = tmp - d->LastRandomNumber [k];
245 d->LastRandomNumber [k] = (int)tmp;
247 val = ROUND64(Sum2) & d->Mask;
250 Sum2 = random_triangular_(d->Dither) - scalar16_(d->DitherHistory[k], d->FilterCoeff + i);
251 Sum += d->DitherHistory [k] [(-1-i)&15] = (float)Sum2;
252 Sum2 = Sum + scalar16_(d->ErrorHistory [k], d->FilterCoeff + i);
253 val = ROUND64(Sum2) & d->Mask;
254 d->ErrorHistory [k] [(-1-i)&15] = (float)(Sum - val);
273 peak is in the range -32768.0 .. 32767.0
275 /* calculate factors for ReplayGain and ClippingPrevention */
276 *track_gain = GetTitleGain() + settings->man_gain;
277 scale = (float) pow(10., *track_gain * 0.05);
278 if(settings->clip_prev) {
279 factor_clip = (float) (32767./( peak + 1));
280 if(scale < factor_clip)
283 factor_clip /= scale;
284 scale *= factor_clip;
286 new_peak = (float) peak * scale;
288 dB = 20. * log10(scale);
289 *track_gain = (float) dB;
291 const double scale = pow(10., (double)gain * 0.05);
295 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)
297 static const FLAC__int64 hard_clip_factors_[33] = {
298 0, /* 0 bits-per-sample (not supported) */
299 0, /* 1 bits-per-sample (not supported) */
300 0, /* 2 bits-per-sample (not supported) */
301 0, /* 3 bits-per-sample (not supported) */
302 -8, /* 4 bits-per-sample */
303 -16, /* 5 bits-per-sample */
304 -32, /* 6 bits-per-sample */
305 -64, /* 7 bits-per-sample */
306 -128, /* 8 bits-per-sample */
307 -256, /* 9 bits-per-sample */
308 -512, /* 10 bits-per-sample */
309 -1024, /* 11 bits-per-sample */
310 -2048, /* 12 bits-per-sample */
311 -4096, /* 13 bits-per-sample */
312 -8192, /* 14 bits-per-sample */
313 -16384, /* 15 bits-per-sample */
314 -32768, /* 16 bits-per-sample */
315 -65536, /* 17 bits-per-sample */
316 -131072, /* 18 bits-per-sample */
317 -262144, /* 19 bits-per-sample */
318 -524288, /* 20 bits-per-sample */
319 -1048576, /* 21 bits-per-sample */
320 -2097152, /* 22 bits-per-sample */
321 -4194304, /* 23 bits-per-sample */
322 -8388608, /* 24 bits-per-sample */
323 -16777216, /* 25 bits-per-sample */
324 -33554432, /* 26 bits-per-sample */
325 -67108864, /* 27 bits-per-sample */
326 -134217728, /* 28 bits-per-sample */
327 -268435456, /* 29 bits-per-sample */
328 -536870912, /* 30 bits-per-sample */
329 -1073741824, /* 31 bits-per-sample */
330 (FLAC__int64)(-1073741824) * 2 /* 32 bits-per-sample */
332 const FLAC__int32 conv_shift = 32 - target_bps;
333 const FLAC__int64 hard_clip_factor = hard_clip_factors_[target_bps];
335 * The integer input coming in has a varying range based on the
336 * source_bps. We want to normalize it to [-1.0, 1.0) so instead
337 * of doing two multiplies on each sample, we just multiple
338 * 'scale' by 1/(2^(source_bps-1))
340 const double multi_scale = scale / (double)(1u << (source_bps-1));
342 FLAC__byte * const start = data_out;
344 const FLAC__int32 *input_;
346 const unsigned bytes_per_sample = target_bps / 8;
347 const unsigned last_history_index = dither_context->LastHistoryIndex;
348 NoiseShaping noise_shaping = dither_context->ShapingType;
352 const FLAC__uint32 twiggle = 1u << (target_bps - 1);
354 FLAC__ASSERT(channels > 0 && channels <= FLAC_SHARE__MAX_SUPPORTED_CHANNELS);
355 FLAC__ASSERT(source_bps >= 4);
356 FLAC__ASSERT(target_bps >= 4);
357 FLAC__ASSERT(source_bps <= 32);
358 FLAC__ASSERT(target_bps < 32);
359 FLAC__ASSERT((target_bps & 7) == 0);
361 for(channel = 0; channel < channels; channel++) {
362 const unsigned incr = bytes_per_sample * channels;
363 data_out = start + bytes_per_sample * channel;
364 input_ = input[channel];
365 for(i = 0; i < wide_samples; i++, data_out += incr) {
366 sample = (double)input_[i] * multi_scale;
369 /* hard 6dB limiting */
371 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
372 else if(sample > 0.5)
373 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
375 sample *= 2147483647.;
377 val64 = dither_output_(dither_context, do_dithering, noise_shaping, (i + last_history_index) % 32, sample, channel) >> conv_shift;
379 val32 = (FLAC__int32)val64;
380 if(val64 >= -hard_clip_factor)
381 val32 = (FLAC__int32)(-(hard_clip_factor+1));
382 else if(val64 < hard_clip_factor)
383 val32 = (FLAC__int32)hard_clip_factor;
385 uval32 = (FLAC__uint32)val32;
386 if (unsigned_data_out)
389 if (little_endian_data_out) {
392 data_out[2] = (FLAC__byte)(uval32 >> 16);
395 data_out[1] = (FLAC__byte)(uval32 >> 8);
398 data_out[0] = (FLAC__byte)uval32;
405 data_out[0] = (FLAC__byte)(uval32 >> 16);
406 data_out[1] = (FLAC__byte)(uval32 >> 8);
407 data_out[2] = (FLAC__byte)uval32;
410 data_out[0] = (FLAC__byte)(uval32 >> 8);
411 data_out[1] = (FLAC__byte)uval32;
414 data_out[0] = (FLAC__byte)uval32;
420 dither_context->LastHistoryIndex = (last_history_index + wide_samples) % 32;
422 return wide_samples * channels * (target_bps/8);