[pcsx_rearmed.git] / deps / flac-1.3.2 / src / share / replaygain_synthesis / replaygain_synthesis.c

/* replaygain_synthesis - Routines for applying ReplayGain to a signal
 * Copyright (C) 2002-2009  Josh Coalson
 * Copyright (C) 2011-2016  Xiph.Org Foundation
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
/*
 * This is an aggregation of pieces of code from John Edwards' WaveGain
 * program.  Mostly cosmetic changes were made; otherwise, the dithering
 * code is almost untouched and the gain processing was converted from
 * processing a whole file to processing chunks of samples.
 *
 * The original copyright notices for WaveGain's dither.c and wavegain.c
 * appear below:
 */
/*
 * (c) 2002 John Edwards
 * mostly lifted from work by Frank Klemm
 * random functions for dithering.
 */
/*
 * Copyright (C) 2002 John Edwards
 * Additional code by Magnus Holmgren and Gian-Carlo Pascutto
 */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

#include <string.h> /* for memset() */
#include <math.h>
#include "share/replaygain_synthesis.h"
#include "FLAC/assert.h"

#define FLAC__I64L(x) x##LL


/*
 * the following is based on parts of dither.c
 */


/*
 *  This is a simple random number generator with good quality for audio purposes.
 *  It consists of two polycounters with opposite rotation direction and different
 *  periods. The periods are coprime, so the total period is the product of both.
 *
 *     -------------------------------------------------------------------------------------------------
 * +-> |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|
 * |   -------------------------------------------------------------------------------------------------
 * |                                                                          |  |  |  |     |        |
 * |                                                                          +--+--+--+-XOR-+--------+
 * |                                                                                      |
 * +--------------------------------------------------------------------------------------+
 *
 *     -------------------------------------------------------------------------------------------------
 *     |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| <-+
 *     -------------------------------------------------------------------------------------------------   |
 *       |  |           |  |                                                                               |
 *       +--+----XOR----+--+                                                                               |
 *                |                                                                                        |
 *                +----------------------------------------------------------------------------------------+
 *
 *
 *  The first has an period of 3*5*17*257*65537, the second of 7*47*73*178481,
 *  which gives a period of 18.410.713.077.675.721.215. The result is the
 *  XORed values of both generators.
 */

static unsigned int random_int_(void)
{
	static const unsigned char parity_[256] = {
		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,
		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,
		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,
		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,
		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,
		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,
		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,
		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
	};
	static unsigned int r1_ = 1;
	static unsigned int r2_ = 1;

	unsigned int t1, t2, t3, t4;

	/* Parity calculation is done via table lookup, this is also available
	 * on CPUs without parity, can be implemented in C and avoid unpredictable
	 * jumps and slow rotate through the carry flag operations.
	 */
	t3   = t1 = r1_;    t4   = t2 = r2_;
	t1  &= 0xF5;        t2 >>= 25;
	t1   = parity_[t1]; t2  &= 0x63;
	t1 <<= 31;          t2   = parity_[t2];

	return (r1_ = (t3 >> 1) | t1 ) ^ (r2_ = (t4 + t4) | t2 );
}

/* gives a equal distributed random number */
/* between -2^31*mult and +2^31*mult */
static double random_equi_(double mult)
{
	return mult * (int) random_int_();
}

/* gives a triangular distributed random number */
/* between -2^32*mult and +2^32*mult */
static double random_triangular_(double mult)
{
	return mult * ( (double) (int) random_int_() + (double) (int) random_int_() );
}


