1/* replaygain_synthesis - Routines for applying ReplayGain to a signal 2 * Copyright (C) 2002,2003,2004,2005,2006,2007 Josh Coalson 3 * 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. 8 * 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. 13 * 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 17 */ 18/* 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. 23 * 24 * The original copyright notices for WaveGain's dither.c and wavegain.c 25 * appear below: 26 */ 27/* 28 * (c) 2002 John Edwards 29 * mostly lifted from work by Frank Klemm 30 * random functions for dithering. 31 */ 32/* 33 * Copyright (C) 2002 John Edwards 34 * Additional code by Magnus Holmgren and Gian-Carlo Pascutto 35 */ 36 37#if HAVE_CONFIG_H 38# include <config.h> 39#endif 40 41#include <string.h> /* for memset() */ 42#include <math.h> 43#include "private/fast_float_math_hack.h" 44#include "replaygain_synthesis.h" 45#include "FLAC/assert.h" 46 47#ifndef FLaC__INLINE 48#define FLaC__INLINE 49#endif 50 51/* adjust for compilers that can't understand using LL suffix for int64_t literals */ 52#ifdef _MSC_VER 53#define FLAC__I64L(x) x 54#else 55#define FLAC__I64L(x) x##LL 56#endif 57 58 59/* 60 * the following is based on parts of dither.c 61 */ 62 63 64/* 65 * This is a simple random number generator with good quality for audio purposes. 66 * It consists of two polycounters with opposite rotation direction and different 67 * periods. The periods are coprime, so the total period is the product of both. 68 * 69 * ------------------------------------------------------------------------------------------------- 70 * +-> |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| 71 * | ------------------------------------------------------------------------------------------------- 72 * | | | | | | | 73 * | +--+--+--+-XOR-+--------+ 74 * | | 75 * +--------------------------------------------------------------------------------------+ 76 * 77 * ------------------------------------------------------------------------------------------------- 78 * |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| <-+ 79 * ------------------------------------------------------------------------------------------------- | 80 * | | | | | 81 * +--+----XOR----+--+ | 82 * | | 83 * +----------------------------------------------------------------------------------------+ 84 * 85 * 86 * The first has an period of 3*5*17*257*65537, the second of 7*47*73*178481, 87 * which gives a period of 18.410.713.077.675.721.215. The result is the 88 * XORed values of both generators. 89 */ 90 91static unsigned int random_int_(void) 92{ 93 static const unsigned char parity_[256] = { 94 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, 95 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, 96 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, 97 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, 98 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 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, 100 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, 101 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 102 }; 103 static unsigned int r1_ = 1; 104 static unsigned int r2_ = 1; 105 106 unsigned int t1, t2, t3, t4; 107 108 /* Parity calculation is done via table lookup, this is also available 109 * on CPUs without parity, can be implemented in C and avoid unpredictable 110 * jumps and slow rotate through the carry flag operations. 111 */ 112 t3 = t1 = r1_; t4 = t2 = r2_; 113 t1 &= 0xF5; t2 >>= 25; 114 t1 = parity_[t1]; t2 &= 0x63; 115 t1 <<= 31; t2 = parity_[t2]; 116 117 return (r1_ = (t3 >> 1) | t1 ) ^ (r2_ = (t4 + t4) | t2 ); 118} 119 120/* gives a equal distributed random number */ 121/* between -2^31*mult and +2^31*mult */ 122static double random_equi_(double mult) 123{ 124 return mult * (int) random_int_(); 125} 126 127/* gives a triangular distributed random number */ 128/* between -2^32*mult and +2^32*mult */ 129static double random_triangular_(double mult) 130{ 131 return mult * ( (double) (int) random_int_() + (double) (int) random_int_() ); 132} 133 134 135static const float F44_0 [16 + 32] = { 136 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, 137 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, 138 139 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, 140 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, 141 142 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, 143 (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0 144}; 145 146 147static const float F44_1 [16 + 32] = { /* SNR(w) = 4.843163 dB, SNR = -3.192134 dB */ 148 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833, 149 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967, 150 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116, 151 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024, 152 153 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833, 154 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967, 155 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116, 156 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024, 157 158 (float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833, 159 (float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967, 160 (float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116, 161 (float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024, 162}; 163 164 165static const float F44_2 [16 + 32] = { /* SNR(w) = 10.060213 dB, SNR = -12.766730 dB */ 166 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437, 167 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264, 168 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562, 169 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816, 170 171 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437, 172 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264, 173 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562, 174 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816, 175 176 (float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437, 177 (float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264, 178 (float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562, 179 (float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816, 180}; 181 182 183static const float F44_3 [16 + 32] = { /* SNR(w) = 15.382598 dB, SNR = -29.402334 dB */ 184 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515, 185 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785, 186 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927, 187 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099, 188 189 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515, 190 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785, 191 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927, 192 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099, 193 194 (float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515, 195 (float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785, 196 (float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927, 197 (float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099 198}; 199 200 201static double scalar16_(const float* x, const float* y) 202{ 203 return 204 x[ 0]*y[ 0] + x[ 1]*y[ 1] + x[ 2]*y[ 2] + x[ 3]*y[ 3] + 205 x[ 4]*y[ 4] + x[ 5]*y[ 5] + x[ 6]*y[ 6] + x[ 7]*y[ 7] + 206 x[ 8]*y[ 8] + x[ 9]*y[ 9] + x[10]*y[10] + x[11]*y[11] + 207 x[12]*y[12] + x[13]*y[13] + x[14]*y[14] + x[15]*y[15]; 208} 209 210 211void FLAC__replaygain_synthesis__init_dither_context(DitherContext *d, int bits, int shapingtype) 212{ 213 static unsigned char default_dither [] = { 92, 92, 88, 84, 81, 78, 74, 67, 0, 0 }; 214 static const float* F [] = { F44_0, F44_1, F44_2, F44_3 }; 215 216 int index; 217 218 if (shapingtype < 0) shapingtype = 0; 219 if (shapingtype > 3) shapingtype = 3; 220 d->ShapingType = (NoiseShaping)shapingtype; 221 index = bits - 11 - shapingtype; 222 if (index < 0) index = 0; 223 if (index > 9) index = 9; 224 225 memset ( d->ErrorHistory , 0, sizeof (d->ErrorHistory ) ); 226 memset ( d->DitherHistory, 0, sizeof (d->DitherHistory) ); 227 228 d->FilterCoeff = F [shapingtype]; 229 d->Mask = ((FLAC__uint64)-1) << (32 - bits); 230 d->Add = 0.5 * ((1L << (32 - bits)) - 1); 231 d->Dither = 0.01f*default_dither[index] / (((FLAC__int64)1) << bits); 232 d->LastHistoryIndex = 0; 233} 234 235/* 236 * the following is based on parts of wavegain.c 237 */ 238 239static FLaC__INLINE FLAC__int64 dither_output_(DitherContext *d, FLAC__bool do_dithering, int shapingtype, int i, double Sum, int k) 240{ 241 union { 242 double d; 243 FLAC__int64 i; 244 } doubletmp; 245 double Sum2; 246 FLAC__int64 val; 247 248#define ROUND64(x) ( doubletmp.d = (x) + d->Add + (FLAC__int64)FLAC__I64L(0x001FFFFD80000000), doubletmp.i - (FLAC__int64)FLAC__I64L(0x433FFFFD80000000) ) 249 250 if(do_dithering) { 251 if(shapingtype == 0) { 252 double tmp = random_equi_(d->Dither); 253 Sum2 = tmp - d->LastRandomNumber [k]; 254 d->LastRandomNumber [k] = (int)tmp; 255 Sum2 = Sum += Sum2; 256 val = ROUND64(Sum2) & d->Mask; 257 } 258 else { 259 Sum2 = random_triangular_(d->Dither) - scalar16_(d->DitherHistory[k], d->FilterCoeff + i); 260 Sum += d->DitherHistory [k] [(-1-i)&15] = (float)Sum2; 261 Sum2 = Sum + scalar16_(d->ErrorHistory [k], d->FilterCoeff + i); 262 val = ROUND64(Sum2) & d->Mask; 263 d->ErrorHistory [k] [(-1-i)&15] = (float)(Sum - val); 264 } 265 return val; 266 } 267 else 268 return ROUND64(Sum); 269 270#undef ROUND64 271} 272 273#if 0 274 float peak = 0.f, 275 new_peak, 276 factor_clip 277 double scale, 278 dB; 279 280 ... 281 282 peak is in the range -32768.0 .. 32767.0 283 284 /* calculate factors for ReplayGain and ClippingPrevention */ 285 *track_gain = GetTitleGain() + settings->man_gain; 286 scale = (float) pow(10., *track_gain * 0.05); 287 if(settings->clip_prev) { 288 factor_clip = (float) (32767./( peak + 1)); 289 if(scale < factor_clip) 290 factor_clip = 1.f; 291 else 292 factor_clip /= scale; 293 scale *= factor_clip; 294 } 295 new_peak = (float) peak * scale; 296 297 dB = 20. * log10(scale); 298 *track_gain = (float) dB; 299 300 const double scale = pow(10., (double)gain * 0.05); 301#endif 302 303 304size_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) 305{ 306 static const FLAC__int32 conv_factors_[33] = { 307 -1, /* 0 bits-per-sample (not supported) */ 308 -1, /* 1 bits-per-sample (not supported) */ 309 -1, /* 2 bits-per-sample (not supported) */ 310 -1, /* 3 bits-per-sample (not supported) */ 311 268435456, /* 4 bits-per-sample */ 312 134217728, /* 5 bits-per-sample */ 313 67108864, /* 6 bits-per-sample */ 314 33554432, /* 7 bits-per-sample */ 315 16777216, /* 8 bits-per-sample */ 316 8388608, /* 9 bits-per-sample */ 317 4194304, /* 10 bits-per-sample */ 318 2097152, /* 11 bits-per-sample */ 319 1048576, /* 12 bits-per-sample */ 320 524288, /* 13 bits-per-sample */ 321 262144, /* 14 bits-per-sample */ 322 131072, /* 15 bits-per-sample */ 323 65536, /* 16 bits-per-sample */ 324 32768, /* 17 bits-per-sample */ 325 16384, /* 18 bits-per-sample */ 326 8192, /* 19 bits-per-sample */ 327 4096, /* 20 bits-per-sample */ 328 2048, /* 21 bits-per-sample */ 329 1024, /* 22 bits-per-sample */ 330 512, /* 23 bits-per-sample */ 331 256, /* 24 bits-per-sample */ 332 128, /* 25 bits-per-sample */ 333 64, /* 26 bits-per-sample */ 334 32, /* 27 bits-per-sample */ 335 16, /* 28 bits-per-sample */ 336 8, /* 29 bits-per-sample */ 337 4, /* 30 bits-per-sample */ 338 2, /* 31 bits-per-sample */ 339 1 /* 32 bits-per-sample */ 340 }; 341 static const FLAC__int64 hard_clip_factors_[33] = { 342 0, /* 0 bits-per-sample (not supported) */ 343 0, /* 1 bits-per-sample (not supported) */ 344 0, /* 2 bits-per-sample (not supported) */ 345 0, /* 3 bits-per-sample (not supported) */ 346 -8, /* 4 bits-per-sample */ 347 -16, /* 5 bits-per-sample */ 348 -32, /* 6 bits-per-sample */ 349 -64, /* 7 bits-per-sample */ 350 -128, /* 8 bits-per-sample */ 351 -256, /* 9 bits-per-sample */ 352 -512, /* 10 bits-per-sample */ 353 -1024, /* 11 bits-per-sample */ 354 -2048, /* 12 bits-per-sample */ 355 -4096, /* 13 bits-per-sample */ 356 -8192, /* 14 bits-per-sample */ 357 -16384, /* 15 bits-per-sample */ 358 -32768, /* 16 bits-per-sample */ 359 -65536, /* 17 bits-per-sample */ 360 -131072, /* 18 bits-per-sample */ 361 -262144, /* 19 bits-per-sample */ 362 -524288, /* 20 bits-per-sample */ 363 -1048576, /* 21 bits-per-sample */ 364 -2097152, /* 22 bits-per-sample */ 365 -4194304, /* 23 bits-per-sample */ 366 -8388608, /* 24 bits-per-sample */ 367 -16777216, /* 25 bits-per-sample */ 368 -33554432, /* 26 bits-per-sample */ 369 -67108864, /* 27 bits-per-sample */ 370 -134217728, /* 28 bits-per-sample */ 371 -268435456, /* 29 bits-per-sample */ 372 -536870912, /* 30 bits-per-sample */ 373 -1073741824, /* 31 bits-per-sample */ 374 (FLAC__int64)(-1073741824) * 2 /* 32 bits-per-sample */ 375 }; 376 const FLAC__int32 conv_factor = conv_factors_[target_bps]; 377 const FLAC__int64 hard_clip_factor = hard_clip_factors_[target_bps]; 378 /* 379 * The integer input coming in has a varying range based on the 380 * source_bps. We want to normalize it to [-1.0, 1.0) so instead 381 * of doing two multiplies on each sample, we just multiple 382 * 'scale' by 1/(2^(source_bps-1)) 383 */ 384 const double multi_scale = scale / (double)(1u << (source_bps-1)); 385 386 FLAC__byte * const start = data_out; 387 unsigned i, channel; 388 const FLAC__int32 *input_; 389 double sample; 390 const unsigned bytes_per_sample = target_bps / 8; 391 const unsigned last_history_index = dither_context->LastHistoryIndex; 392 NoiseShaping noise_shaping = dither_context->ShapingType; 393 FLAC__int64 val64; 394 FLAC__int32 val32; 395 FLAC__int32 uval32; 396 const FLAC__uint32 twiggle = 1u << (target_bps - 1); 397 398 FLAC__ASSERT(channels > 0 && channels <= FLAC_SHARE__MAX_SUPPORTED_CHANNELS); 399 FLAC__ASSERT(source_bps >= 4); 400 FLAC__ASSERT(target_bps >= 4); 401 FLAC__ASSERT(source_bps <= 32); 402 FLAC__ASSERT(target_bps < 32); 403 FLAC__ASSERT((target_bps & 7) == 0); 404 405 for(channel = 0; channel < channels; channel++) { 406 const unsigned incr = bytes_per_sample * channels; 407 data_out = start + bytes_per_sample * channel; 408 input_ = input[channel]; 409 for(i = 0; i < wide_samples; i++, data_out += incr) { 410 sample = (double)input_[i] * multi_scale; 411 412 if(hard_limit) { 413 /* hard 6dB limiting */ 414 if(sample < -0.5) 415 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5; 416 else if(sample > 0.5) 417 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5; 418 } 419 sample *= 2147483647.f; 420 421 val64 = dither_output_(dither_context, do_dithering, noise_shaping, (i + last_history_index) % 32, sample, channel) / conv_factor; 422 423 val32 = (FLAC__int32)val64; 424 if(val64 >= -hard_clip_factor) 425 val32 = (FLAC__int32)(-(hard_clip_factor+1)); 426 else if(val64 < hard_clip_factor) 427 val32 = (FLAC__int32)hard_clip_factor; 428 429 uval32 = (FLAC__uint32)val32; 430 if (unsigned_data_out) 431 uval32 ^= twiggle; 432 433 if (little_endian_data_out) { 434 switch(target_bps) { 435 case 24: 436 data_out[2] = (FLAC__byte)(uval32 >> 16); 437 /* fall through */ 438 case 16: 439 data_out[1] = (FLAC__byte)(uval32 >> 8); 440 /* fall through */ 441 case 8: 442 data_out[0] = (FLAC__byte)uval32; 443 break; 444 } 445 } 446 else { 447 switch(target_bps) { 448 case 24: 449 data_out[0] = (FLAC__byte)(uval32 >> 16); 450 data_out[1] = (FLAC__byte)(uval32 >> 8); 451 data_out[2] = (FLAC__byte)uval32; 452 break; 453 case 16: 454 data_out[0] = (FLAC__byte)(uval32 >> 8); 455 data_out[1] = (FLAC__byte)uval32; 456 break; 457 case 8: 458 data_out[0] = (FLAC__byte)uval32; 459 break; 460 } 461 } 462 } 463 } 464 dither_context->LastHistoryIndex = (last_history_index + wide_samples) % 32; 465 466 return wide_samples * channels * (target_bps/8); 467} 468