Remove _MSC_VER specific FLAC__I64L definition.
[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-2013  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 #if HAVE_CONFIG_H
39 #  include <config.h>
40 #endif
41
42 #include <string.h> /* for memset() */
43 #include <math.h>
44 #include "private/fast_float_math_hack.h"
45 #include "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 unsigned int random_int_(void)
84 {
85         static const unsigned char 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 unsigned int r1_ = 1;
96         static unsigned int r2_ = 1;
97
98         unsigned int 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 unsigned char 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 /*
228  * the following is based on parts of wavegain.c
229  */
230
231 static FLAC__int64 dither_output_(DitherContext *d, FLAC__bool do_dithering, int shapingtype, int i, double Sum, int k)
232 {
233         union {
234                 double d;
235                 FLAC__int64 i;
236         } doubletmp;
237         double Sum2;
238         FLAC__int64 val;
239
240 #define ROUND64(x)   ( doubletmp.d = (x) + d->Add + (FLAC__int64)FLAC__I64L(0x001FFFFD80000000), doubletmp.i - (FLAC__int64)FLAC__I64L(0x433FFFFD80000000) )
241
242         if(do_dithering) {
243                 if(shapingtype == 0) {
244                         double  tmp = random_equi_(d->Dither);
245                         Sum2 = tmp - d->LastRandomNumber [k];
246                         d->LastRandomNumber [k] = (int)tmp;
247                         Sum2 = Sum += Sum2;
248                         val = ROUND64(Sum2) & d->Mask;
249                 }
250                 else {
251                         Sum2 = random_triangular_(d->Dither) - scalar16_(d->DitherHistory[k], d->FilterCoeff + i);
252                         Sum += d->DitherHistory [k] [(-1-i)&15] = (float)Sum2;
253                         Sum2 = Sum + scalar16_(d->ErrorHistory [k], d->FilterCoeff + i);
254                         val = ROUND64(Sum2) & d->Mask;
255                         d->ErrorHistory [k] [(-1-i)&15] = (float)(Sum - val);
256                 }
257                 return val;
258         }
259         else
260                 return ROUND64(Sum);
261
262 #undef ROUND64
263 }
264
265 #if 0
266         float        peak = 0.f,
267                      new_peak,
268                      factor_clip
269         double       scale,
270                      dB;
271
272         ...
273
274         peak is in the range -32768.0 .. 32767.0
275
276         /* calculate factors for ReplayGain and ClippingPrevention */
277         *track_gain = GetTitleGain() + settings->man_gain;
278         scale = (float) pow(10., *track_gain * 0.05);
279         if(settings->clip_prev) {
280                 factor_clip  = (float) (32767./( peak + 1));
281                 if(scale < factor_clip)
282                         factor_clip = 1.f;
283                 else
284                         factor_clip /= scale;
285                 scale *= factor_clip;
286         }
287         new_peak = (float) peak * scale;
288
289         dB = 20. * log10(scale);
290         *track_gain = (float) dB;
291
292         const double scale = pow(10., (double)gain * 0.05);
293 #endif
294
295
296 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 {
298         static const FLAC__int32 conv_factors_[33] = {
299                 -1, /* 0 bits-per-sample (not supported) */
300                 -1, /* 1 bits-per-sample (not supported) */
301                 -1, /* 2 bits-per-sample (not supported) */
302                 -1, /* 3 bits-per-sample (not supported) */
303                 268435456, /* 4 bits-per-sample */
304                 134217728, /* 5 bits-per-sample */
305                 67108864, /* 6 bits-per-sample */
306                 33554432, /* 7 bits-per-sample */
307                 16777216, /* 8 bits-per-sample */
308                 8388608, /* 9 bits-per-sample */
309                 4194304, /* 10 bits-per-sample */
310                 2097152, /* 11 bits-per-sample */
311                 1048576, /* 12 bits-per-sample */
312                 524288, /* 13 bits-per-sample */
313                 262144, /* 14 bits-per-sample */
314                 131072, /* 15 bits-per-sample */
315                 65536, /* 16 bits-per-sample */
316                 32768, /* 17 bits-per-sample */
317                 16384, /* 18 bits-per-sample */
318                 8192, /* 19 bits-per-sample */
319                 4096, /* 20 bits-per-sample */
320                 2048, /* 21 bits-per-sample */
321                 1024, /* 22 bits-per-sample */
322                 512, /* 23 bits-per-sample */
323                 256, /* 24 bits-per-sample */
324                 128, /* 25 bits-per-sample */
325                 64, /* 26 bits-per-sample */
326                 32, /* 27 bits-per-sample */
327                 16, /* 28 bits-per-sample */
328                 8, /* 29 bits-per-sample */
329                 4, /* 30 bits-per-sample */
330                 2, /* 31 bits-per-sample */
331                 1 /* 32 bits-per-sample */
332         };
333         static const FLAC__int64 hard_clip_factors_[33] = {
334                 0, /* 0 bits-per-sample (not supported) */
335                 0, /* 1 bits-per-sample (not supported) */
336                 0, /* 2 bits-per-sample (not supported) */
337                 0, /* 3 bits-per-sample (not supported) */
338                 -8, /* 4 bits-per-sample */
339                 -16, /* 5 bits-per-sample */
340                 -32, /* 6 bits-per-sample */
341                 -64, /* 7 bits-per-sample */
342                 -128, /* 8 bits-per-sample */
343                 -256, /* 9 bits-per-sample */
344                 -512, /* 10 bits-per-sample */
345                 -1024, /* 11 bits-per-sample */
346                 -2048, /* 12 bits-per-sample */
347                 -4096, /* 13 bits-per-sample */
348                 -8192, /* 14 bits-per-sample */
349                 -16384, /* 15 bits-per-sample */
350                 -32768, /* 16 bits-per-sample */
351                 -65536, /* 17 bits-per-sample */
352                 -131072, /* 18 bits-per-sample */
353                 -262144, /* 19 bits-per-sample */
354                 -524288, /* 20 bits-per-sample */
355                 -1048576, /* 21 bits-per-sample */
356                 -2097152, /* 22 bits-per-sample */
357                 -4194304, /* 23 bits-per-sample */
358                 -8388608, /* 24 bits-per-sample */
359                 -16777216, /* 25 bits-per-sample */
360                 -33554432, /* 26 bits-per-sample */
361                 -67108864, /* 27 bits-per-sample */
362                 -134217728, /* 28 bits-per-sample */
363                 -268435456, /* 29 bits-per-sample */
364                 -536870912, /* 30 bits-per-sample */
365                 -1073741824, /* 31 bits-per-sample */
366                 (FLAC__int64)(-1073741824) * 2 /* 32 bits-per-sample */
367         };
368         const FLAC__int32 conv_factor = conv_factors_[target_bps];
369         const FLAC__int64 hard_clip_factor = hard_clip_factors_[target_bps];
370         /*
371          * The integer input coming in has a varying range based on the
372          * source_bps.  We want to normalize it to [-1.0, 1.0) so instead
373          * of doing two multiplies on each sample, we just multiple
374          * 'scale' by 1/(2^(source_bps-1))
375          */
376         const double multi_scale = scale / (double)(1u << (source_bps-1));
377
378         FLAC__byte * const start = data_out;
379         unsigned i, channel;
380         const FLAC__int32 *input_;
381         double sample;
382         const unsigned bytes_per_sample = target_bps / 8;
383         const unsigned last_history_index = dither_context->LastHistoryIndex;
384         NoiseShaping noise_shaping = dither_context->ShapingType;
385         FLAC__int64 val64;
386         FLAC__int32 val32;
387         FLAC__int32 uval32;
388         const FLAC__uint32 twiggle = 1u << (target_bps - 1);
389
390         FLAC__ASSERT(channels > 0 && channels <= FLAC_SHARE__MAX_SUPPORTED_CHANNELS);
391         FLAC__ASSERT(source_bps >= 4);
392         FLAC__ASSERT(target_bps >= 4);
393         FLAC__ASSERT(source_bps <= 32);
394         FLAC__ASSERT(target_bps < 32);
395         FLAC__ASSERT((target_bps & 7) == 0);
396
397         for(channel = 0; channel < channels; channel++) {
398                 const unsigned incr = bytes_per_sample * channels;
399                 data_out = start + bytes_per_sample * channel;
400                 input_ = input[channel];
401                 for(i = 0; i < wide_samples; i++, data_out += incr) {
402                         sample = (double)input_[i] * multi_scale;
403
404                         if(hard_limit) {
405                                 /* hard 6dB limiting */
406                                 if(sample < -0.5)
407                                         sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
408                                 else if(sample > 0.5)
409                                         sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
410                         }
411                         sample *= 2147483647.;
412
413                         val64 = dither_output_(dither_context, do_dithering, noise_shaping, (i + last_history_index) % 32, sample, channel) / conv_factor;
414
415                         val32 = (FLAC__int32)val64;
416                         if(val64 >= -hard_clip_factor)
417                                 val32 = (FLAC__int32)(-(hard_clip_factor+1));
418                         else if(val64 < hard_clip_factor)
419                                 val32 = (FLAC__int32)hard_clip_factor;
420
421                         uval32 = (FLAC__uint32)val32;
422                         if (unsigned_data_out)
423                                 uval32 ^= twiggle;
424
425                         if (little_endian_data_out) {
426                                 switch(target_bps) {
427                                         case 24:
428                                                 data_out[2] = (FLAC__byte)(uval32 >> 16);
429                                                 /* fall through */
430                                         case 16:
431                                                 data_out[1] = (FLAC__byte)(uval32 >> 8);
432                                                 /* fall through */
433                                         case 8:
434                                                 data_out[0] = (FLAC__byte)uval32;
435                                                 break;
436                                 }
437                         }
438                         else {
439                                 switch(target_bps) {
440                                         case 24:
441                                                 data_out[0] = (FLAC__byte)(uval32 >> 16);
442                                                 data_out[1] = (FLAC__byte)(uval32 >> 8);
443                                                 data_out[2] = (FLAC__byte)uval32;
444                                                 break;
445                                         case 16:
446                                                 data_out[0] = (FLAC__byte)(uval32 >> 8);
447                                                 data_out[1] = (FLAC__byte)uval32;
448                                                 break;
449                                         case 8:
450                                                 data_out[0] = (FLAC__byte)uval32;
451                                                 break;
452                                 }
453                         }
454                 }
455         }
456         dither_context->LastHistoryIndex = (last_history_index + wide_samples) % 32;
457
458         return wide_samples * channels * (target_bps/8);
459 }