History

Further fix UB in PNS (#150)

And simplify code. Given the nature of noise we do not have to match noise energy exactly. Matching the expected energy seems to be good enough, and have much cleaner code. Now we do not need FP division and square root.

Eugene Kliuchnikov committed 3 years ago

Commit 75a85f81 parent d062deb2 files · +47 −61Browse files

Unified · Split

2 files
  1. libfaad/decoder.c +3 −2
  2. libfaad/pns.c +44 −59
libfaad/decoder.c+3−2
@@ -148,8 +148,9 @@ NeAACDecHandle NeAACDecOpen(void) · 147 unmodified lines
148148 hDecoder->frame = 0;
149149 hDecoder->sample_buffer = NULL;
150150
151− hDecoder->__r1 = 1;
152− hDecoder->__r2 = 1;
151+ // Same as (1, 1) after 1024 iterations; otherwise first values does not look random at all.
152+ hDecoder->__r1 = 0x2bb431ea;
153+ hDecoder->__r2 = 0x206155b7;
153154
154155 for (i = 0; i < MAX_CHANNELS; i++)
155156 {
libfaad/pns.c+44−59
@@ -42,36 +42,6 @@ static void gen_rand_vector(real_t *spec, int16_t scale_factor, uint16_t size, · 41 unmodified lines
4242
4343 #ifdef FIXED_POINT
4444
45−#define DIV(A, B) (((int64_t)A << REAL_BITS)/B)
46−
47−#define step(shift) \
48− if ((0x40000000l >> shift) + root <= value) \
49− { \
50− value -= (0x40000000l >> shift) + root; \
51− root = (root >> 1) | (0x40000000l >> shift); \
52− } else { \
53− root = root >> 1; \
54− }
55−
56−/* fixed point square root approximation */
57−/* !!!! ONLY WORKS FOR EVEN %REAL_BITS% !!!! */
58−real_t fp_sqrt(real_t value)
59−{
60− real_t root = 0;
61−
62− step( 0); step( 2); step( 4); step( 6);
63− step( 8); step(10); step(12); step(14);
64− step(16); step(18); step(20); step(22);
65− step(24); step(26); step(28); step(30);
66−
67− if (root < value)
68− ++root;
69−
70− root <<= (REAL_BITS/2);
71−
72− return root;
73−}
74−
7545 static real_t const pow2_table[] =
7646 {
7747 COEF_CONST(1.0),
@@ -79,6 +49,27 @@ static real_t const pow2_table[] = · 1 unmodified line
7949 COEF_CONST(1.41421356237310),
8050 COEF_CONST(1.68179283050743)
8151 };
52+
53+// mean_energy_table[x] == sqrt(3 / x)
54+static real_t const mean_energy_table[] =
55+{
56+ COEF_CONST(0.0), // should not happen
57+ COEF_CONST(1.7320508075688772),
58+ COEF_CONST(1.224744871391589),
59+ COEF_CONST(1.0), // sqrt(3/3)
60+ COEF_CONST(0.8660254037844386),
61+ COEF_CONST(0.7745966692414834),
62+ COEF_CONST(0.7071067811865476),
63+ COEF_CONST(0.6546536707079771),
64+ COEF_CONST(0.6123724356957945),
65+ COEF_CONST(0.5773502691896257),
66+ COEF_CONST(0.5477225575051661),
67+ COEF_CONST(0.5222329678670935),
68+ COEF_CONST(0.5), // sqrt(3/12)
69+ COEF_CONST(0.4803844614152614),
70+ COEF_CONST(0.4629100498862757),
71+ COEF_CONST(0.4472135954999579),
72+};
8273 #endif
8374
8475 /* The function gen_rand_vector(addr, size) generates a vector of length
@@ -114,49 +105,43 @@ static INLINE void gen_rand_vector(real_t *spec, int16_t scale_factor, uint16_t · 29 unmodified lines
114105 }
115106 #else
116107 uint16_t i;
117− real_t energy = 0, scale;
108+ real_t scale;
118109 int32_t exp, frac;
110+ int32_t idx, mask;
119111
120112 /* IMDCT pre-scaling */
121113 scale_factor -= 4 * sub;
122114
123− scale_factor = min(max(scale_factor, -60), 60);
115+ // 52 stands for 2**13 == 8192 factor; larger factor causes overflows later (in cfft).
116+ scale_factor = min(max(scale_factor, -(REAL_BITS * 4)), 52);
124117
125118 exp = scale_factor >> 2;
126119 frac = scale_factor & 3;
127120
128− for (i = 0; i < size; i++)
129− {
130− /* this can be replaced by a 16 bit random generator!!!! */
131− real_t tmp = (int32_t)ne_rng(__r1, __r2);
132− if (tmp < 0)
133− tmp = -(tmp & ((1<<(REAL_BITS-1))-1));
134− else
135− tmp = (tmp & ((1<<(REAL_BITS-1))-1));
136−
137− energy += MUL_R(tmp,tmp);
121+ /* 29 <= REAL_BITS + exp <= 0 */
122+ mask = (1 << (REAL_BITS + exp)) - 1;
138123
139− spec[i] = tmp;
124+ idx = size;
125+ scale = COEF_CONST(1);
126+ // At most 2 iterations.
127+ while (idx >= 16)
128+ {
129+ idx >>= 2;
130+ scale >>= 1;
140131 }
141− /* energy ~= size / 3 */
132+ scale = MUL_C(scale, mean_energy_table[idx]);
133+ if (frac)
134+ scale = MUL_C(scale, pow2_table[frac]);
135+ // scale is less than 4.0 now.
142136
143− energy = fp_sqrt(energy);
144− if (energy > 0)
137+ for (i = 0; i < size; i++)
145138 {
146− scale = DIV(REAL_CONST(1),energy);
147−
148− if (exp < 0)
149− scale >>= -exp;
139+ real_t tmp = (int32_t)ne_rng(__r1, __r2);
140+ if (tmp < 0)
141+ tmp = -(tmp & mask);
150142 else
151− scale <<= exp;
152−
153− if (frac)
154− scale = MUL_C(scale, pow2_table[frac]);
155−
156− for (i = 0; i < size; i++)
157− {
158− spec[i] = MUL_R(spec[i], scale);
159− }
143+ tmp = (tmp & mask);
144+ spec[i] = MUL_C(tmp, scale);
160145 }
161146 #endif
162147 }

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t