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