1 /* $Id: float32.c,v 1.11 2007/09/19 11:49:46 toad32767 Exp $ */
2
3 /*-
4 * Copyright (c) 2005, 2006, 2007 Marco Trillo
5 *
6 * Permission is hereby granted, free of charge, to any
7 * person obtaining a copy of this software and associated
8 * documentation files (the "Software"), to deal in the
9 * Software without restriction, including without limitation
10 * the rights to use, copy, modify, merge, publish,
11 * distribute, sublicense, and/or sell copies of the
12 * Software, and to permit persons to whom the Software is
13 * furnished to do so, subject to the following conditions:
14 *
15 * The above copyright notice and this permission notice
16 * shall be included in all copies or substantial portions of
17 * the Software.
18 *
19 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY
20 * KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
21 * WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
22 * PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
23 * OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
24 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
25 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
26 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
27 */
28
29 #define LIBAIFF 1
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <math.h>
33 #include <libaiff/libaiff.h>
34 #include <libaiff/endian.h>
35 #include "private.h"
36
37 /*
38 * private flags for this module
39 */
40 #define F_IEEE754_CHECKED (1<<27)
41 #define F_IEEE754_NATIVE (1<<28)
42
43
44 /*
45 * IEEE-754 32-bit single-precision floating point
46 *
47 * This is a conversor to linear 32-bit integer PCM
48 * using only integer arithmetic.
49 *
50 * Some IEEE-754 documentation and references on the WWW:
51 * -------------------------------------------------------------
52 * http://stevehollasch.com/cgindex/coding/ieeefloat.html
53 * http://www.answers.com/topic/ieee-floating-point-standard
54 * http://www.cs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF
55 * http://www.psc.edu/general/software/packages/ieee/ieee.html
56 * http://www.duke.edu/~twf/cps104/floating.html
57 * -------------------------------------------------------------
58 */
59
60 static int32_t
float32dec(uint32_t in)61 float32dec(uint32_t in)
62 {
63 int sgn, exp;
64 uint32_t mantissa;
65
66 if (in == 0 || in == 0x80000000)
67 return (0);
68
69 sgn = (int) (in >> 31);
70 exp = (int) ((in >> 23) & 0xFF);
71 mantissa = in & 0x7FFFFF;
72
73 if (exp == 255) {
74 /*
75 * Infinity and NaN.
76 * XXX: what we should do with these?
77 */
78 if (mantissa == 0) /* infinity */
79 return (int32_t) (0x7FFFFFFF);
80 else /* NaN */
81 return (0);
82 } else if (exp == 0) {
83 exp = -126; /* denormalized number */
84 } else {
85 /* Normalized numbers */
86 mantissa |= 0x800000; /* hidden bit */
87 exp -= 127; /* unbias exponent */
88 }
89
90 /*
91 * Quantization to linear PCM 31 bits.
92 *
93 * Multiply the mantissa by 2^(-23+exp) to get the floating
94 * point value, and then multiply by 2^31 to quantize
95 * the floating point to 31 bits. So:
96 *
97 * 2^(-23+exp) * 2^31 = 2^(8 + exp)
98 */
99 exp += 8;
100 if (exp < 0) {
101 exp = -exp;
102 mantissa >>= exp;
103 } else {
104 mantissa <<= exp;
105 }
106
107 if (sgn) {
108 int32_t val;
109
110 if (mantissa == 0x80000000)
111 return (-2147483647 - 1); /* -(2^31) */
112 else {
113 val = (int32_t) mantissa;
114 return (-val);
115 }
116 } else {
117 /*
118 * if mantissa is 0x80000000 (for a 1.0 float),
119 * we should return +2147483648 (2^31), but in
120 * two's complement arithmetic the max. positive value
121 * is +2147483647, so clip the result.
122 */
123 if (mantissa == 0x80000000)
124 return (2147483647);
125 else
126 return ((int32_t) mantissa);
127 }
128 }
129
130 static void
float32_decode(void * dst,void * src,int n)131 float32_decode(void *dst, void *src, int n)
132 {
133 uint32_t *srcSamples = (uint32_t *) src;
134 int32_t *dstSamples = (int32_t *) dst;
135
136 while (n-- > 0)
137 dstSamples[n] = float32dec(srcSamples[n]);
138 }
139
140 static void
float32_swap_samples(void * buffer,int n)141 float32_swap_samples(void *buffer, int n)
142 {
143 uint32_t *streams = (uint32_t *) buffer;
144
145 while (n-- > 0)
146 streams[n] = ARRANGE_ENDIAN_32(streams[n]);
147 }
148
149 static size_t
float32_read_lpcm(AIFF_Ref r,void * buffer,size_t len)150 float32_read_lpcm(AIFF_Ref r, void *buffer, size_t len)
151 {
152 int n;
153 uint32_t clen;
154 size_t slen;
155 size_t bytes_in;
156 size_t bytesToRead;
157
158 n = (int) len;
159 /* 'n' should be divisible by 4 (32 bits) */
160 while (n >= 0 && ((n & 3) != 0)) {
161 --n;
162 --len;
163 }
164 n >>= 2;
165
166 slen = (size_t) (r->soundLen) - (size_t) (r->pos);
167 bytesToRead = MIN(len, slen);
168
169 if (bytesToRead == 0)
170 return 0;
171 if (r->buffer2 == NULL || r->buflen2 < bytesToRead) {
172 if (r->buffer2 != NULL)
173 free(r->buffer2);
174 r->buffer2 = malloc(bytesToRead);
175 if (r->buffer2 == NULL) {
176 r->buflen2 = 0;
177 return 0;
178 }
179 r->buflen2 = bytesToRead;
180 }
181
182 bytes_in = fread(r->buffer2, 1, bytesToRead, r->fd);
183 if (bytes_in > 0)
184 clen = (uint32_t) bytes_in;
185 else
186 clen = 0;
187 r->pos += clen;
188
189 if (r->flags & LPCM_NEED_SWAP)
190 float32_swap_samples(r->buffer2, n);
191 float32_decode(buffer, r->buffer2, n);
192
193 return bytes_in;
194 }
195
196 static int
float32_seek(AIFF_Ref r,uint64_t pos)197 float32_seek(AIFF_Ref r, uint64_t pos)
198 {
199 long of;
200 uint32_t b;
201
202 b = (uint32_t) pos * r->nChannels * 4;
203 if (b >= r->soundLen)
204 return 0;
205 of = (long) b;
206
207 if (fseek(r->fd, of, SEEK_CUR) < 0) {
208 return -1;
209 }
210 r->pos = b;
211 return 1;
212 }
213
214 /*
215 * Infinite & NAN values
216 * for non-IEEE systems
217 */
218 #ifndef HUGE_VAL
219 # ifdef HUGE
220 # define INFINITE_VALUE HUGE
221 # define NAN_VALUE HUGE
222 # endif
223 #else
224 # define INFINITE_VALUE HUGE_VAL
225 # define NAN_VALUE HUGE_VAL
226 #endif
227
228 /*
229 * Write IEEE Single Precision Numbers
230 */
231 static uint32_t
ieee754_write_single(float in)232 ieee754_write_single(float in)
233 {
234 uint32_t sgn ; /* sign */
235 double fraction ;
236 uint32_t mantissa ;
237 int exp ;
238 uint32_t bitstream ;
239
240 /* Zero`s can be processed quickly */
241 if (in == 0.0)
242 {
243 return 0;
244 }
245
246 /*
247 * For negative values we set the sign to 1
248 * and use the absolute value
249 */
250 if (in < 0.0)
251 {
252 in = (float)fabs(in);
253 sgn = 1;
254 }
255 else
256 {
257 sgn = 0;
258 }
259
260 /*
261 * Grab the exponent & mantissa.
262 * This function will return
263 * mantissas of type 0.1F .
264 *
265 * These mantissas are valid only
266 * for 'denormalised numbers'.
267 */
268 fraction = ldexp( frexp( in, &exp ), 24 ) ; /* extract 24-bits of
269 mantissa.
270 (we have only 23-bits
271 of precision, but
272 we multiply it by 2
273 to use in normalised
274 numbers)
275 */
276 mantissa = (uint32_t)(floor( fraction ));
277
278 /*
279 * Check for special numbers (NaN or infinity)
280 * and for out-of-range exponents (> 128)
281 */
282 if (exp == 0 || exp > 128)
283 {
284 if (exp > 128)
285 mantissa = 0 ; /* infinity have a mantissa of 0 */
286 else
287 mantissa = 0x400000 ; /* non-zero */
288
289 exp = 255 ;
290 goto done ;
291 }
292
293 /*
294 * Using '0.F' mantissas, if the exponent is -126 or lesser,
295 * we use <denormalised numbers>
296 */
297 if (exp <= -126)
298 {
299 /* Mantissa is multiplied by 2,
300 * so divide it.
301 *
302 * Convert also exponents lesser than -126
303 * to -126 (by dividing the mantissa).
304 */
305 mantissa >>= (1 - exp - 126) ;
306 exp = 0 ;
307 goto done ;
308 }
309
310 /* Any other number are a stored
311 * as <normalised number>.
312 *
313 * But for this numbers
314 * we need a mantissa of '1.F',
315 * and we have '0.1F',
316 * so decrement the 2's exponent,
317 * and multiply by 2 the mantissa
318 * (which is already done)
319 */
320
321 mantissa &= 0x7FFFFF ; /* hidden bit */
322 exp += 127-1 ; /* biased exponent */
323
324 done:
325 /* Construct the bitstream */
326 bitstream = (uint32_t)exp ;
327 bitstream &= 0xFF ; /* assert that exp (biased) is in [0,255] */
328 bitstream <<= 23 ; /* exponent is on offset {1,8} */
329 sgn <<= 31 ; /* sign is on offset {0} */
330 bitstream |= sgn ;
331 bitstream |= mantissa ; /* mantissa is on offset {9,31} */
332
333 return bitstream ;
334 }
335
336 #define IEEE754_TEST_VALUE -1.12877
337
338 /*
339 * ieee754_native(): check if host supports single-precision IEEE-754 natively
340 */
341 static int
ieee754_native(void)342 ieee754_native(void)
343 {
344 float f = IEEE754_TEST_VALUE, *ptr;
345 uint32_t sf, *hw;
346
347 if (sizeof(float) != sizeof(uint32_t))
348 return (0);
349
350 sf = ieee754_write_single(f);
351 ptr = &f;
352 hw = (uint32_t *) ptr;
353
354 if (sf == *hw)
355 return (1);
356 else
357 return (0);
358
359 return (0); /* NOTREACHED */
360 }
361
362
363 /*
364 * Read IEEE Single Precision Numbers
365 */
366 static float
ieee754_read_single(uint32_t in)367 ieee754_read_single(uint32_t in)
368 {
369 uint32_t sgn ;
370 int exp ;
371 uint32_t mantissa ;
372 double fraction ;
373 float out ;
374
375 if( in == 0 )
376 return 0.0 ;
377 else if( in == 0x80000000 )
378 return -0.0 ;
379
380 sgn = in >> 31 ;
381 exp = (int)( (in >> 23) & 0xFF ) ;
382 mantissa = in & 0x7FFFFF ;
383
384 switch( exp )
385 {
386 case 255:
387 switch( mantissa )
388 {
389 case 0:
390 return ( sgn ? -INFINITE_VALUE : INFINITE_VALUE ) ;
391 default:
392 return ( sgn ? -NAN_VALUE : NAN_VALUE ) ;
393 }
394 case 0:
395 /* Denormalised numbers */
396 exp = -126 ;
397 break ;
398 default:
399 /* Normalised numbers */
400 mantissa |= 0x800000 ; /* hidden bit */
401 exp -= 127 ; /* unbias exponent */
402 }
403
404 fraction = (double)mantissa ;
405 out = (float)ldexp( fraction, -23+exp ) ;
406
407 return ( sgn ? -out : out ) ;
408 }
409
410
411 static int
float32_read_float32(AIFF_Ref r,float * buffer,int n)412 float32_read_float32(AIFF_Ref r, float *buffer, int n)
413 {
414 int nSamplesRead;
415 uint32_t clen;
416 size_t len, slen;
417 size_t bytes_in;
418 size_t bytesToRead;
419
420 len = (size_t) n << 2; /* n * 4 */
421 slen = (size_t) (r->soundLen) - (size_t) (r->pos);
422 bytesToRead = MIN(len, slen);
423 if (bytesToRead == 0)
424 return 0;
425
426 /*
427 * Check if this host supports 32-bit IEEE floats
428 * natively and take note about it to avoid doing
429 * the test each time a buffer is to be read...
430 */
431 if (!(r->flags & F_IEEE754_CHECKED))
432 {
433 r->flags |= F_IEEE754_CHECKED;
434 if (ieee754_native())
435 r->flags |= F_IEEE754_NATIVE;
436 }
437
438 if (r->flags & F_IEEE754_NATIVE)
439 {
440 bytes_in = fread((void *) buffer, 1, bytesToRead, r->fd);
441
442 if (bytes_in > 0)
443 {
444 nSamplesRead = (int) bytes_in >> 2;
445 if (r->flags & LPCM_NEED_SWAP)
446 float32_swap_samples((void *) buffer, nSamplesRead);
447 }
448 else
449 {
450 nSamplesRead = 0;
451 }
452 }
453 else
454 {
455 if (r->buffer2 == NULL || r->buflen2 < bytesToRead) {
456 if (r->buffer2 != NULL)
457 free(r->buffer2);
458 r->buffer2 = malloc(bytesToRead);
459 if (r->buffer2 == NULL) {
460 r->buflen2 = 0;
461 return 0;
462 }
463 r->buflen2 = bytesToRead;
464 }
465
466 bytes_in = fread(r->buffer2, 1, bytesToRead, r->fd);
467 if (bytes_in > 0)
468 {
469 uint32_t *dwords = (uint32_t *) (r->buffer2);
470
471 nSamplesRead = (int) bytes_in >> 2;
472 if (r->flags & LPCM_NEED_SWAP)
473 float32_swap_samples(dwords, nSamplesRead);
474
475 while (nSamplesRead-- > 0)
476 {
477 buffer[nSamplesRead] = ieee754_read_single(dwords[nSamplesRead]);
478 }
479 }
480 else
481 {
482 nSamplesRead = 0;
483 }
484 }
485
486 if (bytes_in > 0)
487 clen = (uint32_t) bytes_in;
488 else
489 clen = 0;
490 r->pos += clen;
491
492 return nSamplesRead;
493 }
494
495
496 struct decoder float32 = {
497 AUDIO_FORMAT_FL32,
498 NULL,
499 float32_read_lpcm,
500 float32_read_float32,
501 float32_seek,
502 NULL
503 };
504
505