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