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