Mercurial > audlegacy-plugins

comparison src/flac/replaygain_synthesis.c @ 12:3da1b8942b8b trunk

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[svn] - remove src/Input src/Output src/Effect src/General src/Visualization src/Container
author nenolod
date Mon, 18 Sep 2006 03:14:20 -0700
parents src/Input/flac/replaygain_synthesis.c@13389e613d67
children
comparison
equal deleted inserted replaced
11:cff1d04026ae 12:3da1b8942b8b
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 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 "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 }