-
Notifications
You must be signed in to change notification settings - Fork 5
/
dither.c
184 lines (146 loc) · 10.3 KB
/
dither.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
/* This program is licensed under the GNU Library General Public License, version 2,
* a copy of which is included with this program (with filename LICENSE.LGPL).
*
* (c) 2002 John Edwards
* mostly lifted from work by Frank Klemm
* random functions for dithering.
*
* last modified: $ID:$
*/
#include "dither.h"
#include <string.h>
static const unsigned char Parity [256] = { // parity
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;
/*
* 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.
*/
unsigned int
random_int ( void )
{
unsigned int t1, t2, t3, t4;
t3 = t1 = __r1; t4 = t2 = __r2; // Parity calculation is done via table lookup, this is also available
t1 &= 0xF5; t2 >>= 25; // on CPUs without parity, can be implemented in C and avoid unpredictable
t1 = Parity [t1]; t2 &= 0x63; // jumps and slow rotate through the carry flag operations.
t1 <<= 31; t2 = Parity [t2];
return (__r1 = (t3 >> 1) | t1 ) ^ (__r2 = (t4 + t4) | t2 );
}
double
Random_Equi ( double mult ) // gives a equal distributed random number
{ // between -2^31*mult and +2^31*mult
return mult * (int) random_int ();
}
double
Random_Triangular ( double mult ) // gives a triangular distributed random number
{ // between -2^32*mult and +2^32*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
};
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
Init_Dither ( 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 index;
if (shapingtype < 0) shapingtype = 0;
if (shapingtype > 3) shapingtype = 3;
index = bits - 11 - shapingtype;
if (index < 0) index = 0;
if (index > 9) index = 9;
memset ( Dither.ErrorHistory , 0, sizeof (Dither.ErrorHistory ) );
memset ( Dither.DitherHistory, 0, sizeof (Dither.DitherHistory) );
Dither.FilterCoeff = F [shapingtype];
Dither.Mask = ((Uint64_t)-1) << (32 - bits);
Dither.Add = 0.5 * ((1L << (32 - bits)) - 1);
Dither.Dither = 0.01*default_dither[index] / (((Int64_t)1) << bits);
}