static const float  F44_0 [16 + 32] = {
	(float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
	(float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,

	(float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
	(float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,

	(float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0,
	(float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0, (float)0
};


static const float  F44_1 [16 + 32] = {  /* SNR(w) = 4.843163 dB, SNR = -3.192134 dB */
	(float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
	(float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
	(float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
	(float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,

	(float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
	(float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
	(float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
	(float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,

	(float) 0.85018292704024355931, (float) 0.29089597350995344721, (float)-0.05021866022121039450, (float)-0.23545456294599161833,
	(float)-0.58362726442227032096, (float)-0.67038978965193036429, (float)-0.38566861572833459221, (float)-0.15218663390367969967,
	(float)-0.02577543084864530676, (float) 0.14119295297688728127, (float) 0.22398848581628781612, (float) 0.15401727203382084116,
	(float) 0.05216161232906000929, (float)-0.00282237820999675451, (float)-0.03042794608323867363, (float)-0.03109780942998826024,
};


static const float  F44_2 [16 + 32] = {  /* SNR(w) = 10.060213 dB, SNR = -12.766730 dB */
	(float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
	(float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
	(float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
	(float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,

	(float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
	(float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
	(float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
	(float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,

	(float) 1.78827593892108555290, (float) 0.95508210637394326553, (float)-0.18447626783899924429, (float)-0.44198126506275016437,
	(float)-0.88404052492547413497, (float)-1.42218907262407452967, (float)-1.02037566838362314995, (float)-0.34861755756425577264,
	(float)-0.11490230170431934434, (float) 0.12498899339968611803, (float) 0.38065885268563131927, (float) 0.31883491321310506562,
	(float) 0.10486838686563442765, (float)-0.03105361685110374845, (float)-0.06450524884075370758, (float)-0.02939198261121969816,
};


static const float  F44_3 [16 + 32] = {  /* SNR(w) = 15.382598 dB, SNR = -29.402334 dB */
	(float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
	(float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
	(float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
	(float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099,

	(float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
	(float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
	(float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
	(float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099,

	(float) 2.89072132015058161445, (float) 2.68932810943698754106, (float) 0.21083359339410251227, (float)-0.98385073324997617515,
	(float)-1.11047823227097316719, (float)-2.18954076314139673147, (float)-2.36498032881953056225, (float)-0.95484132880101140785,
	(float)-0.23924057925542965158, (float)-0.13865235703915925642, (float) 0.43587843191057992846, (float) 0.65903257226026665927,
	(float) 0.24361815372443152787, (float)-0.00235974960154720097, (float) 0.01844166574603346289, (float) 0.01722945988740875099
};


static double scalar16_(const float* x, const float* y)
{
	return
		x[ 0]*y[ 0] + x[ 1]*y[ 1] + x[ 2]*y[ 2] + x[ 3]*y[ 3] +
		x[ 4]*y[ 4] + x[ 5]*y[ 5] + x[ 6]*y[ 6] + x[ 7]*y[ 7] +
		x[ 8]*y[ 8] + x[ 9]*y[ 9] + x[10]*y[10] + x[11]*y[11] +
		x[12]*y[12] + x[13]*y[13] + x[14]*y[14] + x[15]*y[15];
}


void FLAC__replaygain_synthesis__init_dither_context(DitherContext *d, int bits, int shapingtype)
{
	static unsigned char default_dither [] = { 92, 92, 88, 84, 81, 78, 74, 67,  0,  0 };
	static const float*               F [] = { F44_0, F44_1, F44_2, F44_3 };

	int indx;

	if (shapingtype < 0) shapingtype = 0;
	if (shapingtype > 3) shapingtype = 3;
	d->ShapingType = (NoiseShaping)shapingtype;
	indx = bits - 11 - shapingtype;
	if (indx < 0) indx = 0;
	if (indx > 9) indx = 9;

	memset ( d->ErrorHistory , 0, sizeof (d->ErrorHistory ) );
	memset ( d->DitherHistory, 0, sizeof (d->DitherHistory) );

	d->FilterCoeff = F [shapingtype];
	d->Mask   = ((FLAC__uint64)-1) << (32 - bits);
	d->Add    = 0.5     * ((1L << (32 - bits)) - 1);
	d->Dither = 0.01f*default_dither[indx] / (((FLAC__int64)1) << bits);
	d->LastHistoryIndex = 0;
}

/*
 * the following is based on parts of wavegain.c
 */

static FLAC__int64 dither_output_(DitherContext *d, FLAC__bool do_dithering, int shapingtype, int i, double Sum, int k)
{
	union {
		double d;
		FLAC__int64 i;
	} doubletmp;
	double Sum2;
	FLAC__int64 val;

#define ROUND64(x)   ( doubletmp.d = (x) + d->Add + (FLAC__int64)FLAC__I64L(0x001FFFFD80000000), doubletmp.i - (FLAC__int64)FLAC__I64L(0x433FFFFD80000000) )

	if(do_dithering) {
		if(shapingtype == 0) {
			double  tmp = random_equi_(d->Dither);
			Sum2 = tmp - d->LastRandomNumber [k];
			d->LastRandomNumber [k] = (int)tmp;
			Sum2 = Sum += Sum2;
			val = ROUND64(Sum2) & d->Mask;
		}
		else {
			Sum2 = random_triangular_(d->Dither) - scalar16_(d->DitherHistory[k], d->FilterCoeff + i);
			Sum += d->DitherHistory [k] [(-1-i)&15] = (float)Sum2;
			Sum2 = Sum + scalar16_(d->ErrorHistory [k], d->FilterCoeff + i);
			val = ROUND64(Sum2) & d->Mask;
			d->ErrorHistory [k] [(-1-i)&15] = (float)(Sum - val);
		}
		return val;
	}
	else
		return ROUND64(Sum);

#undef ROUND64
}

#if 0
	float        peak = 0.f,
	             new_peak,
	             factor_clip
	double       scale,
	             dB;

	...

	peak is in the range -32768.0 .. 32767.0

	/* calculate factors for ReplayGain and ClippingPrevention */
	*track_gain = GetTitleGain() + settings->man_gain;
	scale = (float) pow(10., *track_gain * 0.05);
	if(settings->clip_prev) {
		factor_clip  = (float) (32767./( peak + 1));
		if(scale < factor_clip)
			factor_clip = 1.f;
		else
			factor_clip /= scale;
		scale *= factor_clip;
	}
	new_peak = (float) peak * scale;

	dB = 20. * log10(scale);
	*track_gain = (float) dB;

 	const double scale = pow(10., (double)gain * 0.05);
#endif


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)
{
	static const FLAC__int64 hard_clip_factors_[33] = {
		0, /* 0 bits-per-sample (not supported) */
		0, /* 1 bits-per-sample (not supported) */
		0, /* 2 bits-per-sample (not supported) */
		0, /* 3 bits-per-sample (not supported) */
		-8, /* 4 bits-per-sample */
		-16, /* 5 bits-per-sample */
		-32, /* 6 bits-per-sample */
		-64, /* 7 bits-per-sample */
		-128, /* 8 bits-per-sample */
		-256, /* 9 bits-per-sample */
		-512, /* 10 bits-per-sample */
		-1024, /* 11 bits-per-sample */
		-2048, /* 12 bits-per-sample */
		-4096, /* 13 bits-per-sample */
		-8192, /* 14 bits-per-sample */
		-16384, /* 15 bits-per-sample */
		-32768, /* 16 bits-per-sample */
		-65536, /* 17 bits-per-sample */
		-131072, /* 18 bits-per-sample */
		-262144, /* 19 bits-per-sample */
		-524288, /* 20 bits-per-sample */
		-1048576, /* 21 bits-per-sample */
		-2097152, /* 22 bits-per-sample */
		-4194304, /* 23 bits-per-sample */
		-8388608, /* 24 bits-per-sample */
		-16777216, /* 25 bits-per-sample */
		-33554432, /* 26 bits-per-sample */
		-67108864, /* 27 bits-per-sample */
		-134217728, /* 28 bits-per-sample */
		-268435456, /* 29 bits-per-sample */
		-536870912, /* 30 bits-per-sample */
		-1073741824, /* 31 bits-per-sample */
		(FLAC__int64)(-1073741824) * 2 /* 32 bits-per-sample */
	};
	const FLAC__int32 conv_shift = 32 - target_bps;
	const FLAC__int64 hard_clip_factor = hard_clip_factors_[target_bps];
	/*
	 * The integer input coming in has a varying range based on the
	 * source_bps.  We want to normalize it to [-1.0, 1.0) so instead
	 * of doing two multiplies on each sample, we just multiple
	 * 'scale' by 1/(2^(source_bps-1))
	 */
	const double multi_scale = scale / (double)(1u << (source_bps-1));

	FLAC__byte * const start = data_out;
	unsigned i, channel;
	const FLAC__int32 *input_;
	double sample;
	const unsigned bytes_per_sample = target_bps / 8;
	const unsigned last_history_index = dither_context->LastHistoryIndex;
	NoiseShaping noise_shaping = dither_context->ShapingType;
	FLAC__int64 val64;
	FLAC__int32 val32;
	FLAC__int32 uval32;
	const FLAC__uint32 twiggle = 1u << (target_bps - 1);

	FLAC__ASSERT(channels > 0 && channels <= FLAC_SHARE__MAX_SUPPORTED_CHANNELS);
	FLAC__ASSERT(source_bps >= 4);
	FLAC__ASSERT(target_bps >= 4);
	FLAC__ASSERT(source_bps <= 32);
	FLAC__ASSERT(target_bps < 32);
	FLAC__ASSERT((target_bps & 7) == 0);

	for(channel = 0; channel < channels; channel++) {
		const unsigned incr = bytes_per_sample * channels;
		data_out = start + bytes_per_sample * channel;
		input_ = input[channel];
		for(i = 0; i < wide_samples; i++, data_out += incr) {
			sample = (double)input_[i] * multi_scale;

			if(hard_limit) {
				/* hard 6dB limiting */
				if(sample < -0.5)
					sample = tanh((sample + 0.5) / (1-0.5)) * (1-0.5) - 0.5;
				else if(sample > 0.5)
					sample = tanh((sample - 0.5) / (1-0.5)) * (1-0.5) + 0.5;
			}
			sample *= 2147483647.;

			val64 = dither_output_(dither_context, do_dithering, noise_shaping, (i + last_history_index) % 32, sample, channel) >> conv_shift;

			val32 = (FLAC__int32)val64;
			if(val64 >= -hard_clip_factor)
				val32 = (FLAC__int32)(-(hard_clip_factor+1));
			else if(val64 < hard_clip_factor)
				val32 = (FLAC__int32)hard_clip_factor;

			uval32 = (FLAC__uint32)val32;
			if (unsigned_data_out)
				uval32 ^= twiggle;

			if (little_endian_data_out) {
				switch(target_bps) {
					case 24:
						data_out[2] = (FLAC__byte)(uval32 >> 16);
						/* fall through */
					case 16:
						data_out[1] = (FLAC__byte)(uval32 >> 8);
						/* fall through */
					case 8:
						data_out[0] = (FLAC__byte)uval32;
						break;
				}
			}
			else {
				switch(target_bps) {
					case 24:
						data_out[0] = (FLAC__byte)(uval32 >> 16);
						data_out[1] = (FLAC__byte)(uval32 >> 8);
						data_out[2] = (FLAC__byte)uval32;
						break;
					case 16:
						data_out[0] = (FLAC__byte)(uval32 >> 8);
						data_out[1] = (FLAC__byte)uval32;
						break;
					case 8:
						data_out[0] = (FLAC__byte)uval32;
						break;
				}
			}
		}
	}
	dither_context->LastHistoryIndex = (last_history_index + wide_samples) % 32;

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