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