1 /* plugin_common - Routines common to several plugins
2 * Copyright (C) 2002 Josh Coalson
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.
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.
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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, 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
37 #include <string.h> /* for memset() */
39 #include "private/fast_float_math_hack.h"
40 #include "replaygain_synthesis.h"
41 #include "FLAC/assert.h"
44 #define FLAC__INLINE __inline
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_()
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__plugin_common__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 index = bits - 11 - shapingtype;
212 if (index < 0) index = 0;
213 if (index > 9) index = 9;
215 memset ( d->ErrorHistory , 0, sizeof (d->ErrorHistory ) );
216 memset ( d->DitherHistory, 0, sizeof (d->DitherHistory) );
218 d->FilterCoeff = F [shapingtype];
219 d->Mask = ((FLAC__uint64)-1) << (32 - bits);
220 d->Add = 0.5 * ((1L << (32 - bits)) - 1);
221 d->Dither = 0.01*default_dither[index] / (((FLAC__int64)1) << bits);
225 * the following is based on parts of wavegain.c
228 static FLAC__INLINE FLAC__int64 dither_output_(DitherContext *d, FLAC__bool do_dithering, int shapingtype, int i, double Sum, int k)
230 double doubletmp, Sum2;
233 #define ROUND64(x) ( doubletmp = (x) + d->Add + (FLAC__int64)0x001FFFFD80000000L, *(FLAC__int64*)(&doubletmp) - (FLAC__int64)0x433FFFFD80000000L )
236 if(shapingtype == 0) {
237 double tmp = random_equi_(d->Dither);
238 Sum2 = tmp - d->LastRandomNumber [k];
239 d->LastRandomNumber [k] = tmp;
241 val = ROUND64(Sum2) & d->Mask;
244 Sum2 = random_triangular_(d->Dither) - scalar16_(d->DitherHistory[k], d->FilterCoeff + i);
245 Sum += d->DitherHistory [k] [(-1-i)&15] = Sum2;
246 Sum2 = Sum + scalar16_(d->ErrorHistory [k], d->FilterCoeff + i);
247 val = ROUND64(Sum2) & d->Mask;
248 d->ErrorHistory [k] [(-1-i)&15] = (float)(Sum - val);
267 peak is in the range -32768.0 .. 32767.0
269 /* calculate factors for ReplayGain and ClippingPrevention */
270 *track_gain = GetTitleGain() + settings->man_gain;
271 scale = (float) pow(10., *track_gain * 0.05);
272 if(settings->clip_prev) {
273 factor_clip = (float) (32767./( peak + 1));
274 if(scale < factor_clip)
277 factor_clip /= scale;
278 scale *= factor_clip;
280 new_peak = (float) peak * scale;
282 dB = 20. * log10(scale);
283 *track_gain = (float) dB;
285 const double scale = (float) pow(10., (double)gain * 0.05); /*@@@@ why downcast pow() output to float? */
289 PLUGIN_COMMON_API int FLAC__plugin_common__apply_gain(FLAC__byte *data_out, FLAC__int32 *input, unsigned wide_samples, unsigned channels, const unsigned source_bps, const unsigned target_bps, const float scale, const FLAC__bool hard_limit, FLAC__bool do_dithering, NoiseShaping noise_shaping, DitherContext *dither_context)
291 static const FLAC__int32 conv_factors_[33] = {
292 -1, /* 0 bits-per-sample (not supported) */
293 -1, /* 1 bits-per-sample (not supported) */
294 -1, /* 2 bits-per-sample (not supported) */
295 -1, /* 3 bits-per-sample (not supported) */
296 268435456, /* 4 bits-per-sample */
297 134217728, /* 5 bits-per-sample */
298 67108864, /* 6 bits-per-sample */
299 33554432, /* 7 bits-per-sample */
300 16777216, /* 8 bits-per-sample */
301 8388608, /* 9 bits-per-sample */
302 4194304, /* 10 bits-per-sample */
303 2097152, /* 11 bits-per-sample */
304 1048576, /* 12 bits-per-sample */
305 524288, /* 13 bits-per-sample */
306 262144, /* 14 bits-per-sample */
307 131072, /* 15 bits-per-sample */
308 65536, /* 16 bits-per-sample */
309 32768, /* 17 bits-per-sample */
310 16384, /* 18 bits-per-sample */
311 8192, /* 19 bits-per-sample */
312 4096, /* 20 bits-per-sample */
313 2048, /* 21 bits-per-sample */
314 1024, /* 22 bits-per-sample */
315 512, /* 23 bits-per-sample */
316 256, /* 24 bits-per-sample */
317 128, /* 25 bits-per-sample */
318 64, /* 26 bits-per-sample */
319 32, /* 27 bits-per-sample */
320 16, /* 28 bits-per-sample */
321 8, /* 29 bits-per-sample */
322 4, /* 30 bits-per-sample */
323 2, /* 31 bits-per-sample */
324 1 /* 32 bits-per-sample */
326 static const FLAC__int64 hard_clip_factors_[33] = {
327 0, /* 0 bits-per-sample (not supported) */
328 0, /* 1 bits-per-sample (not supported) */
329 0, /* 2 bits-per-sample (not supported) */
330 0, /* 3 bits-per-sample (not supported) */
331 -8, /* 4 bits-per-sample */
332 -16, /* 5 bits-per-sample */
333 -32, /* 6 bits-per-sample */
334 -64, /* 7 bits-per-sample */
335 -128, /* 8 bits-per-sample */
336 -256, /* 9 bits-per-sample */
337 -512, /* 10 bits-per-sample */
338 -1024, /* 11 bits-per-sample */
339 -2048, /* 12 bits-per-sample */
340 -4096, /* 13 bits-per-sample */
341 -8192, /* 14 bits-per-sample */
342 -16384, /* 15 bits-per-sample */
343 -32768, /* 16 bits-per-sample */
344 -65536, /* 17 bits-per-sample */
345 -131072, /* 18 bits-per-sample */
346 -262144, /* 19 bits-per-sample */
347 -524288, /* 20 bits-per-sample */
348 -1048576, /* 21 bits-per-sample */
349 -2097152, /* 22 bits-per-sample */
350 -4194304, /* 23 bits-per-sample */
351 -8388608, /* 24 bits-per-sample */
352 -16777216, /* 25 bits-per-sample */
353 -33554432, /* 26 bits-per-sample */
354 -67108864, /* 27 bits-per-sample */
355 -134217728, /* 28 bits-per-sample */
356 -268435456, /* 29 bits-per-sample */
357 -536870912, /* 30 bits-per-sample */
358 -1073741824, /* 31 bits-per-sample */
359 (FLAC__int64)(-1073741824) * 2 /* 32 bits-per-sample */
361 const FLAC__int32 conv_factor = conv_factors_[source_bps];
362 const FLAC__int64 hard_clip_factor = hard_clip_factors_[source_bps];
364 * The integer input coming in has a varying range based on the
365 * source_bps. We want to normalize it to [-1.0, 1.0) so instead
366 * of doing two multiplies on each sample, we just multiple
367 * 'scale' by 1/(2^(source_bps-1))
369 const double multi_scale = scale / (double)(1u << (source_bps-1));
371 FLAC__byte * const start = data_out;
372 const unsigned samples = wide_samples * channels;
373 #ifdef FLAC__PLUGIN_COMMON__DONT_UNROLL
374 const unsigned dither_twiggle = channels - 1;
375 unsigned dither_source = 0;
381 FLAC__ASSERT(FLAC_PLUGIN__MAX_SUPPORTED_CHANNELS == 2);
382 FLAC__ASSERT(channels > 0 && channels <= FLAC_PLUGIN__MAX_SUPPORTED_CHANNELS);
383 FLAC__ASSERT(source_bps >= 4);
384 FLAC__ASSERT(target_bps >= 4);
385 FLAC__ASSERT(source_bps <= 32);
386 FLAC__ASSERT(target_bps < 32);
387 FLAC__ASSERT((target_bps & 7) == 0);
389 #ifdef FLAC__PLUGIN_COMMON__DONT_UNROLL
391 * This flavor handles 1 or 2 channels with the same code
394 for(i = 0; i < samples; i++, coeff++) {
395 sample = (double)input[i] * multi_scale;
398 /* hard 6dB limiting */
400 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
401 else if(sample > 0.5)
402 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
404 sample *= 2147483647.f;
410 if(coeff >= (32<<dither_twiggle))
413 /* 'coeff>>dither_twiggle' is the same as 'coeff/channels' */
414 val64 = dither_output_(dither_context, do_dithering, noise_shaping, coeff>>dither_twiggle, sample, dither_source) / conv_factor;
416 dither_source ^= dither_twiggle;
418 val32 = (FLAC__int32)val64;
419 if(val64 >= -hard_clip_factor)
420 val32 = (FLAC__int32)(-(hard_clip_factor+1));
421 else if(val64 < hard_clip_factor)
422 val32 = (FLAC__int32)hard_clip_factor;
426 data_out[0] = val32 ^ 0x80;
429 data_out[2] = (FLAC__byte)(val32 >> 16);
432 data_out[1] = (FLAC__byte)(val32 >> 8);
433 data_out[0] = (FLAC__byte)val32;
437 data_out += target_bps/8;
441 * This flavor has optimized versions for 1 or 2 channels
448 for(i = 0; i < samples; ) {
449 sample = (double)input[i] * multi_scale;
452 /* hard 6dB limiting */
454 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
455 else if(sample > 0.5)
456 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
458 sample *= 2147483647.f;
460 val64 = dither_output_(dither_context, do_dithering, noise_shaping, coeff, sample, 0) / conv_factor;
462 val32 = (FLAC__int32)val64;
463 if(val64 >= -hard_clip_factor)
464 val32 = (FLAC__int32)(-(hard_clip_factor+1));
465 else if(val64 < hard_clip_factor)
466 val32 = (FLAC__int32)hard_clip_factor;
470 data_out[0] = val32 ^ 0x80;
473 data_out[2] = (FLAC__byte)(val32 >> 16);
476 data_out[1] = (FLAC__byte)(val32 >> 8);
477 data_out[0] = (FLAC__byte)val32;
480 data_out += target_bps/8;
484 sample = (double)input[i] * multi_scale;
487 /* hard 6dB limiting */
489 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
490 else if(sample > 0.5)
491 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
493 sample *= 2147483647.f;
495 val64 = dither_output_(dither_context, do_dithering, noise_shaping, coeff, sample, 1) / conv_factor;
497 val32 = (FLAC__int32)val64;
498 if(val64 >= -hard_clip_factor)
499 val32 = (FLAC__int32)(-(hard_clip_factor+1));
500 else if(val64 < hard_clip_factor)
501 val32 = (FLAC__int32)hard_clip_factor;
505 data_out[0] = val32 ^ 0x80;
508 data_out[2] = (FLAC__byte)(val32 >> 16);
511 data_out[1] = (FLAC__byte)(val32 >> 8);
512 data_out[0] = (FLAC__byte)val32;
515 data_out += target_bps/8;
528 for(i = 0; i < samples; i++, coeff++) {
532 sample = (double)input[i] * multi_scale;
535 /* hard 6dB limiting */
537 sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
538 else if(sample > 0.5)
539 sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
541 sample *= 2147483647.f;
543 val64 = dither_output_(dither_context, do_dithering, noise_shaping, coeff, sample, 0) / conv_factor;
545 val32 = (FLAC__int32)val64;
546 if(val64 >= -hard_clip_factor)
547 val32 = (FLAC__int32)(-(hard_clip_factor+1));
548 else if(val64 < hard_clip_factor)
549 val32 = (FLAC__int32)hard_clip_factor;
553 data_out[0] = val32 ^ 0x80;
556 data_out[2] = (FLAC__byte)(val32 >> 16);
559 data_out[1] = (FLAC__byte)(val32 >> 8);
560 data_out[0] = (FLAC__byte)val32;
563 data_out += target_bps/8;
568 return data_out - start;