Edinburgh Speech Tools 2.4-release
 
Loading...
Searching...
No Matches
EST_wave_utils.cc
1/*************************************************************************/
2/* */
3/* Centre for Speech Technology Research */
4/* University of Edinburgh, UK */
5/* Copyright (c) 1996 */
6/* All Rights Reserved. */
7/* */
8/* Permission is hereby granted, free of charge, to use and distribute */
9/* this software and its documentation without restriction, including */
10/* without limitation the rights to use, copy, modify, merge, publish, */
11/* distribute, sublicense, and/or sell copies of this work, and to */
12/* permit persons to whom this work is furnished to do so, subject to */
13/* the following conditions: */
14/* 1. The code must retain the above copyright notice, this list of */
15/* conditions and the following disclaimer. */
16/* 2. Any modifications must be clearly marked as such. */
17/* 3. Original authors' names are not deleted. */
18/* 4. The authors' names are not used to endorse or promote products */
19/* derived from this software without specific prior written */
20/* permission. */
21/* */
22/* THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK */
23/* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
24/* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
25/* SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE */
26/* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
27/* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
28/* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
29/* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
30/* THIS SOFTWARE. */
31/* */
32/*************************************************************************/
33/* Author : Alan Black and Paul Taylor */
34/* Date : June 1996 */
35/*-----------------------------------------------------------------------*/
36/* Various low-level waveform conversion routines and file format */
37/* independent i/o functions */
38/* */
39/* Acknowledgements */
40/* ulaw conversion code provided by */
41/* Craig Reese: IDA/Supercomputing Research Center */
42/* IEEE extended conversion */
43/* Apple Computer, Inc. */
44/* */
45/*=======================================================================*/
46#include <cstdio>
47#include <cstdlib>
48#include "EST_unix.h"
49#include <cstring>
50#include <cmath>
51#include "EST_wave_utils.h"
52#include "EST_wave_aux.h"
53#include "EST_error.h"
54
55static short st_ulaw_to_short(unsigned char ulawbyte);
56static short st_alaw_to_short(unsigned char alawbyte);
57static unsigned char st_short_to_ulaw(short sample);
58static unsigned char st_short_to_alaw(short sample);
59
60/*
61 * This table is
62 * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
63 * Universitaet Berlin. See the accompanying file "COPYRIGHT" for
64 * details. THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
65 * take from speak_freely-7.2/gsm/src/toast_alaw.c
66 */
67static unsigned short a2s[] = {
68 5120,60160, 320,65200,20480,44032, 1280,64192,
69 2560,62848, 64,65456,10240,54784, 640,64864,
70 7168,58112, 448,65072,28672,35840, 1792,63680,
71 3584,61824, 192,65328,14336,50688, 896,64608,
72 4096,61184, 256,65264,16384,48128, 1024,64448,
73 2048,63360, 0,65520, 8192,56832, 512,64992,
74 6144,59136, 384,65136,24576,39936, 1536,63936,
75 3072,62336, 128,65392,12288,52736, 768,64736,
76 5632,59648, 352,65168,22528,41984, 1408,64064,
77 2816,62592, 96,65424,11264,53760, 704,64800,
78 7680,57600, 480,65040,30720,33792, 1920,63552,
79 3840,61568, 224,65296,15360,49664, 960,64544,
80 4608,60672, 288,65232,18432,46080, 1152,64320,
81 2304,63104, 32,65488, 9216,55808, 576,64928,
82 6656,58624, 416,65104,26624,37888, 1664,63808,
83 3328,62080, 160,65360,13312,51712, 832,64672,
84 5376,59904, 336,65184,21504,43008, 1344,64128,
85 2688,62720, 80,65440,10752,54272, 672,64832,
86 7424,57856, 464,65056,29696,34816, 1856,63616,
87 3712,61696, 208,65312,14848,50176, 928,64576,
88 4352,60928, 272,65248,17408,47104, 1088,64384,
89 2176,63232, 16,65504, 8704,56320, 544,64960,
90 6400,58880, 400,65120,25600,38912, 1600,63872,
91 3200,62208, 144,65376,12800,52224, 800,64704,
92 5888,59392, 368,65152,23552,40960, 1472,64000,
93 2944,62464, 112,65408,11776,53248, 736,64768,
94 7936,57344, 496,65024,31744,32768, 1984,63488,
95 3968,61440, 240,65280,15872,49152, 992,64512,
96 4864,60416, 304,65216,19456,45056, 1216,64256,
97 2432,62976, 48,65472, 9728,55296, 608,64896,
98 6912,58368, 432,65088,27648,36864, 1728,63744,
99 3456,61952, 176,65344,13824,51200, 864,64640
100};
101
102#define st_alaw_to_short(a) (a2s[(unsigned char)a])
103
104void ulaw_to_short(const unsigned char *ulaw,short *data,int length)
105{
106 /* Convert ulaw to shorts */
107 int i;
108
109 for (i=0; i<length; i++)
110 data[i] = st_ulaw_to_short(ulaw[i]); /* ulaw convert */
111
112}
113
114void alaw_to_short(const unsigned char *alaw,short *data,int length)
115{
116 /* Convert ulaw to shorts */
117 int i;
118
119 for (i=0; i<length; i++)
120 {
121 data[i] = st_alaw_to_short(alaw[i])-32768; /* alaw convert */
122 }
123}
124
125void shorten_to_short(unsigned char *ulaw, short *data, int length)
126{
127 /* Convert Tony Robinson's shorten format to shorts */
128 (void)ulaw;
129 (void)data;
130 (void)length;
131
132}
133
134void uchar_to_short(const unsigned char *chars,short *data,int length)
135{
136 /* Convert 8 bit data to shorts UNSIGNED CHAR */
137 int i;
138
139 for (i=0; i<length; i++)
140 {
141 data[i] = (((int)chars[i])-128)*256;
142 }
143
144}
145
146void schar_to_short(const unsigned char *chars,short *data,int length)
147{
148 /* Convert 8 bit data to shorts SIGNED CHAR */
149 int i;
150
151 for (i=0; i<length; i++)
152 data[i] = (((unsigned char)chars[i]))*256;
153
154}
155
156void short_to_uchar(const short *data,unsigned char *chars,int length)
157{
158 /* Convert shorts to 8 bit UNSIGNED CHAR */
159 int i;
160
161 for (i=0; i<length; i++)
162 chars[i] = (data[i]/256)+128;
163
164}
165
166void short_to_schar(const short *data,unsigned char *chars,int length)
167{
168 /* Convert shorts to 8 bit SIGNED CHAR */
169 int i;
170
171 for (i=0; i<length; i++)
172 chars[i] = (data[i]/256);
173
174}
175
176#if 0
177void short_to_adpcm(short *data,signed char *chars,int length)
178{
179 struct adpcm_state state;
180 state.valprev = 0;
181 state.index = 0;
182
183 adpcm_coder(data,chars,length,&state);
184
185}
186
187void adpcm_to_short(signed char *chars, short *data,int length)
188{
189 struct adpcm_state state;
190 state.valprev = 0;
191 state.index = 0;
192
193 adpcm_decoder(chars,data,length/2,&state);
194}
195#endif
196
197void short_to_ulaw(const short *data,unsigned char *ulaw,int length)
198{
199 /* Convert ulaw to shorts */
200 int i;
201
202 for (i=0; i<length; i++)
203 ulaw[i] = st_short_to_ulaw(data[i]); /* ulaw convert */
204
205}
206
207void short_to_alaw(const short *data,unsigned char *alaw,int length)
208{
209 /* Convert alaw to shorts */
210 int i;
211
212 for (i=0; i<length; i++)
213 alaw[i] = st_short_to_alaw(data[i]); /* alaw convert */
214
215}
216
217short *convert_raw_data(unsigned char *file_data,int data_length,
218 enum EST_sample_type_t sample_type,int bo)
219{
220 /* converts data in file_data to native byte order shorts */
221 /* frees file_data if its not returned */
222 short *d=NULL;
223
224 if (sample_type == st_short)
225 {
226 // d = walloc(short,data_length);
227 if (bo != EST_NATIVE_BO)
228 swap_bytes_short((short *)file_data,data_length);
229 return (short *)file_data;
230 }
231 else if (sample_type == st_mulaw)
232 {
233 d = walloc(short,data_length);
234 ulaw_to_short(file_data,d,data_length);
235 wfree(file_data);
236 return d;
237 }
238 else if (sample_type == st_alaw)
239 {
240 d = walloc(short,data_length);
241 alaw_to_short(file_data,d,data_length);
242 wfree(file_data);
243 return d;
244 }
245 else if (sample_type == st_alaw)
246 {
247 d = walloc(short,data_length);
248 alaw_to_short(file_data,d,data_length);
249 wfree(file_data);
250 return d;
251 }
252#if 0
253 else if (sample_type == st_adpcm)
254 { /* Not really checked yet */
255 d = walloc(short,data_length);
256 adpcm_to_short((signed char *)file_data,d,data_length);
257 wfree(file_data);
258 return d;
259 }
260#endif
261 else if (sample_type == st_schar)
262 {
263 d = walloc(short,data_length);
264 schar_to_short((unsigned char *)file_data,d,data_length);
265 wfree(file_data);
266 return d;
267 }
268 else if (sample_type == st_uchar)
269 {
270 d = walloc(short,data_length);
271 uchar_to_short((unsigned char *)file_data,d,data_length);
272 wfree(file_data);
273 return d;
274 }
275 else
276 EST_error("Convert raw data: unsupported sample type %s(%d)",
277 EST_sample_type_map.name(sample_type), sample_type);
278
279 /* NOTREACHED */
280 return NULL;
281}
282
283enum EST_write_status save_raw_data(FILE *fp, const short *data, int offset,
284 int num_samples, int num_channels,
285 enum EST_sample_type_t sample_type,
286 int bo)
287{
288 int i;
289 int n;
290
291 if (sample_type == st_mulaw)
292 {
293 unsigned char *ulaw = walloc(unsigned char,num_samples*num_channels);
294 short_to_ulaw(data+(offset*num_channels),
295 ulaw,num_samples*num_channels);
296 n = fwrite(ulaw,1,num_channels * num_samples,fp);
297 wfree(ulaw);
298 if (n != (num_channels * num_samples))
299 return misc_write_error;
300 }
301 else if (sample_type == st_alaw)
302 {
303 unsigned char *alaw = walloc(unsigned char,num_samples*num_channels);
304 short_to_alaw(data+(offset*num_channels),
305 alaw,num_samples*num_channels);
306 n = fwrite(alaw,1,num_channels * num_samples,fp);
307 wfree(alaw);
308 if (n != (num_channels * num_samples))
309 return misc_write_error;
310 }
311 else if (sample_type == st_ascii)
312 {
313 for (i=offset*num_channels; i < num_samples*num_channels; i++)
314 fprintf(fp,"%d\n",data[i]);
315 }
316 else if (sample_type == st_schar)
317 {
318 unsigned char *chars = walloc(unsigned char,num_samples*num_channels);
319 short_to_schar(data+(offset*num_channels),
320 chars,num_samples*num_channels);
321 n = fwrite(chars,1,num_channels * num_samples,fp);
322 wfree(chars);
323 if (n != (num_channels * num_samples))
324 return misc_write_error;
325 }
326 else if (sample_type == st_uchar)
327 {
328 unsigned char *chars = walloc(unsigned char,num_samples*num_channels);
329 short_to_uchar(data+(offset*num_channels),
330 chars,num_samples*num_channels);
331 n = fwrite(chars,1,num_channels * num_samples,fp);
332 wfree(chars);
333 if ( n != (num_channels * num_samples))
334 return misc_write_error;
335 }
336#if 0
337 else if (sample_type == st_adpcm)
338 {
339 signed char *chars = walloc(signed char,num_samples*num_channels);
340 short_to_adpcm(data+(offset*num_channels),
341 chars,num_samples*num_channels);
342 n = fwrite(chars,1,num_channels * num_samples,fp);
343 wfree(chars);
344 if ( n != (num_channels * num_samples))
345 return misc_write_error;
346 }
347#endif
348 else if (sample_type == st_short)
349 {
350 if (bo != EST_NATIVE_BO)
351 {
352 short *xdata = walloc(short,num_channels*num_samples);
353 memmove(xdata,data+(offset*num_channels),
354 num_channels*num_samples*sizeof(short));
355 swap_bytes_short(xdata,num_channels*num_samples);
356 n = fwrite(xdata, sizeof(short),num_channels * num_samples, fp);
357 wfree(xdata);
358 }
359 else
360 n = fwrite(&data[offset], sizeof(short),
361 num_channels * num_samples, fp);
362 if (n != (num_channels * num_samples))
363 return misc_write_error;
364 }
365 else
366 {
367 fprintf(stderr,"save data file: unsupported sample type\n");
368 return misc_write_error;
369 }
370 return write_ok;
371}
372
373int get_word_size(enum EST_sample_type_t sample_type)
374{
375 /* Returns word size from type */
376 int word_size;
377
378 switch (sample_type)
379 {
380 case st_unknown:
381 word_size = 2; break;
382 case st_uchar:
383 case st_schar:
384 word_size = 1; break;
385 case st_mulaw:
386 case st_alaw:
387 word_size = 1; break;
388#if 0
389 case st_adpcm: /* maybe I mean 0.5 */
390 word_size = 1; break;
391#endif
392 case st_short:
393 word_size = 2; break;
394 case st_int:
395 /* Yes I mean 4 not sizeof(int) these are machine independent defs */
396 word_size = 4; break;
397 case st_float:
398 word_size = 4; break;
399 case st_double:
400 word_size = 8; break;
401 default:
402 fprintf(stderr,"Unknown encoding format error\n");
403 word_size = 2;
404 }
405 return word_size;
406}
407
408enum EST_sample_type_t str_to_sample_type(const char *type)
409{
410 /* string to internal value */
411
412 if (streq(type,"short"))
413 return st_short;
414 if (streq(type,"shorten"))
415 return st_shorten;
416 else if ((streq(type,"ulaw")) || (streq(type,"mulaw")))
417 return st_mulaw;
418 else if ((streq(type,"char")) || (streq(type,"byte")) ||
419 (streq(type,"8bit")))
420 return st_schar;
421 else if ((streq(type,"unsignedchar")) || (streq(type,"unsignedbyte")) ||
422 (streq(type,"unsigned8bit")))
423 return st_uchar;
424 else if (streq(type,"int"))
425 return st_int;
426#if 0
427 else if (streq(type,"adpcm"))
428 return st_adpcm;
429#endif
430 else if ((streq(type,"real")) || (streq(type,"float")) ||
431 (streq(type,"real4")))
432 return st_float;
433 else if ((streq(type,"real8")) || (streq(type,"double")))
434 return st_double;
435 else if (streq(type,"alaw"))
436 return st_alaw;
437 else if (streq(type,"ascii"))
438 return st_ascii;
439 else
440 {
441 fprintf(stderr,"Unknown sample type: \"%s\"\n",type);
442 return st_unknown;
443 }
444}
445
446const char *sample_type_to_str(enum EST_sample_type_t type)
447{
448 switch (type)
449 {
450 case st_short: return "short";
451 case st_shorten: return "shorten";
452 case st_mulaw: return "ulaw";
453 case st_alaw: return "alaw";
454 case st_schar: return "char";
455 case st_uchar: return "unsignedchar";
456 case st_int: return "int";
457#if 0
458 case st_adpcm: return "adpcm";
459#endif
460 case st_float: return "float";
461 case st_double: return "double";
462 case st_ascii: return "ascii";
463 case st_unknown: return "unknown";
464 default:
465 fprintf(stderr,"Unknown sample_type %d\n",type);
466 return "very_unknown";
467 }
468}
469
470/*
471** This routine converts from linear to ulaw.
472**
473** Craig Reese: IDA/Supercomputing Research Center
474** Joe Campbell: Department of Defense
475** 29 September 1989
476**
477** References:
478** 1) CCITT Recommendation G.711 (very difficult to follow)
479** 2) "A New Digital Technique for Implementation of Any
480** Continuous PCM Companding Law," Villeret, Michel,
481** et al. 1973 IEEE Int. Conf. on Communications, Vol 1,
482** 1973, pg. 11.12-11.17
483** 3) MIL-STD-188-113,"Interoperability and Performance Standards
484** for Analog-to_Digital Conversion Techniques,"
485** 17 February 1987
486**
487** Input: Signed 16 bit linear sample
488** Output: 8 bit ulaw sample
489*/
490
491#define ZEROTRAP /* turn on the trap as per the MIL-STD */
492#define BIAS 0x84 /* define the add-in bias for 16 bit samples */
493#define CLIP 32635
494
495static unsigned char st_short_to_ulaw(short sample)
496{
497 static int exp_lut[256] = {0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
498 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
499 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
500 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
501 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
502 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
503 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
504 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
505 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
506 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
507 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
508 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
509 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
510 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
511 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
512 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7};
513 int sign, exponent, mantissa;
514 unsigned char ulawbyte;
515
516 /* Get the sample into sign-magnitude. */
517 sign = (sample >> 8) & 0x80; /* set aside the sign */
518 if ( sign != 0 ) sample = -sample; /* get magnitude */
519 if ( sample > CLIP ) sample = CLIP; /* clip the magnitude */
520
521 /* Convert from 16 bit linear to ulaw. */
522 sample = sample + BIAS;
523 exponent = exp_lut[( sample >> 7 ) & 0xFF];
524 mantissa = ( sample >> ( exponent + 3 ) ) & 0x0F;
525 ulawbyte = ~ ( sign | ( exponent << 4 ) | mantissa );
526#ifdef ZEROTRAP
527 if ( ulawbyte == 0 ) ulawbyte = 0x02; /* optional CCITT trap */
528#endif
529
530 return ulawbyte;
531}
532
533/*
534** This routine converts from ulaw to 16 bit linear.
535**
536** Craig Reese: IDA/Supercomputing Research Center
537** 29 September 1989
538**
539** References:
540** 1) CCITT Recommendation G.711 (very difficult to follow)
541** 2) MIL-STD-188-113,"Interoperability and Performance Standards
542** for Analog-to_Digital Conversion Techniques,"
543** 17 February 1987
544**
545** Input: 8 bit ulaw sample
546** Output: signed 16 bit linear sample
547*/
548
549static short st_ulaw_to_short( unsigned char ulawbyte )
550{
551 static int exp_lut[8] = { 0, 132, 396, 924, 1980, 4092, 8316, 16764 };
552 int sign, exponent, mantissa;
553 short sample;
554
555 ulawbyte = ~ ulawbyte;
556 sign = ( ulawbyte & 0x80 );
557 exponent = ( ulawbyte >> 4 ) & 0x07;
558 mantissa = ulawbyte & 0x0F;
559 sample = exp_lut[exponent] + ( mantissa << ( exponent + 3 ) );
560 if ( sign != 0 ) sample = -sample;
561
562 return sample;
563}
564
565
566/* The following is copied from sox:g711.c */
567/*
568 * This source code is a product of Sun Microsystems, Inc. and is provided
569 * for unrestricted use. Users may copy or modify this source code without
570 * charge.
571 * [ no warranties or liabilities whatsoever; see there for details. ]
572 */
573/* copy from CCITT G.711 specifications */
574static unsigned char _u2a[128] = { /* u- to A-law conversions */
575 1, 1, 2, 2, 3, 3, 4, 4,
576 5, 5, 6, 6, 7, 7, 8, 8,
577 9, 10, 11, 12, 13, 14, 15, 16,
578 17, 18, 19, 20, 21, 22, 23, 24,
579 25, 27, 29, 31, 33, 34, 35, 36,
580 37, 38, 39, 40, 41, 42, 43, 44,
581 46, 48, 49, 50, 51, 52, 53, 54,
582 55, 56, 57, 58, 59, 60, 61, 62,
583 64, 65, 66, 67, 68, 69, 70, 71,
584 72, 73, 74, 75, 76, 77, 78, 79,
585 80, 82, 83, 84, 85, 86, 87, 88,
586 89, 90, 91, 92, 93, 94, 95, 96,
587 97, 98, 99, 100, 101, 102, 103, 104,
588 105, 106, 107, 108, 109, 110, 111, 112,
589 113, 114, 115, 116, 117, 118, 119, 120,
590 121, 122, 123, 124, 125, 126, 127, 128};
591
592static unsigned char _a2u[128] = { /* A- to u-law conversions */
593 1, 3, 5, 7, 9, 11, 13, 15,
594 16, 17, 18, 19, 20, 21, 22, 23,
595 24, 25, 26, 27, 28, 29, 30, 31,
596 32, 32, 33, 33, 34, 34, 35, 35,
597 36, 37, 38, 39, 40, 41, 42, 43,
598 44, 45, 46, 47, 48, 48, 49, 49,
599 50, 51, 52, 53, 54, 55, 56, 57,
600 58, 59, 60, 61, 62, 63, 64, 64,
601 65, 66, 67, 68, 69, 70, 71, 72,
602 73, 74, 75, 76, 77, 78, 79, 80,
603 80, 81, 82, 83, 84, 85, 86, 87,
604 88, 89, 90, 91, 92, 93, 94, 95,
605 96, 97, 98, 99, 100, 101, 102, 103,
606 104, 105, 106, 107, 108, 109, 110, 111,
607 112, 113, 114, 115, 116, 117, 118, 119,
608 120, 121, 122, 123, 124, 125, 126, 127};
609
610/* A-law to u-law conversion */
611static inline unsigned char st_alaw2ulaw(
612 unsigned char aval)
613{
614 aval &= 0xff;
615 return (unsigned char) ((aval & 0x80) ? (0xFF ^ _a2u[aval ^ 0xD5]) :
616 (0x7F ^ _a2u[aval ^ 0x55]));
617}
618
619/* u-law to A-law conversion */
620static inline unsigned char st_ulaw2alaw(
621 unsigned char uval)
622{
623 uval &= 0xff;
624 return (unsigned char) ((uval & 0x80) ? (0xD5 ^ (_u2a[0xFF ^ uval] - 1)) :
625 (unsigned char) (0x55 ^ (_u2a[0x7F ^ uval] - 1)));
626}
627/* end of Sun code */
628
629/* This is somewhat simple-minded, but ... */
630static unsigned char st_short_to_alaw(short sample)
631{
632 return st_ulaw2alaw(st_short_to_ulaw(sample));
633}
634
635
636/*
637 * C O N V E R T T O I E E E E X T E N D E D
638 */
639
640/* Copyright (C) 1988-1991 Apple Computer, Inc.
641 * All rights reserved.
642 *
643 * Machine-independent I/O routines for IEEE floating-point numbers.
644 *
645 * NaN's and infinities are converted to HUGE_VAL or HUGE, which
646 * happens to be infinity on IEEE machines. Unfortunately, it is
647 * impossible to preserve NaN's in a machine-independent way.
648 * Infinities are, however, preserved on IEEE machines.
649 *
650 * These routines have been tested on the following machines:
651 * Apple Macintosh, MPW 3.1 C compiler
652 * Apple Macintosh, THINK C compiler
653 * Silicon Graphics IRIS, MIPS compiler
654 * Cray X/MP and Y/MP
655 * Digital Equipment VAX
656 *
657 *
658 * Implemented by Malcolm Slaney and Ken Turkowski.
659 *
660 * Malcolm Slaney contributions during 1988-1990 include big- and little-
661 * endian file I/O, conversion to and from Motorola's extended 80-bit
662 * floating-point format, and conversions to and from IEEE single-
663 * precision floating-point format.
664 *
665 * In 1991, Ken Turkowski implemented the conversions to and from
666 * IEEE double-precision format, added more precision to the extended
667 * conversions, and accommodated conversions involving +/- infinity,
668 * NaN's, and denormalized numbers.
669 */
670
671#ifndef HUGE_VAL
672# define HUGE_VAL HUGE
673#endif /*HUGE_VAL*/
674
675# define FloatToUnsigned(f) ((unsigned long)(((long)(f - 2147483648.0)) + 2147483647L) + 1)
676
677void ConvertToIeeeExtended(double num,unsigned char *bytes)
678{
679 int sign;
680 int expon;
681 double fMant, fsMant;
682 unsigned long hiMant, loMant;
683
684 if (num < 0) {
685 sign = 0x8000;
686 num *= -1;
687 } else {
688 sign = 0;
689 }
690
691 if (num == 0) {
692 expon = 0; hiMant = 0; loMant = 0;
693 }
694 else {
695 fMant = frexp(num, &expon);
696 if ((expon > 16384) || !(fMant < 1)) { /* Infinity or NaN */
697 expon = sign|0x7FFF; hiMant = 0; loMant = 0; /* infinity */
698 }
699 else { /* Finite */
700 expon += 16382;
701 if (expon < 0) { /* denormalized */
702 fMant = ldexp(fMant, expon);
703 expon = 0;
704 }
705 expon |= sign;
706 fMant = ldexp(fMant, 32);
707 fsMant = floor(fMant);
708 hiMant = FloatToUnsigned(fsMant);
709 fMant = ldexp(fMant - fsMant, 32);
710 fsMant = floor(fMant);
711 loMant = FloatToUnsigned(fsMant);
712 }
713 }
714
715 bytes[0] = expon >> 8;
716 bytes[1] = expon;
717 bytes[2] = hiMant >> 24;
718 bytes[3] = hiMant >> 16;
719 bytes[4] = hiMant >> 8;
720 bytes[5] = hiMant;
721 bytes[6] = loMant >> 24;
722 bytes[7] = loMant >> 16;
723 bytes[8] = loMant >> 8;
724 bytes[9] = loMant;
725}
726
727
728/*
729 * C O N V E R T F R O M I E E E E X T E N D E D
730 */
731
732/*
733 * Copyright (C) 1988-1991 Apple Computer, Inc.
734 * All rights reserved.
735 *
736 * Machine-independent I/O routines for IEEE floating-point numbers.
737 *
738 * NaN's and infinities are converted to HUGE_VAL or HUGE, which
739 * happens to be infinity on IEEE machines. Unfortunately, it is
740 * impossible to preserve NaN's in a machine-independent way.
741 * Infinities are, however, preserved on IEEE machines.
742 *
743 * These routines have been tested on the following machines:
744 * Apple Macintosh, MPW 3.1 C compiler
745 * Apple Macintosh, THINK C compiler
746 * Silicon Graphics IRIS, MIPS compiler
747 * Cray X/MP and Y/MP
748 * Digital Equipment VAX
749 *
750 *
751 * Implemented by Malcolm Slaney and Ken Turkowski.
752 *
753 * Malcolm Slaney contributions during 1988-1990 include big- and little-
754 * endian file I/O, conversion to and from Motorola's extended 80-bit
755 * floating-point format, and conversions to and from IEEE single-
756 * precision floating-point format.
757 *
758 * In 1991, Ken Turkowski implemented the conversions to and from
759 * IEEE double-precision format, added more precision to the extended
760 * conversions, and accommodated conversions involving +/- infinity,
761 * NaN's, and denormalized numbers.
762 */
763
764#ifndef HUGE_VAL
765# define HUGE_VAL HUGE
766#endif /*HUGE_VAL*/
767
768# define UnsignedToFloat(u) (((double)((long)(u - 2147483647L - 1))) + 2147483648.0)
769
770/****************************************************************
771 * Extended precision IEEE floating-point conversion routine.
772 ****************************************************************/
773
774double ConvertFromIeeeExtended(unsigned char *bytes)
775{
776 double f;
777 int expon;
778 unsigned long hiMant, loMant;
779
780 expon = ((bytes[0] & 0x7F) << 8) | (bytes[1] & 0xFF);
781 hiMant = ((unsigned long)(bytes[2] & 0xFF) << 24)
782 | ((unsigned long)(bytes[3] & 0xFF) << 16)
783 | ((unsigned long)(bytes[4] & 0xFF) << 8)
784 | ((unsigned long)(bytes[5] & 0xFF));
785 loMant = ((unsigned long)(bytes[6] & 0xFF) << 24)
786 | ((unsigned long)(bytes[7] & 0xFF) << 16)
787 | ((unsigned long)(bytes[8] & 0xFF) << 8)
788 | ((unsigned long)(bytes[9] & 0xFF));
789
790 if (expon == 0 && hiMant == 0 && loMant == 0) {
791 f = 0;
792 }
793 else {
794 if (expon == 0x7FFF) { /* Infinity or NaN */
795 f = HUGE_VAL;
796 }
797 else {
798 expon -= 16383;
799 f = ldexp((double)(hiMant), expon-=31);
800 f += ldexp((double)(loMant), expon-=32);
801 }
802 }
803
804 if (bytes[0] & 0x80)
805 return -f;
806 else
807 return f;
808}
809