zenilib: zenilib/jni/external/openal-soft/utils/makehrtf.c Source File

zenilib 0.5.3.0

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
makehrtf.c
Go to the documentation of this file.
1 /*
2 * HRTF utility for producing and demonstrating the process of creating an
3 * OpenAL Soft compatible HRIR data set.
4 *
5 * Copyright (C) 2011-2012 Christopher Fitzgerald
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License along
18 * with this program; if not, write to the Free Software Foundation, Inc.,
19 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Or visit: http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
22 *
23 * --------------------------------------------------------------------------
24 *
25 * A big thanks goes out to all those whose work done in the field of
26 * binaural sound synthesis using measured HRTFs makes this utility and the
27 * OpenAL Soft implementation possible.
28 *
29 * The algorithm for diffuse-field equalization was adapted from the work
30 * done by Rio Emmanuel and Larcher Veronique of IRCAM and Bill Gardner of
31 * MIT Media Laboratory. It operates as follows:
32 *
33 * 1. Take the FFT of each HRIR and only keep the magnitude responses.
34 * 2. Calculate the diffuse-field power-average of all HRIRs weighted by
35 * their contribution to the total surface area covered by their
36 * measurement.
37 * 3. Take the diffuse-field average and limit its magnitude range.
38 * 4. Equalize the responses by using the inverse of the diffuse-field
39 * average.
40 * 5. Reconstruct the minimum-phase responses.
41 * 5. Zero the DC component.
42 * 6. IFFT the result and truncate to the desired-length minimum-phase FIR.
43 *
44 * The spherical head algorithm for calculating propagation delay was adapted
45 * from the paper:
46 *
47 * Modeling Interaural Time Difference Assuming a Spherical Head
48 * Joel David Miller
49 * Music 150, Musical Acoustics, Stanford University
50 * December 2, 2001
51 *
52 * The formulae for calculating the Kaiser window metrics are from the
53 * the textbook:
54 *
55 * Discrete-Time Signal Processing
56 * Alan V. Oppenheim and Ronald W. Schafer
57 * Prentice-Hall Signal Processing Series
58 * 1999
59 */
60
61 /* Needed for 64-bit unsigned integer. */
62 #include "config.h"
63
64 #include <stdio.h>
65 #include <stdlib.h>
66 #include <stdarg.h>
67 #include <string.h>
68 #include <ctype.h>
69 #include <math.h>
70
71 // Rely (if naively) on OpenAL's header for the types used for serialization.
72 #include "AL/al.h"
73
74 #ifndef M_PI
75 #define M_PI (3.14159265358979323846)
76 #endif
77
78 #ifndef HUGE_VAL
79 #define HUGE_VAL (1.0 / 0.0)
80 #endif
81
82 // The epsilon used to maintain signal stability.
83 #define EPSILON (1e-15)
84
85 // Constants for accessing the token reader's ring buffer.
86 #define TR_RING_BITS (16)
87 #define TR_RING_SIZE (1 << TR_RING_BITS)
88 #define TR_RING_MASK (TR_RING_SIZE - 1)
89
90 // The token reader's load interval in bytes.
91 #define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
92
93 // The maximum identifier length used when processing the data set
94 // definition.
95 #define MAX_IDENT_LEN (16)
96
97 // The maximum path length used when processing filenames.
98 #define MAX_PATH_LEN (256)
99
100 // The limits for the sample 'rate' metric in the data set definition and for
101 // resampling.
102 #define MIN_RATE (32000)
103 #define MAX_RATE (96000)
104
105 // The limits for the HRIR 'points' metric in the data set definition.
106 #define MIN_POINTS (16)
107 #define MAX_POINTS (8192)
108
109 // The limits to the number of 'azimuths' listed in the data set definition.
110 #define MIN_EV_COUNT (5)
111 #define MAX_EV_COUNT (128)
112
113 // The limits for each of the 'azimuths' listed in the data set definition.
114 #define MIN_AZ_COUNT (1)
115 #define MAX_AZ_COUNT (128)
116
117 // The limits for the listener's head 'radius' in the data set definition.
118 #define MIN_RADIUS (0.05)
119 #define MAX_RADIUS (0.15)
120
121 // The limits for the 'distance' from source to listener in the definition
122 // file.
123 #define MIN_DISTANCE (0.5)
124 #define MAX_DISTANCE (2.5)
125
126 // The maximum number of channels that can be addressed for a WAVE file
127 // source listed in the data set definition.
128 #define MAX_WAVE_CHANNELS (65535)
129
130 // The limits to the byte size for a binary source listed in the definition
131 // file.
132 #define MIN_BIN_SIZE (2)
133 #define MAX_BIN_SIZE (4)
134
135 // The minimum number of significant bits for binary sources listed in the
136 // data set definition. The maximum is calculated from the byte size.
137 #define MIN_BIN_BITS (16)
138
139 // The limits to the number of significant bits for an ASCII source listed in
140 // the data set definition.
141 #define MIN_ASCII_BITS (16)
142 #define MAX_ASCII_BITS (32)
143
144 // The limits to the FFT window size override on the command line.
145 #define MIN_FFTSIZE (512)
146 #define MAX_FFTSIZE (16384)
147
148 // The limits to the equalization range limit on the command line.
149 #define MIN_LIMIT (2.0)
150 #define MAX_LIMIT (120.0)
151
152 // The limits to the truncation window size on the command line.
153 #define MIN_TRUNCSIZE (8)
154 #define MAX_TRUNCSIZE (128)
155
156 // The truncation window size must be a multiple of the below value to allow
157 // for vectorized convolution.
158 #define MOD_TRUNCSIZE (8)
159
160 // The defaults for the command line options.
161 #define DEFAULT_EQUALIZE (1)
162 #define DEFAULT_SURFACE (1)
163 #define DEFAULT_LIMIT (24.0)
164 #define DEFAULT_TRUNCSIZE (32)
165
166 // The four-character-codes for RIFF/RIFX WAVE file chunks.
167 #define FOURCC_RIFF (0x46464952) // 'RIFF'
168 #define FOURCC_RIFX (0x58464952) // 'RIFX'
169 #define FOURCC_WAVE (0x45564157) // 'WAVE'
170 #define FOURCC_FMT (0x20746D66) // 'fmt '
171 #define FOURCC_DATA (0x61746164) // 'data'
172 #define FOURCC_LIST (0x5453494C) // 'LIST'
173 #define FOURCC_WAVL (0x6C766177) // 'wavl'
174 #define FOURCC_SLNT (0x746E6C73) // 'slnt'
175
176 // The supported wave formats.
177 #define WAVE_FORMAT_PCM (0x0001)
178 #define WAVE_FORMAT_IEEE_FLOAT (0x0003)
179 #define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
180
181 // The maximum propagation delay value supported by OpenAL Soft.
182 #define MAX_HRTD (63.0)
183
184 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
185 // response protocol 01.
186 #define MHR_FORMAT ("MinPHR01")
187
188 // Byte order for the serialization routines.
189 enum ByteOrderT {
190 BO_NONE = 0,
191 BO_LITTLE ,
192 BO_BIG
193 };
194
195 // Source format for the references listed in the data set definition.
196 enum SourceFormatT {
197 SF_NONE = 0,
198 SF_WAVE , // RIFF/RIFX WAVE file.
199 SF_BIN_LE , // Little-endian binary file.
200 SF_BIN_BE , // Big-endian binary file.
201 SF_ASCII // ASCII text file.
202 };
203
204 // Element types for the references listed in the data set definition.
205 enum ElementTypeT {
206 ET_NONE = 0,
207 ET_INT , // Integer elements.
208 ET_FP // Floating-point elements.
209 };
210
211 // Desired output format from the command line.
212 enum OutputFormatT {
213 OF_NONE = 0,
214 OF_MHR , // OpenAL Soft MHR data set file.
215 OF_TABLE // OpenAL Soft built-in table file (used when compiling).
216 };
217
218 // Unsigned integer type.
219 typedef unsigned int uint;
220
221 // Serialization types. The trailing digit indicates the number of bytes.
222 typedef ALubyte uint1;
223
224 typedef ALint int4;
225 typedef ALuint uint4;
226
227 #if defined (HAVE_STDINT_H)
228 #include <stdint.h>
229
230 typedef uint64_t uint8;
231 #elif defined (HAVE___INT64)
232 typedef unsigned __int64 uint8;
233 #elif (SIZEOF_LONG == 8)
234 typedef unsigned long uint8;
235 #elif (SIZEOF_LONG_LONG == 8)
236 typedef unsigned long long uint8;
237 #endif
238
239 typedef enum ByteOrderT ByteOrderT;
240 typedef enum SourceFormatT SourceFormatT;
241 typedef enum ElementTypeT ElementTypeT;
242 typedef enum OutputFormatT OutputFormatT;
243
244 typedef struct TokenReaderT TokenReaderT;
245 typedef struct SourceRefT SourceRefT;
246 typedef struct HrirDataT HrirDataT;
247 typedef struct ResamplerT ResamplerT;
248
249 // Token reader state for parsing the data set definition.
250 struct TokenReaderT {
251 FILE * mFile;
252 const char * mName;
253 uint mLine,
254 mColumn;
255 char mRing [TR_RING_SIZE];
256 size_t mIn,
257 mOut;
258 };
259
260 // Source reference state used when loading sources.
261 struct SourceRefT {
262 SourceFormatT mFormat;
263 ElementTypeT mType;
264 uint mSize;
265 int mBits;
266 uint mChannel,
267 mSkip,
268 mOffset;
269 char mPath [MAX_PATH_LEN + 1];
270 };
271
272 // The HRIR metrics and data set used when loading, processing, and storing
273 // the resulting HRTF.
274 struct HrirDataT {
275 uint mIrRate,
276 mIrCount,
277 mIrSize,
278 mIrPoints,
279 mFftSize,
280 mEvCount,
281 mEvStart,
282 mAzCount [MAX_EV_COUNT],
283 mEvOffset [MAX_EV_COUNT];
284 double mRadius,
285 mDistance,
286 * mHrirs,
287 * mHrtds,
288 mMaxHrtd;
289 };
290
291 // The resampler metrics and FIR filter.
292 struct ResamplerT {
293 uint mP,
294 mQ,
295 mM,
296 mL;
297 double * mF;
298 };
299
300 /* Token reader routines for parsing text files. Whitespace is not
301 * significant. It can process tokens as identifiers, numbers (integer and
302 * floating-point), strings, and operators. Strings must be encapsulated by
303 * double-quotes and cannot span multiple lines.
304 */
305
306 // Setup the reader on the given file. The filename can be NULL if no error
307 // output is desired.
308 static void TrSetup (FILE * fp, const char * filename, TokenReaderT * tr) {
309 const char * name = NULL;
310 char ch;
311
312 tr -> mFile = fp;
313 name = filename;
314 // If a filename was given, store a pointer to the base name.
315 if (filename != NULL) {
316 while ((ch = (* filename)) != '\0') {
317 if ((ch == '/') || (ch == '\\'))
318 name = filename + 1;
319 filename ++;
320 }
321 }
322 tr -> mName = name;
323 tr -> mLine = 1;
324 tr -> mColumn = 1;
325 tr -> mIn = 0;
326 tr -> mOut = 0;
327 }
328
329 // Prime the reader's ring buffer, and return a result indicating that there
330 // is text to process.
331 static int TrLoad (TokenReaderT * tr) {
332 size_t toLoad, in, count;
333
334 toLoad = TR_RING_SIZE - (tr -> mIn - tr -> mOut);
335 if ((toLoad >= TR_LOAD_SIZE) && (! feof (tr -> mFile))) {
336 // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
337 toLoad = TR_LOAD_SIZE;
338 in = tr -> mIn & TR_RING_MASK;
339 count = TR_RING_SIZE - in;
340 if (count < toLoad) {
341 tr -> mIn += fread (& tr -> mRing [in], 1, count, tr -> mFile);
342 tr -> mIn += fread (& tr -> mRing [0], 1, toLoad - count, tr -> mFile);
343 } else {
344 tr -> mIn += fread (& tr -> mRing [in], 1, toLoad, tr -> mFile);
345 }
346 if (tr -> mOut >= TR_RING_SIZE) {
347 tr -> mOut -= TR_RING_SIZE;
348 tr -> mIn -= TR_RING_SIZE;
349 }
350 }
351 if (tr -> mIn > tr -> mOut)
352 return (1);
353 return (0);
354 }
355
356 // Error display routine. Only displays when the base name is not NULL.
357 static void TrErrorVA (const TokenReaderT * tr, uint line, uint column, const char * format, va_list argPtr) {
358 if (tr -> mName != NULL) {
359 fprintf (stderr, "Error (%s:%u:%u): ", tr -> mName, line, column);
360 vfprintf (stderr, format, argPtr);
361 }
362 }
363
364 // Used to display an error at a saved line/column.
365 static void TrErrorAt (const TokenReaderT * tr, uint line, uint column, const char * format, ...) {
366 va_list argPtr;
367
368 va_start (argPtr, format);
369 TrErrorVA (tr, line, column, format, argPtr);
370 va_end (argPtr);
371 }
372
373 // Used to display an error at the current line/column.
374 static void TrError (const TokenReaderT * tr, const char * format, ...) {
375 va_list argPtr;
376
377 va_start (argPtr, format);
378 TrErrorVA (tr, tr -> mLine, tr -> mColumn, format, argPtr);
379 va_end (argPtr);
380 }
381
382 // Skips to the next line.
383 static void TrSkipLine (TokenReaderT * tr) {
384 char ch;
385
386 while (TrLoad (tr)) {
387 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
388 tr -> mOut ++;
389 if (ch == '\n') {
390 tr -> mLine ++;
391 tr -> mColumn = 1;
392 break;
393 }
394 tr -> mColumn ++;
395 }
396 }
397
398 // Skips to the next token.
399 static int TrSkipWhitespace (TokenReaderT * tr) {
400 char ch;
401
402 while (TrLoad (tr)) {
403 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
404 if (isspace (ch)) {
405 tr -> mOut ++;
406 if (ch == '\n') {
407 tr -> mLine ++;
408 tr -> mColumn = 1;
409 } else {
410 tr -> mColumn ++;
411 }
412 } else if (ch == '#') {
413 TrSkipLine (tr);
414 } else {
415 return (1);
416 }
417 }
418 return (0);
419 }
420
421 // Get the line and/or column of the next token (or the end of input).
422 static void TrIndication (TokenReaderT * tr, uint * line, uint * column) {
423 TrSkipWhitespace (tr);
424 if (line != NULL)
425 (* line) = tr -> mLine;
426 if (column != NULL)
427 (* column) = tr -> mColumn;
428 }
429
430 // Checks to see if a token is the given operator. It does not display any
431 // errors and will not proceed to the next token.
432 static int TrIsOperator (TokenReaderT * tr, const char * op) {
433 size_t out, len;
434 char ch;
435
436 if (! TrSkipWhitespace (tr))
437 return (0);
438 out = tr -> mOut;
439 len = 0;
440 while ((op [len] != '\0') && (out < tr -> mIn)) {
441 ch = tr -> mRing [out & TR_RING_MASK];
442 if (ch != op [len])
443 break;
444 len ++;
445 out ++;
446 }
447 if (op [len] == '\0')
448 return (1);
449 return (0);
450 }
451
452 /* The TrRead*() routines obtain the value of a matching token type. They
453 * display type, form, and boundary errors and will proceed to the next
454 * token.
455 */
456
457 // Reads and validates an identifier token.
458 static int TrReadIdent (TokenReaderT * tr, const uint maxLen, char * ident) {
459 uint col, len;
460 char ch;
461
462 col = tr -> mColumn;
463 if (TrSkipWhitespace (tr)) {
464 col = tr -> mColumn;
465 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
466 if ((ch == '_') || isalpha (ch)) {
467 len = 0;
468 do {
469 if (len < maxLen)
470 ident [len] = ch;
471 len ++;
472 tr -> mOut ++;
473 if (! TrLoad (tr))
474 break;
475 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
476 } while ((ch == '_') || isdigit (ch) || isalpha (ch));
477 tr -> mColumn += len;
478 if (len > maxLen) {
479 TrErrorAt (tr, tr -> mLine, col, "Identifier is too long.\n");
480 return (0);
481 }
482 ident [len] = '\0';
483 return (1);
484 }
485 }
486 TrErrorAt (tr, tr -> mLine, col, "Expected an identifier.\n");
487 return (0);
488 }
489
490 // Reads and validates (including bounds) an integer token.
491 static int TrReadInt (TokenReaderT * tr, const int loBound, const int hiBound, int * value) {
492 uint col, digis, len;
493 char ch, temp [64 + 1];
494
495 col = tr -> mColumn;
496 if (TrSkipWhitespace (tr)) {
497 col = tr -> mColumn;
498 len = 0;
499 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
500 if ((ch == '+') || (ch == '-')) {
501 temp [len] = ch;
502 len ++;
503 tr -> mOut ++;
504 }
505 digis = 0;
506 while (TrLoad (tr)) {
507 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
508 if (! isdigit (ch))
509 break;
510 if (len < 64)
511 temp [len] = ch;
512 len ++;
513 digis ++;
514 tr -> mOut ++;
515 }
516 tr -> mColumn += len;
517 if ((digis > 0) && (ch != '.') && (! isalpha (ch))) {
518 if (len > 64) {
519 TrErrorAt (tr, tr -> mLine, col, "Integer is too long.");
520 return (0);
521 }
522 temp [len] = '\0';
523 (* value) = strtol (temp, NULL, 10);
524 if (((* value) < loBound) || ((* value) > hiBound)) {
525 TrErrorAt (tr, tr -> mLine, col, "Expected a value from %d to %d.\n", loBound, hiBound);
526 return (0);
527 }
528 return (1);
529 }
530 }
531 TrErrorAt (tr, tr -> mLine, col, "Expected an integer.\n");
532 return (0);
533 }
534
535 // Reads and validates (including bounds) a float token.
536 static int TrReadFloat (TokenReaderT * tr, const double loBound, const double hiBound, double * value) {
537 uint col, digis, len;
538 char ch, temp [64 + 1];
539
540 col = tr -> mColumn;
541 if (TrSkipWhitespace (tr)) {
542 col = tr -> mColumn;
543 len = 0;
544 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
545 if ((ch == '+') || (ch == '-')) {
546 temp [len] = ch;
547 len ++;
548 tr -> mOut ++;
549 }
550 digis = 0;
551 while (TrLoad (tr)) {
552 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
553 if (! isdigit (ch))
554 break;
555 if (len < 64)
556 temp [len] = ch;
557 len ++;
558 digis ++;
559 tr -> mOut ++;
560 }
561 if (ch == '.') {
562 if (len < 64)
563 temp [len] = ch;
564 len ++;
565 tr -> mOut ++;
566 }
567 while (TrLoad (tr)) {
568 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
569 if (! isdigit (ch))
570 break;
571 if (len < 64)
572 temp [len] = ch;
573 len ++;
574 digis ++;
575 tr -> mOut ++;
576 }
577 if (digis > 0) {
578 if ((ch == 'E') || (ch == 'e')) {
579 if (len < 64)
580 temp [len] = ch;
581 len ++;
582 digis = 0;
583 tr -> mOut ++;
584 if ((ch == '+') || (ch == '-')) {
585 if (len < 64)
586 temp [len] = ch;
587 len ++;
588 tr -> mOut ++;
589 }
590 while (TrLoad (tr)) {
591 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
592 if (! isdigit (ch))
593 break;
594 if (len < 64)
595 temp [len] = ch;
596 len ++;
597 digis ++;
598 tr -> mOut ++;
599 }
600 }
601 tr -> mColumn += len;
602 if ((digis > 0) && (ch != '.') && (! isalpha (ch))) {
603 if (len > 64) {
604 TrErrorAt (tr, tr -> mLine, col, "Float is too long.");
605 return (0);
606 }
607 temp [len] = '\0';
608 (* value) = strtod (temp, NULL);
609 if (((* value) < loBound) || ((* value) > hiBound)) {
610 TrErrorAt (tr, tr -> mLine, col, "Expected a value from %f to %f.\n", loBound, hiBound);
611 return (0);
612 }
613 return (1);
614 }
615 } else {
616 tr -> mColumn += len;
617 }
618 }
619 TrErrorAt (tr, tr -> mLine, col, "Expected a float.\n");
620 return (0);
621 }
622
623 // Reads and validates a string token.
624 static int TrReadString (TokenReaderT * tr, const uint maxLen, char * text) {
625 uint col, len;
626 char ch;
627
628 col = tr -> mColumn;
629 if (TrSkipWhitespace (tr)) {
630 col = tr -> mColumn;
631 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
632 if (ch == '\"') {
633 tr -> mOut ++;
634 len = 0;
635 while (TrLoad (tr)) {
636 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
637 tr -> mOut ++;
638 if (ch == '\"')
639 break;
640 if (ch == '\n') {
641 TrErrorAt (tr, tr -> mLine, col, "Unterminated string at end of line.\n");
642 return (0);
643 }
644 if (len < maxLen)
645 text [len] = ch;
646 len ++;
647 }
648 if (ch != '\"') {
649 tr -> mColumn += 1 + len;
650 TrErrorAt (tr, tr -> mLine, col, "Unterminated string at end of input.\n");
651 return (0);
652 }
653 tr -> mColumn += 2 + len;
654 if (len > maxLen) {
655 TrErrorAt (tr, tr -> mLine, col, "String is too long.\n");
656 return (0);
657 }
658 text [len] = '\0';
659 return (1);
660 }
661 }
662 TrErrorAt (tr, tr -> mLine, col, "Expected a string.\n");
663 return (0);
664 }
665
666 // Reads and validates the given operator.
667 static int TrReadOperator (TokenReaderT * tr, const char * op) {
668 uint col, len;
669 char ch;
670
671 col = tr -> mColumn;
672 if (TrSkipWhitespace (tr)) {
673 col = tr -> mColumn;
674 len = 0;
675 while ((op [len] != '\0') && TrLoad (tr)) {
676 ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
677 if (ch != op [len])
678 break;
679 len ++;
680 tr -> mOut ++;
681 }
682 tr -> mColumn += len;
683 if (op [len] == '\0')
684 return (1);
685 }
686 TrErrorAt (tr, tr -> mLine, col, "Expected '%s' operator.\n", op);
687 return (0);
688 }
689
690 /* Performs a string substitution. Any case-insensitive occurrences of the
691 * pattern string are replaced with the replacement string. The result is
692 * truncated if necessary.
693 */
694 static int StrSubst (const char * in, const char * pat, const char * rep, const size_t maxLen, char * out) {
695 size_t inLen, patLen, repLen;
696 size_t si, di;
697 int truncated;
698
699 inLen = strlen (in);
700 patLen = strlen (pat);
701 repLen = strlen (rep);
702 si = 0;
703 di = 0;
704 truncated = 0;
705 while ((si < inLen) && (di < maxLen)) {
706 if (patLen <= (inLen - si)) {
707 if (strncasecmp (& in [si], pat, patLen) == 0) {
708 if (repLen > (maxLen - di)) {
709 repLen = maxLen - di;
710 truncated = 1;
711 }
712 strncpy (& out [di], rep, repLen);
713 si += patLen;
714 di += repLen;
715 }
716 }
717 out [di] = in [si];
718 si ++;
719 di ++;
720 }
721 if (si < inLen)
722 truncated = 1;
723 out [di] = '\0';
724 return (! truncated);
725 }
726
727 // Provide missing math routines for MSVC.
728 #ifdef _MSC_VER
729 static double round (double val) {
730 if (val < 0.0)
731 return (ceil (val - 0.5));
732 return (floor (val + 0.5));
733 }
734
735 static double fmin (double a, double b) {
736 return ((a < b) ? a : b);
737 }
738
739 static double fmax (double a, double b) {
740 return ((a > b) ? a : b);
741 }
742 #endif
743
744 // Simple clamp routine.
745 static double Clamp (const double val, const double lower, const double upper) {
746 return (fmin (fmax (val, lower), upper));
747 }
748
749 // Performs linear interpolation.
750 static double Lerp (const double a, const double b, const double f) {
751 return (a + (f * (b - a)));
752 }
753
754 // Performs a high-passed triangular probability density function dither from
755 // a double to an integer. It assumes the input sample is already scaled.
756 static int HpTpdfDither (const double in, int * hpHist) {
757 const double PRNG_SCALE = 1.0 / (RAND_MAX + 1.0);
758 int prn;
759 double out;
760
761 prn = rand ();
762 out = round (in + (PRNG_SCALE * (prn - (* hpHist))));
763 (* hpHist) = prn;
764 return ((int) out);
765 }
766
767 // Allocates an array of doubles.
768 static double * CreateArray (const size_t n) {
769 double * a = NULL;
770
771 a = (double *) calloc (n, sizeof (double));
772 if (a == NULL) {
773 fprintf (stderr, "Error: Out of memory.\n");
774 exit (-1);
775 }
776 return (a);
777 }
778
779 // Frees an array of doubles.
780 static void DestroyArray (const double * a) {
781 free ((void *) a);
782 }
783
784 // Complex number routines. All outputs must be non-NULL.
785
786 // Magnitude/absolute value.
787 static double ComplexAbs (const double r, const double i) {
788 return (sqrt ((r * r) + (i * i)));
789 }
790
791 // Multiply.
792 static void ComplexMul (const double aR, const double aI, const double bR, const double bI, double * outR, double * outI) {
793 (* outR) = (aR * bR) - (aI * bI);
794 (* outI) = (aI * bR) + (aR * bI);
795 }
796
797 // Base-e exponent.
798 static void ComplexExp (const double inR, const double inI, double * outR, double * outI) {
799 double e;
800
801 e = exp (inR);
802 (* outR) = e * cos (inI);
803 (* outI) = e * sin (inI);
804 }
805
806 /* Fast Fourier transform routines. The number of points must be a power of
807 * two. In-place operation is possible only if both the real and imaginary
808 * parts are in-place together.
809 */
810
811 // Performs bit-reversal ordering.
812 static void FftArrange (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
813 uint rk, k, m;
814 double tempR, tempI;
815
816 if ((inR == outR) && (inI == outI)) {
817 // Handle in-place arrangement.
818 rk = 0;
819 for (k = 0; k < n; k ++) {
820 if (rk > k) {
821 tempR = inR [rk];
822 tempI = inI [rk];
823 outR [rk] = inR [k];
824 outI [rk] = inI [k];
825 outR [k] = tempR;
826 outI [k] = tempI;
827 }
828 m = n;
829 while (rk & (m >>= 1))
830 rk &= ~m;
831 rk |= m;
832 }
833 } else {
834 // Handle copy arrangement.
835 rk = 0;
836 for (k = 0; k < n; k ++) {
837 outR [rk] = inR [k];
838 outI [rk] = inI [k];
839 m = n;
840 while (rk & (m >>= 1))
841 rk &= ~m;
842 rk |= m;
843 }
844 }
845 }
846
847 // Performs the summation.
848 static void FftSummation (const uint n, const double s, double * re, double * im) {
849 double pi;
850 uint m, m2;
851 double vR, vI, wR, wI;
852 uint i, k, mk;
853 double tR, tI;
854
855 pi = s * M_PI;
856 for (m = 1, m2 = 2; m < n; m <<= 1, m2 <<= 1) {
857 // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
858 vR = sin (0.5 * pi / m);
859 vR = -2.0 * vR * vR;
860 vI = -sin (pi / m);
861 // w = Complex (1.0, 0.0)
862 wR = 1.0;
863 wI = 0.0;
864 for (i = 0; i < m; i ++) {
865 for (k = i; k < n; k += m2) {
866 mk = k + m;
867 // t = ComplexMul (w, out [km2])
868 tR = (wR * re [mk]) - (wI * im [mk]);
869 tI = (wR * im [mk]) + (wI * re [mk]);
870 // out [mk] = ComplexSub (out [k], t)
871 re [mk] = re [k] - tR;
872 im [mk] = im [k] - tI;
873 // out [k] = ComplexAdd (out [k], t)
874 re [k] += tR;
875 im [k] += tI;
876 }
877 // t = ComplexMul (v, w)
878 tR = (vR * wR) - (vI * wI);
879 tI = (vR * wI) + (vI * wR);
880 // w = ComplexAdd (w, t)
881 wR += tR;
882 wI += tI;
883 }
884 }
885 }
886
887 // Performs a forward FFT.
888 static void FftForward (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
889 FftArrange (n, inR, inI, outR, outI);
890 FftSummation (n, 1.0, outR, outI);
891 }
892
893 // Performs an inverse FFT.
894 static void FftInverse (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
895 double f;
896 uint i;
897
898 FftArrange (n, inR, inI, outR, outI);
899 FftSummation (n, -1.0, outR, outI);
900 f = 1.0 / n;
901 for (i = 0; i < n; i ++) {
902 outR [i] *= f;
903 outI [i] *= f;
904 }
905 }
906
907 /* Calculate the complex helical sequence (or discrete-time analytical
908 * signal) of the given input using the Hilbert transform. Given the
909 * negative natural logarithm of a signal's magnitude response, the imaginary
910 * components can be used as the angles for minimum-phase reconstruction.
911 */
912 static void Hilbert (const uint n, const double * in, double * outR, double * outI) {
913 uint i;
914
915 if (in == outR) {
916 // Handle in-place operation.
917 for (i = 0; i < n; i ++)
918 outI [i] = 0.0;
919 } else {
920 // Handle copy operation.
921 for (i = 0; i < n; i ++) {
922 outR [i] = in [i];
923 outI [i] = 0.0;
924 }
925 }
926 FftForward (n, outR, outI, outR, outI);
927 /* Currently the Fourier routines operate only on point counts that are
928 * powers of two. If that changes and n is odd, the following conditional
929 * should be: i < (n + 1) / 2.
930 */
931 for (i = 1; i < (n / 2); i ++) {
932 outR [i] *= 2.0;
933 outI [i] *= 2.0;
934 }
935 // If n is odd, the following increment should be skipped.
936 i ++;
937 for (; i < n; i ++) {
938 outR [i] = 0.0;
939 outI [i] = 0.0;
940 }
941 FftInverse (n, outR, outI, outR, outI);
942 }
943
944 /* Calculate the magnitude response of the given input. This is used in
945 * place of phase decomposition, since the phase residuals are discarded for
946 * minimum phase reconstruction. The mirrored half of the response is also
947 * discarded.
948 */
949 static void MagnitudeResponse (const uint n, const double * inR, const double * inI, double * out) {
950 const uint m = 1 + (n / 2);
951 uint i;
952
953 for (i = 0; i < m; i ++)
954 out [i] = fmax (ComplexAbs (inR [i], inI [i]), EPSILON);
955 }
956
957 /* Apply a range limit (in dB) to the given magnitude response. This is used
958 * to adjust the effects of the diffuse-field average on the equalization
959 * process.
960 */
961 static void LimitMagnitudeResponse (const uint n, const double limit, const double * in, double * out) {
962 const uint m = 1 + (n / 2);
963 double halfLim;
964 uint i, lower, upper;
965 double ave;
966
967 halfLim = limit / 2.0;
968 // Convert the response to dB.
969 for (i = 0; i < m; i ++)
970 out [i] = 20.0 * log10 (in [i]);
971 // Use six octaves to calculate the average magnitude of the signal.
972 lower = ((uint) ceil (n / pow (2.0, 8.0))) - 1;
973 upper = ((uint) floor (n / pow (2.0, 2.0))) - 1;
974 ave = 0.0;
975 for (i = lower; i <= upper; i ++)
976 ave += out [i];
977 ave /= upper - lower + 1;
978 // Keep the response within range of the average magnitude.
979 for (i = 0; i < m; i ++)
980 out [i] = Clamp (out [i], ave - halfLim, ave + halfLim);
981 // Convert the response back to linear magnitude.
982 for (i = 0; i < m; i ++)
983 out [i] = pow (10.0, out [i] / 20.0);
984 }
985
986 /* Reconstructs the minimum-phase component for the given magnitude response
987 * of a signal. This is equivalent to phase recomposition, sans the missing
988 * residuals (which were discarded). The mirrored half of the response is
989 * reconstructed.
990 */
991 static void MinimumPhase (const uint n, const double * in, double * outR, double * outI) {
992 const uint m = 1 + (n / 2);
993 double * mags = NULL;
994 uint i;
995 double aR, aI;
996
997 mags = CreateArray (n);
998 for (i = 0; i < m; i ++) {
999 mags [i] = fmax (in [i], EPSILON);
1000 outR [i] = -log (mags [i]);
1001 }
1002 for (; i < n; i ++) {
1003 mags [i] = mags [n - i];
1004 outR [i] = outR [n - i];
1005 }
1006 Hilbert (n, outR, outR, outI);
1007 // Remove any DC offset the filter has.
1008 outR [0] = 0.0;
1009 outI [0] = 0.0;
1010 for (i = 1; i < n; i ++) {
1011 ComplexExp (0.0, outI [i], & aR, & aI);
1012 ComplexMul (mags [i], 0.0, aR, aI, & outR [i], & outI [i]);
1013 }
1014 DestroyArray (mags);
1015 }
1016
1017 /* This is the normalized cardinal sine (sinc) function.
1018 *
1019 * sinc(x) = { 1, x = 0
1020 * { sin(pi x) / (pi x), otherwise.
1021 */
1022 static double Sinc (const double x) {
1023 if (fabs (x) < EPSILON)
1024 return (1.0);
1025 return (sin (M_PI * x) / (M_PI * x));
1026 }
1027
1028 /* The zero-order modified Bessel function of the first kind, used for the
1029 * Kaiser window.
1030 *
1031 * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
1032 * = sum_{k=0}^inf ((x / 2)^k / k!)^2
1033 */
1034 static double BesselI_0 (const double x) {
1035 double term, sum, x2, y, last_sum;
1036 int k;
1037
1038 // Start at k=1 since k=0 is trivial.
1039 term = 1.0;
1040 sum = 1.0;
1041 x2 = x / 2.0;
1042 k = 1;
1043 // Let the integration converge until the term of the sum is no longer
1044 // significant.
1045 do {
1046 y = x2 / k;
1047 k ++;
1048 last_sum = sum;
1049 term *= y * y;
1050 sum += term;
1051 } while (sum != last_sum);
1052 return (sum);
1053 }
1054
1055 /* Calculate a Kaiser window from the given beta value and a normalized k
1056 * [-1, 1].
1057 *
1058 * w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
1059 * { 0, elsewhere.
1060 *
1061 * Where k can be calculated as:
1062 *
1063 * k = i / l, where -l <= i <= l.
1064 *
1065 * or:
1066 *
1067 * k = 2 i / M - 1, where 0 <= i <= M.
1068 */
1069 static double Kaiser (const double b, const double k) {
1070 double k2;
1071
1072 k2 = Clamp (k, -1.0, 1.0);
1073 if ((k < -1.0) || (k > 1.0))
1074 return (0.0);
1075 k2 *= k2;
1076 return (BesselI_0 (b * sqrt (1.0 - k2)) / BesselI_0 (b));
1077 }
1078
1079 // Calculates the greatest common divisor of a and b.
1080 static uint Gcd (const uint a, const uint b) {
1081 uint x, y, z;
1082
1083 x = a;
1084 y = b;
1085 while (y > 0) {
1086 z = y;
1087 y = x % y;
1088 x = z;
1089 }
1090 return (x);
1091 }
1092
1093 /* Calculates the size (order) of the Kaiser window. Rejection is in dB and
1094 * the transition width is normalized frequency (0.5 is nyquist).
1095 *
1096 * M = { ceil((r - 7.95) / (2.285 2 pi f_t)), r > 21
1097 * { ceil(5.79 / 2 pi f_t), r <= 21.
1098 *
1099 */
1100 static uint CalcKaiserOrder (const double rejection, const double transition) {
1101 double w_t;
1102
1103 w_t = 2.0 * M_PI * transition;
1104 if (rejection > 21.0)
1105 return ((uint) ceil ((rejection - 7.95) / (2.285 * w_t)));
1106 return ((uint) ceil (5.79 / w_t));
1107 }
1108
1109 // Calculates the beta value of the Kaiser window. Rejection is in dB.
1110 static double CalcKaiserBeta (const double rejection) {
1111 if (rejection > 50.0)
1112 return (0.1102 * (rejection - 8.7));
1113 else if (rejection >= 21.0)
1114 return ((0.5842 * pow (rejection - 21.0, 0.4)) +
1115 (0.07886 * (rejection - 21.0)));
1116 else
1117 return (0.0);
1118 }
1119
1120 /* Calculates a point on the Kaiser-windowed sinc filter for the given half-
1121 * width, beta, gain, and cutoff. The point is specified in non-normalized
1122 * samples, from 0 to M, where M = (2 l + 1).
1123 *
1124 * w(k) 2 p f_t sinc(2 f_t x)
1125 *
1126 * x -- centered sample index (i - l)
1127 * k -- normalized and centered window index (x / l)
1128 * w(k) -- window function (Kaiser)
1129 * p -- gain compensation factor when sampling
1130 * f_t -- normalized center frequency (or cutoff; 0.5 is nyquist)
1131 */
1132 static double SincFilter (const int l, const double b, const double gain, const double cutoff, const int i) {
1133 return (Kaiser (b, ((double) (i - l)) / l) * 2.0 * gain * cutoff * Sinc (2.0 * cutoff * (i - l)));
1134 }
1135
1136 /* This is a polyphase sinc-filtered resampler.
1137 *
1138 * Upsample Downsample
1139 *
1140 * p/q = 3/2 p/q = 3/5
1141 *
1142 * M-+-+-+-> M-+-+-+->
1143 * -------------------+ ---------------------+
1144 * p s * f f f f|f| | p s * f f f f f |
1145 * | 0 * 0 0 0|0|0 | | 0 * 0 0 0 0|0| |
1146 * v 0 * 0 0|0|0 0 | v 0 * 0 0 0|0|0 |
1147 * s * f|f|f f f | s * f f|f|f f |
1148 * 0 * |0|0 0 0 0 | 0 * 0|0|0 0 0 |
1149 * --------+=+--------+ 0 * |0|0 0 0 0 |
1150 * d . d .|d|. d . d ----------+=+--------+
1151 * d . . . .|d|. . . .
1152 * q->
1153 * q-+-+-+->
1154 *
1155 * P_f(i,j) = q i mod p + pj
1156 * P_s(i,j) = floor(q i / p) - j
1157 * d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
1158 * { f[P_f(i,j)] s[P_s(i,j)], P_f(i,j) < M
1159 * { 0, P_f(i,j) >= M. }
1160 */
1161
1162 // Calculate the resampling metrics and build the Kaiser-windowed sinc filter
1163 // that's used to cut frequencies above the destination nyquist.
1164 static void ResamplerSetup (ResamplerT * rs, const uint srcRate, const uint dstRate) {
1165 uint gcd, l;
1166 double cutoff, width, beta;
1167 int i;
1168
1169 gcd = Gcd (srcRate, dstRate);
1170 rs -> mP = dstRate / gcd;
1171 rs -> mQ = srcRate / gcd;
1172 /* The cutoff is adjusted by half the transition width, so the transition
1173 * ends before the nyquist (0.5). Both are scaled by the downsampling
1174 * factor.
1175 */
1176 if (rs -> mP > rs -> mQ) {
1177 cutoff = 0.45 / rs -> mP;
1178 width = 0.1 / rs -> mP;
1179 } else {
1180 cutoff = 0.45 / rs -> mQ;
1181 width = 0.1 / rs -> mQ;
1182 }
1183 // A rejection of -180 dB is used for the stop band.
1184 l = CalcKaiserOrder (180.0, width) / 2;
1185 beta = CalcKaiserBeta (180.0);
1186 rs -> mM = (2 * l) + 1;
1187 rs -> mL = l;
1188 rs -> mF = CreateArray (rs -> mM);
1189 for (i = 0; i < ((int) rs -> mM); i ++)
1190 rs -> mF [i] = SincFilter ((int) l, beta, rs -> mP, cutoff, i);
1191 }
1192
1193 // Clean up after the resampler.
1194 static void ResamplerClear (ResamplerT * rs) {
1195 DestroyArray (rs -> mF);
1196 rs -> mF = NULL;
1197 }
1198
1199 // Perform the upsample-filter-downsample resampling operation using a
1200 // polyphase filter implementation.
1201 static void ResamplerRun (ResamplerT * rs, const uint inN, const double * in, const uint outN, double * out) {
1202 const uint p = rs -> mP, q = rs -> mQ, m = rs -> mM, l = rs -> mL;
1203 const double * f = rs -> mF;
1204 double * work = NULL;
1205 uint i;
1206 double r;
1207 uint j_f, j_s;
1208
1209 // Handle in-place operation.
1210 if (in == out)
1211 work = CreateArray (outN);
1212 else
1213 work = out;
1214 // Resample the input.
1215 for (i = 0; i < outN; i ++) {
1216 r = 0.0;
1217 // Input starts at l to compensate for the filter delay. This will
1218 // drop any build-up from the first half of the filter.
1219 j_f = (l + (q * i)) % p;
1220 j_s = (l + (q * i)) / p;
1221 while (j_f < m) {
1222 // Only take input when 0 <= j_s < inN. This single unsigned
1223 // comparison catches both cases.
1224 if (j_s < inN)
1225 r += f [j_f] * in [j_s];
1226 j_f += p;
1227 j_s --;
1228 }
1229 work [i] = r;
1230 }
1231 // Clean up after in-place operation.
1232 if (in == out) {
1233 for (i = 0; i < outN; i ++)
1234 out [i] = work [i];
1235 DestroyArray (work);
1236 }
1237 }
1238
1239 // Read a binary value of the specified byte order and byte size from a file,
1240 // storing it as a 32-bit unsigned integer.
1241 static int ReadBin4 (FILE * fp, const char * filename, const ByteOrderT order, const uint bytes, uint4 * out) {
1242 uint1 in [4];
1243 uint4 accum;
1244 uint i;
1245
1246 if (fread (in, 1, bytes, fp) != bytes) {
1247 fprintf (stderr, "Error: Bad read from file '%s'.\n", filename);
1248 return (0);
1249 }
1250 accum = 0;
1251 switch (order) {
1252 case BO_LITTLE :
1253 for (i = 0; i < bytes; i ++)
1254 accum = (accum << 8) | in [bytes - i - 1];
1255 break;
1256 case BO_BIG :
1257 for (i = 0; i < bytes; i ++)
1258 accum = (accum << 8) | in [i];
1259 break;
1260 default :
1261 break;
1262 }
1263 (* out) = accum;
1264 return (1);
1265 }
1266
1267 // Read a binary value of the specified byte order from a file, storing it as
1268 // a 64-bit unsigned integer.
1269 static int ReadBin8 (FILE * fp, const char * filename, const ByteOrderT order, uint8 * out) {
1270 uint1 in [8];
1271 uint8 accum;
1272 uint i;
1273
1274 if (fread (in, 1, 8, fp) != 8) {
1275 fprintf (stderr, "Error: Bad read from file '%s'.\n", filename);
1276 return (0);
1277 }
1278 accum = 0ULL;
1279 switch (order) {
1280 case BO_LITTLE :
1281 for (i = 0; i < 8; i ++)
1282 accum = (accum << 8) | in [8 - i - 1];
1283 break;
1284 case BO_BIG :
1285 for (i = 0; i < 8; i ++)
1286 accum = (accum << 8) | in [i];
1287 break;
1288 default :
1289 break;
1290 }
1291 (* out) = accum;
1292 return (1);
1293 }
1294
1295 // Write an ASCII string to a file.
1296 static int WriteAscii (const char * out, FILE * fp, const char * filename) {
1297 size_t len;
1298
1299 len = strlen (out);
1300 if (fwrite (out, 1, len, fp) != len) {
1301 fclose (fp);
1302 fprintf (stderr, "Error: Bad write to file '%s'.\n", filename);
1303 return (0);
1304 }
1305 return (1);
1306 }
1307
1308 // Write a binary value of the given byte order and byte size to a file,
1309 // loading it from a 32-bit unsigned integer.
1310 static int WriteBin4 (const ByteOrderT order, const uint bytes, const uint4 in, FILE * fp, const char * filename) {
1311 uint1 out [4];
1312 uint i;
1313
1314 switch (order) {
1315 case BO_LITTLE :
1316 for (i = 0; i < bytes; i ++)
1317 out [i] = (in >> (i * 8)) & 0x000000FF;
1318 break;
1319 case BO_BIG :
1320 for (i = 0; i < bytes; i ++)
1321 out [bytes - i - 1] = (in >> (i * 8)) & 0x000000FF;
1322 break;
1323 default :
1324 break;
1325 }
1326 if (fwrite (out, 1, bytes, fp) != bytes) {
1327 fprintf (stderr, "Error: Bad write to file '%s'.\n", filename);
1328 return (0);
1329 }
1330 return (1);
1331 }
1332
1333 /* Read a binary value of the specified type, byte order, and byte size from
1334 * a file, converting it to a double. For integer types, the significant
1335 * bits are used to normalize the result. The sign of bits determines
1336 * whether they are padded toward the MSB (negative) or LSB (positive).
1337 * Floating-point types are not normalized.
1338 */
1339 static int ReadBinAsDouble (FILE * fp, const char * filename, const ByteOrderT order, const ElementTypeT type, const uint bytes, const int bits, double * out) {
1340 union {
1341 uint4 ui;
1342 int4 i;
1343 float f;
1344 } v4;
1345 union {
1346 uint8 ui;
1347 double f;
1348 } v8;
1349
1350 (* out) = 0.0;
1351 if (bytes > 4) {
1352 if (! ReadBin8 (fp, filename, order, & v8 . ui))
1353 return (0);
1354 if (type == ET_FP)
1355 (* out) = v8 . f;
1356 } else {
1357 if (! ReadBin4 (fp, filename, order, bytes, & v4 . ui))
1358 return (0);
1359 if (type == ET_FP) {
1360 (* out) = (double) v4 . f;
1361 } else {
1362 if (bits > 0)
1363 v4 . ui >>= (8 * bytes) - ((uint) bits);
1364 else
1365 v4 . ui &= (0xFFFFFFFF >> (32 + bits));
1366 if (v4 . ui & ((uint) (1 << (abs (bits) - 1))))
1367 v4 . ui |= (0xFFFFFFFF << abs (bits));
1368 (* out) = v4 . i / ((double) (1 << (abs (bits) - 1)));
1369 }
1370 }
1371 return (1);
1372 }
1373
1374 /* Read an ascii value of the specified type from a file, converting it to a
1375 * double. For integer types, the significant bits are used to normalize the
1376 * result. The sign of the bits should always be positive. This also skips
1377 * up to one separator character before the element itself.
1378 */
1379 static int ReadAsciiAsDouble (TokenReaderT * tr, const char * filename, const ElementTypeT type, const uint bits, double * out) {
1380 int v;
1381
1382 if (TrIsOperator (tr, ","))
1383 TrReadOperator (tr, ",");
1384 else if (TrIsOperator (tr, ":"))
1385 TrReadOperator (tr, ":");
1386 else if (TrIsOperator (tr, ";"))
1387 TrReadOperator (tr, ";");
1388 else if (TrIsOperator (tr, "|"))
1389 TrReadOperator (tr, "|");
1390 if (type == ET_FP) {
1391 if (! TrReadFloat (tr, -HUGE_VAL, HUGE_VAL, out)) {
1392 fprintf (stderr, "Error: Bad read from file '%s'.\n", filename);
1393 return (0);
1394 }
1395 } else {
1396 if (! TrReadInt (tr, -(1 << (bits - 1)), (1 << (bits - 1)) - 1, & v)) {
1397 fprintf (stderr, "Error: Bad read from file '%s'.\n", filename);
1398 return (0);
1399 }
1400 (* out) = v / ((double) ((1 << (bits - 1)) - 1));
1401 }
1402 return (1);
1403 }
1404
1405 // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
1406 // the source parameters and data set metrics.
1407 static int ReadWaveFormat (FILE * fp, const ByteOrderT order, const uint hrirRate, SourceRefT * src) {
1408 uint4 fourCC, chunkSize;
1409 uint4 format, channels, rate, dummy, block, size, bits;
1410
1411 chunkSize = 0;
1412 do {
1413 if (chunkSize > 0)
1414 fseek (fp, (long) chunkSize, SEEK_CUR);
1415 if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
1416 (! ReadBin4 (fp, src -> mPath, order, 4, & chunkSize)))
1417 return (0);
1418 } while (fourCC != FOURCC_FMT);
1419 if ((! ReadBin4 (fp, src -> mPath, order, 2, & format)) ||
1420 (! ReadBin4 (fp, src -> mPath, order, 2, & channels)) ||
1421 (! ReadBin4 (fp, src -> mPath, order, 4, & rate)) ||
1422 (! ReadBin4 (fp, src -> mPath, order, 4, & dummy)) ||
1423 (! ReadBin4 (fp, src -> mPath, order, 2, & block)))
1424 return (0);
1425 block /= channels;
1426 if (chunkSize > 14) {
1427 if (! ReadBin4 (fp, src -> mPath, order, 2, & size))
1428 return (0);
1429 size /= 8;
1430 if (block > size)
1431 size = block;
1432 } else {
1433 size = block;
1434 }
1435 if (format == WAVE_FORMAT_EXTENSIBLE) {
1436 fseek (fp, 2, SEEK_CUR);
1437 if (! ReadBin4 (fp, src -> mPath, order, 2, & bits))
1438 return (0);
1439 if (bits == 0)
1440 bits = 8 * size;
1441 fseek (fp, 4, SEEK_CUR);
1442 if (! ReadBin4 (fp, src -> mPath, order, 2, & format))
1443 return (0);
1444 fseek (fp, (long) (chunkSize - 26), SEEK_CUR);
1445 } else {
1446 bits = 8 * size;
1447 if (chunkSize > 14)
1448 fseek (fp, (long) (chunkSize - 16), SEEK_CUR);
1449 else
1450 fseek (fp, (long) (chunkSize - 14), SEEK_CUR);
1451 }
1452 if ((format != WAVE_FORMAT_PCM) && (format != WAVE_FORMAT_IEEE_FLOAT)) {
1453 fprintf (stderr, "Error: Unsupported WAVE format in file '%s'.\n", src -> mPath);
1454 return (0);
1455 }
1456 if (src -> mChannel >= channels) {
1457 fprintf (stderr, "Error: Missing source channel in WAVE file '%s'.\n", src -> mPath);
1458 return (0);
1459 }
1460 if (rate != hrirRate) {
1461 fprintf (stderr, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src -> mPath);
1462 return (0);
1463 }
1464 if (format == WAVE_FORMAT_PCM) {
1465 if ((size < 2) || (size > 4)) {
1466 fprintf (stderr, "Error: Unsupported sample size in WAVE file '%s'.\n", src -> mPath);
1467 return (0);
1468 }
1469 if ((bits < 16) || (bits > (8 * size))) {
1470 fprintf (stderr, "Error: Bad significant bits in WAVE file '%s'.\n", src -> mPath);
1471 return (0);
1472 }
1473 src -> mType = ET_INT;
1474 } else {
1475 if ((size != 4) && (size != 8)) {
1476 fprintf (stderr, "Error: Unsupported sample size in WAVE file '%s'.\n", src -> mPath);
1477 return (0);
1478 }
1479 src -> mType = ET_FP;
1480 }
1481 src -> mSize = size;
1482 src -> mBits = (int) bits;
1483 src -> mSkip = channels;
1484 return (1);
1485 }
1486
1487 // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
1488 static int ReadWaveData (FILE * fp, const SourceRefT * src, const ByteOrderT order, const uint n, double * hrir) {
1489 int pre, post, skip;
1490 uint i;
1491
1492 pre = (int) (src -> mSize * src -> mChannel);
1493 post = (int) (src -> mSize * (src -> mSkip - src -> mChannel - 1));
1494 skip = 0;
1495 for (i = 0; i < n; i ++) {
1496 skip += pre;
1497 if (skip > 0)
1498 fseek (fp, skip, SEEK_CUR);
1499 if (! ReadBinAsDouble (fp, src -> mPath, order, src -> mType, src -> mSize, src -> mBits, & hrir [i]))
1500 return (0);
1501 skip = post;
1502 }
1503 if (skip > 0)
1504 fseek (fp, skip, SEEK_CUR);
1505 return (1);
1506 }
1507
1508 // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
1509 // doubles.
1510 static int ReadWaveList (FILE * fp, const SourceRefT * src, const ByteOrderT order, const uint n, double * hrir) {
1511 uint4 fourCC, chunkSize, listSize, count;
1512 uint block, skip, offset, i;
1513 double lastSample;
1514
1515 for (;;) {
1516 if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
1517 (! ReadBin4 (fp, src -> mPath, order, 4, & chunkSize)))
1518 return (0);
1519 if (fourCC == FOURCC_DATA) {
1520 block = src -> mSize * src -> mSkip;
1521 count = chunkSize / block;
1522 if (count < (src -> mOffset + n)) {
1523 fprintf (stderr, "Error: Bad read from file '%s'.\n", src -> mPath);
1524 return (0);
1525 }
1526 fseek (fp, (long) (src -> mOffset * block), SEEK_CUR);
1527 if (! ReadWaveData (fp, src, order, n, & hrir [0]))
1528 return (0);
1529 return (1);
1530 } else if (fourCC == FOURCC_LIST) {
1531 if (! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC))
1532 return (0);
1533 chunkSize -= 4;
1534 if (fourCC == FOURCC_WAVL)
1535 break;
1536 }
1537 if (chunkSize > 0)
1538 fseek (fp, (long) chunkSize, SEEK_CUR);
1539 }
1540 listSize = chunkSize;
1541 block = src -> mSize * src -> mSkip;
1542 skip = src -> mOffset;
1543 offset = 0;
1544 lastSample = 0.0;
1545 while ((offset < n) && (listSize > 8)) {
1546 if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
1547 (! ReadBin4 (fp, src -> mPath, order, 4, & chunkSize)))
1548 return (0);
1549 listSize -= 8 + chunkSize;
1550 if (fourCC == FOURCC_DATA) {
1551 count = chunkSize / block;
1552 if (count > skip) {
1553 fseek (fp, (long) (skip * block), SEEK_CUR);
1554 chunkSize -= skip * block;
1555 count -= skip;
1556 skip = 0;
1557 if (count > (n - offset))
1558 count = n - offset;
1559 if (! ReadWaveData (fp, src, order, count, & hrir [offset]))
1560 return (0);
1561 chunkSize -= count * block;
1562 offset += count;
1563 lastSample = hrir [offset - 1];
1564 } else {
1565 skip -= count;
1566 count = 0;
1567 }
1568 } else if (fourCC == FOURCC_SLNT) {
1569 if (! ReadBin4 (fp, src -> mPath, order, 4, & count))
1570 return (0);
1571 chunkSize -= 4;
1572 if (count > skip) {
1573 count -= skip;
1574 skip = 0;
1575 if (count > (n - offset))
1576 count = n - offset;
1577 for (i = 0; i < count; i ++)
1578 hrir [offset + i] = lastSample;
1579 offset += count;
1580 } else {
1581 skip -= count;
1582 count = 0;
1583 }
1584 }
1585 if (chunkSize > 0)
1586 fseek (fp, (long) chunkSize, SEEK_CUR);
1587 }
1588 if (offset < n) {
1589 fprintf (stderr, "Error: Bad read from file '%s'.\n", src -> mPath);
1590 return (0);
1591 }
1592 return (1);
1593 }
1594
1595 // Load a source HRIR from a RIFF/RIFX WAVE file.
1596 static int LoadWaveSource (FILE * fp, SourceRefT * src, const uint hrirRate, const uint n, double * hrir) {
1597 uint4 fourCC, dummy;
1598 ByteOrderT order;
1599
1600 if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
1601 (! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & dummy)))
1602 return (0);
1603 if (fourCC == FOURCC_RIFF) {
1604 order = BO_LITTLE;
1605 } else if (fourCC == FOURCC_RIFX) {
1606 order = BO_BIG;
1607 } else {
1608 fprintf (stderr, "Error: No RIFF/RIFX chunk in file '%s'.\n", src -> mPath);
1609 return (0);
1610 }
1611 if (! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC))
1612 return (0);
1613 if (fourCC != FOURCC_WAVE) {
1614 fprintf (stderr, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src -> mPath);
1615 return (0);
1616 }
1617 if (! ReadWaveFormat (fp, order, hrirRate, src))
1618 return (0);
1619 if (! ReadWaveList (fp, src, order, n, hrir))
1620 return (0);
1621 return (1);
1622 }
1623
1624 // Load a source HRIR from a binary file.
1625 static int LoadBinarySource (FILE * fp, const SourceRefT * src, const ByteOrderT order, const uint n, double * hrir) {
1626 uint i;
1627
1628 fseek (fp, (long) src -> mOffset, SEEK_SET);
1629 for (i = 0; i < n; i ++) {
1630 if (! ReadBinAsDouble (fp, src -> mPath, order, src -> mType, src -> mSize, src -> mBits, & hrir [i]))
1631 return (0);
1632 if (src -> mSkip > 0)
1633 fseek (fp, (long) src -> mSkip, SEEK_CUR);
1634 }
1635 return (1);
1636 }
1637
1638 // Load a source HRIR from an ASCII text file containing a list of elements
1639 // separated by whitespace or common list operators (',', ';', ':', '|').
1640 static int LoadAsciiSource (FILE * fp, const SourceRefT * src, const uint n, double * hrir) {
1641 TokenReaderT tr;
1642 uint i, j;
1643 double dummy;
1644
1645 TrSetup (fp, NULL, & tr);
1646 for (i = 0; i < src -> mOffset; i ++) {
1647 if (! ReadAsciiAsDouble (& tr, src -> mPath, src -> mType, (uint) src -> mBits, & dummy))
1648 return (0);
1649 }
1650 for (i = 0; i < n; i ++) {
1651 if (! ReadAsciiAsDouble (& tr, src -> mPath, src -> mType, (uint) src -> mBits, & hrir [i]))
1652 return (0);
1653 for (j = 0; j < src -> mSkip; j ++) {
1654 if (! ReadAsciiAsDouble (& tr, src -> mPath, src -> mType, (uint) src -> mBits, & dummy))
1655 return (0);
1656 }
1657 }
1658 return (1);
1659 }
1660
1661 // Load a source HRIR from a supported file type.
1662 static int LoadSource (SourceRefT * src, const uint hrirRate, const uint n, double * hrir) {
1663 FILE * fp = NULL;
1664 int result;
1665
1666 if (src -> mFormat == SF_ASCII)
1667 fp = fopen (src -> mPath, "r");
1668 else
1669 fp = fopen (src -> mPath, "rb");
1670 if (fp == NULL) {
1671 fprintf (stderr, "Error: Could not open source file '%s'.\n", src -> mPath);
1672 return (0);
1673 }
1674 if (src -> mFormat == SF_WAVE)
1675 result = LoadWaveSource (fp, src, hrirRate, n, hrir);
1676 else if (src -> mFormat == SF_BIN_LE)
1677 result = LoadBinarySource (fp, src, BO_LITTLE, n, hrir);
1678 else if (src -> mFormat == SF_BIN_BE)
1679 result = LoadBinarySource (fp, src, BO_BIG, n, hrir);
1680 else
1681 result = LoadAsciiSource (fp, src, n, hrir);
1682 fclose (fp);
1683 return (result);
1684 }
1685
1686 // Calculate the magnitude response of an HRIR and average it with any
1687 // existing responses for its elevation and azimuth.
1688 static void AverageHrirMagnitude (const double * hrir, const double f, const uint ei, const uint ai, const HrirDataT * hData) {
1689 double * re = NULL, * im = NULL;
1690 uint n, m, i, j;
1691
1692 n = hData -> mFftSize;
1693 re = CreateArray (n);
1694 im = CreateArray (n);
1695 for (i = 0; i < hData -> mIrPoints; i ++) {
1696 re [i] = hrir [i];
1697 im [i] = 0.0;
1698 }
1699 for (; i < n; i ++) {
1700 re [i] = 0.0;
1701 im [i] = 0.0;
1702 }
1703 FftForward (n, re, im, re, im);
1704 MagnitudeResponse (n, re, im, re);
1705 m = 1 + (n / 2);
1706 j = (hData -> mEvOffset [ei] + ai) * hData -> mIrSize;
1707 for (i = 0; i < m; i ++)
1708 hData -> mHrirs [j + i] = Lerp (hData -> mHrirs [j + i], re [i], f);
1709 DestroyArray (im);
1710 DestroyArray (re);
1711 }
1712
1713 /* Calculate the contribution of each HRIR to the diffuse-field average based
1714 * on the area of its surface patch. All patches are centered at the HRIR
1715 * coordinates on the unit sphere and are measured by solid angle.
1716 */
1717 static void CalculateDfWeights (const HrirDataT * hData, double * weights) {
1718 uint ei;
1719 double evs, sum, ev, up_ev, down_ev, solidAngle;
1720
1721 evs = 90.0 / (hData -> mEvCount - 1);
1722 sum = 0.0;
1723 for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++) {
1724 // For each elevation, calculate the upper and lower limits of the
1725 // patch band.
1726 ev = -90.0 + (ei * 2.0 * evs);
1727 if (ei < (hData -> mEvCount - 1))
1728 up_ev = (ev + evs) * M_PI / 180.0;
1729 else
1730 up_ev = M_PI / 2.0;
1731 if (ei > 0)
1732 down_ev = (ev - evs) * M_PI / 180.0;
1733 else
1734 down_ev = -M_PI / 2.0;
1735 // Calculate the area of the patch band.
1736 solidAngle = 2.0 * M_PI * (sin (up_ev) - sin (down_ev));
1737 // Each weight is the area of one patch.
1738 weights [ei] = solidAngle / hData -> mAzCount [ei];
1739 // Sum the total surface area covered by the HRIRs.
1740 sum += solidAngle;
1741 }
1742 // Normalize the weights given the total surface coverage.
1743 for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
1744 weights [ei] /= sum;
1745 }
1746
1747 /* Calculate the diffuse-field average from the given magnitude responses of
1748 * the HRIR set. Weighting can be applied to compensate for the varying
1749 * surface area covered by each HRIR. The final average can then be limited
1750 * by the specified magnitude range (in positive dB; 0.0 to skip).
1751 */
1752 static void CalculateDiffuseFieldAverage (const HrirDataT * hData, const int weighted, const double limit, double * dfa) {
1753 double * weights = NULL;
1754 uint ei, ai, count, step, start, end, m, j, i;
1755 double weight;
1756
1757 weights = CreateArray (hData -> mEvCount);
1758 if (weighted) {
1759 // Use coverage weighting to calculate the average.
1760 CalculateDfWeights (hData, weights);
1761 } else {
1762 // If coverage weighting is not used, the weights still need to be
1763 // averaged by the number of HRIRs.
1764 count = 0;
1765 for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
1766 count += hData -> mAzCount [ei];
1767 for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
1768 weights [ei] = 1.0 / count;
1769 }
1770 ei = hData -> mEvStart;
1771 ai = 0;
1772 step = hData -> mIrSize;
1773 start = hData -> mEvOffset [ei] * step;
1774 end = hData -> mIrCount * step;
1775 m = 1 + (hData -> mFftSize / 2);
1776 for (i = 0; i < m; i ++)
1777 dfa [i] = 0.0;
1778 for (j = start; j < end; j += step) {
1779 // Get the weight for this HRIR's contribution.
1780 weight = weights [ei];
1781 // Add this HRIR's weighted power average to the total.
1782 for (i = 0; i < m; i ++)
1783 dfa [i] += weight * hData -> mHrirs [j + i] * hData -> mHrirs [j + i];
1784 // Determine the next weight to use.
1785 ai ++;
1786 if (ai >= hData -> mAzCount [ei]) {
1787 ei ++;
1788 ai = 0;
1789 }
1790 }
1791 // Finish the average calculation and keep it from being too small.
1792 for (i = 0; i < m; i ++)
1793 dfa [i] = fmax (sqrt (dfa [i]), EPSILON);
1794 // Apply a limit to the magnitude range of the diffuse-field average if
1795 // desired.
1796 if (limit > 0.0)
1797 LimitMagnitudeResponse (hData -> mFftSize, limit, dfa, dfa);
1798 DestroyArray (weights);
1799 }
1800
1801 // Perform diffuse-field equalization on the magnitude responses of the HRIR
1802 // set using the given average response.
1803 static void DiffuseFieldEqualize (const double * dfa, const HrirDataT * hData) {
1804 uint step, start, end, m, j, i;
1805
1806 step = hData -> mIrSize;
1807 start = hData -> mEvOffset [hData -> mEvStart] * step;
1808 end = hData -> mIrCount * step;
1809 m = 1 + (hData -> mFftSize / 2);
1810 for (j = start; j < end; j += step) {
1811 for (i = 0; i < m; i ++)
1812 hData -> mHrirs [j + i] /= dfa [i];
1813 }
1814 }
1815
1816 // Perform minimum-phase reconstruction using the magnitude responses of the
1817 // HRIR set.
1818 static void ReconstructHrirs (const HrirDataT * hData) {
1819 double * re = NULL, * im = NULL;
1820 uint step, start, end, n, j, i;
1821
1822 step = hData -> mIrSize;
1823 start = hData -> mEvOffset [hData -> mEvStart] * step;
1824 end = hData -> mIrCount * step;
1825 n = hData -> mFftSize;
1826 re = CreateArray (n);
1827 im = CreateArray (n);
1828 for (j = start; j < end; j += step) {
1829 MinimumPhase (n, & hData -> mHrirs [j], re, im);
1830 FftInverse (n, re, im, re, im);
1831 for (i = 0; i < hData -> mIrPoints; i ++)
1832 hData -> mHrirs [j + i] = re [i];
1833 }
1834 DestroyArray (im);
1835 DestroyArray (re);
1836 }
1837
1838 // Resamples the HRIRs for use at the given sampling rate.
1839 static void ResampleHrirs (const uint rate, HrirDataT * hData) {
1840 ResamplerT rs;
1841 uint n, step, start, end, j;
1842
1843 ResamplerSetup (& rs, hData -> mIrRate, rate);
1844 n = hData -> mIrPoints;
1845 step = hData -> mIrSize;
1846 start = hData -> mEvOffset [hData -> mEvStart] * step;
1847 end = hData -> mIrCount * step;
1848 for (j = start; j < end; j += step)
1849 ResamplerRun (& rs, n, & hData -> mHrirs [j], n, & hData -> mHrirs [j]);
1850 ResamplerClear (& rs);
1851 hData -> mIrRate = rate;
1852 }
1853
1854 /* Given an elevation index and an azimuth, calculate the indices of the two
1855 * HRIRs that bound the coordinate along with a factor for calculating the
1856 * continous HRIR using interpolation.
1857 */
1858 static void CalcAzIndices (const HrirDataT * hData, const uint ei, const double az, uint * j0, uint * j1, double * jf) {
1859 double af;
1860 uint ai;
1861
1862 af = ((2.0 * M_PI) + az) * hData -> mAzCount [ei] / (2.0 * M_PI);
1863 ai = ((uint) af) % hData -> mAzCount [ei];
1864 af -= floor (af);
1865 (* j0) = hData -> mEvOffset [ei] + ai;
1866 (* j1) = hData -> mEvOffset [ei] + ((ai + 1) % hData -> mAzCount [ei]);
1867 (* jf) = af;
1868 }
1869
1870 /* Attempt to synthesize any missing HRIRs at the bottom elevations. Right
1871 * now this just blends the lowest elevation HRIRs together and applies some
1872 * attenuation and high frequency damping. It is a simple, if inaccurate
1873 * model.
1874 */
1875 static void SynthesizeHrirs (HrirDataT * hData) {
1876 uint oi, a, e, step, n, i, j;
1877 double of, b;
1878 uint j0, j1;
1879 double jf;
1880 double lp [4], s0, s1;
1881
1882 if (hData -> mEvStart <= 0)
1883 return;
1884 step = hData -> mIrSize;
1885 oi = hData -> mEvStart;
1886 n = hData -> mIrPoints;
1887 for (i = 0; i < n; i ++)
1888 hData -> mHrirs [i] = 0.0;
1889 for (a = 0; a < hData -> mAzCount [oi]; a ++) {
1890 j = (hData -> mEvOffset [oi] + a) * step;
1891 for (i = 0; i < n; i ++)
1892 hData -> mHrirs [i] += hData -> mHrirs [j + i] / hData -> mAzCount [oi];
1893 }
1894 for (e = 1; e < hData -> mEvStart; e ++) {
1895 of = ((double) e) / hData -> mEvStart;
1896 b = (1.0 - of) * (3.5e-6 * hData -> mIrRate);
1897 for (a = 0; a < hData -> mAzCount [e]; a ++) {
1898 j = (hData -> mEvOffset [e] + a) * step;
1899 CalcAzIndices (hData, oi, a * 2.0 * M_PI / hData -> mAzCount [e], & j0, & j1, & jf);
1900 j0 *= step;
1901 j1 *= step;
1902 lp [0] = 0.0;
1903 lp [1] = 0.0;
1904 lp [2] = 0.0;
1905 lp [3] = 0.0;
1906 for (i = 0; i < n; i ++) {
1907 s0 = hData -> mHrirs [i];
1908 s1 = Lerp (hData -> mHrirs [j0 + i], hData -> mHrirs [j1 + i], jf);
1909 s0 = Lerp (s0, s1, of);
1910 lp [0] = Lerp (s0, lp [0], b);
1911 lp [1] = Lerp (lp [0], lp [1], b);
1912 lp [2] = Lerp (lp [1], lp [2], b);
1913 lp [3] = Lerp (lp [2], lp [3], b);
1914 hData -> mHrirs [j + i] = lp [3];
1915 }
1916 }
1917 }
1918 b = 3.5e-6 * hData -> mIrRate;
1919 lp [0] = 0.0;
1920 lp [1] = 0.0;
1921 lp [2] = 0.0;
1922 lp [3] = 0.0;
1923 for (i = 0; i < n; i ++) {
1924 s0 = hData -> mHrirs [i];
1925 lp [0] = Lerp (s0, lp [0], b);
1926 lp [1] = Lerp (lp [0], lp [1], b);
1927 lp [2] = Lerp (lp [1], lp [2], b);
1928 lp [3] = Lerp (lp [2], lp [3], b);
1929 hData -> mHrirs [i] = lp [3];
1930 }
1931 hData -> mEvStart = 0;
1932 }
1933
1934 // The following routines assume a full set of HRIRs for all elevations.
1935
1936 // Normalize the HRIR set and slightly attenuate the result.
1937 static void NormalizeHrirs (const HrirDataT * hData) {
1938 uint step, end, n, j, i;
1939 double maxLevel;
1940
1941 step = hData -> mIrSize;
1942 end = hData -> mIrCount * step;
1943 n = hData -> mIrPoints;
1944 maxLevel = 0.0;
1945 for (j = 0; j < end; j += step) {
1946 for (i = 0; i < n; i ++)
1947 maxLevel = fmax (fabs (hData -> mHrirs [j + i]), maxLevel);
1948 }
1949 maxLevel = 1.01 * maxLevel;
1950 for (j = 0; j < end; j += step) {
1951 for (i = 0; i < n; i ++)
1952 hData -> mHrirs [j + i] /= maxLevel;
1953 }
1954 }
1955
1956 // Calculate the left-ear time delay using a spherical head model.
1957 static double CalcLTD (const double ev, const double az, const double rad, const double dist) {
1958 double azp, dlp, l, al;
1959
1960 azp = asin (cos (ev) * sin (az));
1961 dlp = sqrt ((dist * dist) + (rad * rad) + (2.0 * dist * rad * sin (azp)));
1962 l = sqrt ((dist * dist) - (rad * rad));
1963 al = (0.5 * M_PI) + azp;
1964 if (dlp > l)
1965 dlp = l + (rad * (al - acos (rad / dist)));
1966 return (dlp / 343.3);
1967 }
1968
1969 // Calculate the effective head-related time delays for the each HRIR, now
1970 // that they are minimum-phase.
1971 static void CalculateHrtds (HrirDataT * hData) {
1972 double minHrtd, maxHrtd;
1973 uint e, a, j;
1974 double t;
1975
1976 minHrtd = 1000.0;
1977 maxHrtd = -1000.0;
1978 for (e = 0; e < hData -> mEvCount; e ++) {
1979 for (a = 0; a < hData -> mAzCount [e]; a ++) {
1980 j = hData -> mEvOffset [e] + a;
1981 t = CalcLTD ((-90.0 + (e * 180.0 / (hData -> mEvCount - 1))) * M_PI / 180.0,
1982 (a * 360.0 / hData -> mAzCount [e]) * M_PI / 180.0,
1983 hData -> mRadius, hData -> mDistance);
1984 hData -> mHrtds [j] = t;
1985 maxHrtd = fmax (t, maxHrtd);
1986 minHrtd = fmin (t, minHrtd);
1987 }
1988 }
1989 maxHrtd -= minHrtd;
1990 for (j = 0; j < hData -> mIrCount; j ++)
1991 hData -> mHrtds [j] -= minHrtd;
1992 hData -> mMaxHrtd = maxHrtd;
1993 }
1994
1995 // Store the OpenAL Soft HRTF data set.
1996 static int StoreMhr (const HrirDataT * hData, const char * filename) {
1997 FILE * fp = NULL;
1998 uint e, step, end, n, j, i;
1999 int hpHist, v;
2000
2001 if ((fp = fopen (filename, "wb")) == NULL) {
2002 fprintf (stderr, "Error: Could not open MHR file '%s'.\n", filename);
2003 return (0);
2004 }
2005 if (! WriteAscii (MHR_FORMAT, fp, filename))
2006 return (0);
2007 if (! WriteBin4 (BO_LITTLE, 4, (uint4) hData -> mIrRate, fp, filename))
2008 return (0);
2009 if (! WriteBin4 (BO_LITTLE, 1, (uint4) hData -> mIrPoints, fp, filename))
2010 return (0);
2011 if (! WriteBin4 (BO_LITTLE, 1, (uint4) hData -> mEvCount, fp, filename))
2012 return (0);
2013 for (e = 0; e < hData -> mEvCount; e ++) {
2014 if (! WriteBin4 (BO_LITTLE, 1, (uint4) hData -> mAzCount [e], fp, filename))
2015 return (0);
2016 }
2017 step = hData -> mIrSize;
2018 end = hData -> mIrCount * step;
2019 n = hData -> mIrPoints;
2020 srand (0x31DF840C);
2021 for (j = 0; j < end; j += step) {
2022 hpHist = 0;
2023 for (i = 0; i < n; i ++) {
2024 v = HpTpdfDither (32767.0 * hData -> mHrirs [j + i], & hpHist);
2025 if (! WriteBin4 (BO_LITTLE, 2, (uint4) v, fp, filename))
2026 return (0);
2027 }
2028 }
2029 for (j = 0; j < hData -> mIrCount; j ++) {
2030 v = (int) fmin (round (hData -> mIrRate * hData -> mHrtds [j]), MAX_HRTD);
2031 if (! WriteBin4 (BO_LITTLE, 1, (uint4) v, fp, filename))
2032 return (0);
2033 }
2034 fclose (fp);
2035 return (1);
2036 }
2037
2038 // Store the OpenAL Soft built-in table.
2039 static int StoreTable (const HrirDataT * hData, const char * filename) {
2040 FILE * fp = NULL;
2041 uint step, end, n, j, i;
2042 int hpHist, v;
2043 char text [128 + 1];
2044
2045 if ((fp = fopen (filename, "wb")) == NULL) {
2046 fprintf (stderr, "Error: Could not open table file '%s'.\n", filename);
2047 return (0);
2048 }
2049 snprintf (text, 128, "/* Elevation metrics */\n"
2050 "static const ALubyte defaultAzCount[%u] = { ", hData -> mEvCount);
2051 if (! WriteAscii (text, fp, filename))
2052 return (0);
2053 for (i = 0; i < hData -> mEvCount; i ++) {
2054 snprintf (text, 128, "%u, ", hData -> mAzCount [i]);
2055 if (! WriteAscii (text, fp, filename))
2056 return (0);
2057 }
2058 snprintf (text, 128, "};\n"
2059 "static const ALushort defaultEvOffset[%u] = { ", hData -> mEvCount);
2060 if (! WriteAscii (text, fp, filename))
2061 return (0);
2062 for (i = 0; i < hData -> mEvCount; i ++) {
2063 snprintf (text, 128, "%u, ", hData -> mEvOffset [i]);
2064 if (! WriteAscii (text, fp, filename))
2065 return (0);
2066 }
2067 step = hData -> mIrSize;
2068 end = hData -> mIrCount * step;
2069 n = hData -> mIrPoints;
2070 snprintf (text, 128, "};\n\n"
2071 "/* HRIR Coefficients */\n"
2072 "static const ALshort defaultCoeffs[%u] =\n{\n", hData -> mIrCount * n);
2073 if (! WriteAscii (text, fp, filename))
2074 return (0);
2075 srand (0x31DF840C);
2076 for (j = 0; j < end; j += step) {
2077 if (! WriteAscii (" ", fp, filename))
2078 return (0);
2079 hpHist = 0;
2080 for (i = 0; i < n; i ++) {
2081 v = HpTpdfDither (32767.0 * hData -> mHrirs [j + i], & hpHist);
2082 snprintf (text, 128, " %+d,", v);
2083 if (! WriteAscii (text, fp, filename))
2084 return (0);
2085 }
2086 if (! WriteAscii ("\n", fp, filename))
2087 return (0);
2088 }
2089 snprintf (text, 128, "};\n\n"
2090 "/* HRIR Delays */\n"
2091 "static const ALubyte defaultDelays[%u] =\n{\n"
2092 " ", hData -> mIrCount);
2093 if (! WriteAscii (text, fp, filename))
2094 return (0);
2095 for (j = 0; j < hData -> mIrCount; j ++) {
2096 v = (int) fmin (round (hData -> mIrRate * hData -> mHrtds [j]), MAX_HRTD);
2097 snprintf (text, 128, " %d,", v);
2098 if (! WriteAscii (text, fp, filename))
2099 return (0);
2100 }
2101 if (! WriteAscii ("\n};\n\n"
2102 "/* Default HRTF Definition */\n", fp, filename))
2103 return (0);
2104 snprintf (text, 128, "static const struct Hrtf DefaultHrtf = {\n"
2105 " %u, %u, %u, defaultAzCount, defaultEvOffset,\n",
2106 hData -> mIrRate, hData -> mIrPoints, hData -> mEvCount);
2107 if (! WriteAscii (text, fp, filename))
2108 return (0);
2109 if (! WriteAscii (" defaultCoeffs, defaultDelays, NULL\n"
2110 "};\n", fp, filename))
2111 return (0);
2112 fclose (fp);
2113 return (1);
2114 }
2115
2116 // Process the data set definition to read and validate the data set metrics.
2117 static int ProcessMetrics (TokenReaderT * tr, const uint fftSize, const uint truncSize, HrirDataT * hData) {
2118 char ident [MAX_IDENT_LEN + 1];
2119 uint line, col;
2120 int intVal;
2121 uint points;
2122 double fpVal;
2123 int hasRate = 0, hasPoints = 0, hasAzimuths = 0;
2124 int hasRadius = 0, hasDistance = 0;
2125
2126 while (! (hasRate && hasPoints && hasAzimuths && hasRadius && hasDistance)) {
2127 TrIndication (tr, & line, & col);
2128 if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
2129 return (0);
2130 if (strcasecmp (ident, "rate") == 0) {
2131 if (hasRate) {
2132 TrErrorAt (tr, line, col, "Redefinition of 'rate'.\n");
2133 return (0);
2134 }
2135 if (! TrReadOperator (tr, "="))
2136 return (0);
2137 if (! TrReadInt (tr, MIN_RATE, MAX_RATE, & intVal))
2138 return (0);
2139 hData -> mIrRate = (uint) intVal;
2140 hasRate = 1;
2141 } else if (strcasecmp (ident, "points") == 0) {
2142 if (hasPoints) {
2143 TrErrorAt (tr, line, col, "Redefinition of 'points'.\n");
2144 return (0);
2145 }
2146 if (! TrReadOperator (tr, "="))
2147 return (0);
2148 TrIndication (tr, & line, & col);
2149 if (! TrReadInt (tr, MIN_POINTS, MAX_POINTS, & intVal))
2150 return (0);
2151 points = (uint) intVal;
2152 if ((fftSize > 0) && (points > fftSize)) {
2153 TrErrorAt (tr, line, col, "Value exceeds the overriden FFT size.\n");
2154 return (0);
2155 }
2156 if (points < truncSize) {
2157 TrErrorAt (tr, line, col, "Value is below the truncation size.\n");
2158 return (0);
2159 }
2160 hData -> mIrPoints = points;
2161 hData -> mFftSize = fftSize;
2162 if (fftSize <= 0) {
2163 points = 1;
2164 while (points < (4 * hData -> mIrPoints))
2165 points <<= 1;
2166 hData -> mFftSize = points;
2167 hData -> mIrSize = 1 + (points / 2);
2168 } else {
2169 hData -> mFftSize = fftSize;
2170 hData -> mIrSize = 1 + (fftSize / 2);
2171 if (points > hData -> mIrSize)
2172 hData -> mIrSize = points;
2173 }
2174 hasPoints = 1;
2175 } else if (strcasecmp (ident, "azimuths") == 0) {
2176 if (hasAzimuths) {
2177 TrErrorAt (tr, line, col, "Redefinition of 'azimuths'.\n");
2178 return (0);
2179 }
2180 if (! TrReadOperator (tr, "="))
2181 return (0);
2182 hData -> mIrCount = 0;
2183 hData -> mEvCount = 0;
2184 hData -> mEvOffset [0] = 0;
2185 for (;;) {
2186 if (! TrReadInt (tr, MIN_AZ_COUNT, MAX_AZ_COUNT, & intVal))
2187 return (0);
2188 hData -> mAzCount [hData -> mEvCount] = (uint) intVal;
2189 hData -> mIrCount += (uint) intVal;
2190 hData -> mEvCount ++;
2191 if (! TrIsOperator (tr, ","))
2192 break;
2193 if (hData -> mEvCount >= MAX_EV_COUNT) {
2194 TrError (tr, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT);
2195 return (0);
2196 }
2197 hData -> mEvOffset [hData -> mEvCount] = hData -> mEvOffset [hData -> mEvCount - 1] + ((uint) intVal);
2198 TrReadOperator (tr, ",");
2199 }
2200 if (hData -> mEvCount < MIN_EV_COUNT) {
2201 TrErrorAt (tr, line, col, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT);
2202 return (0);
2203 }
2204 hasAzimuths = 1;
2205 } else if (strcasecmp (ident, "radius") == 0) {
2206 if (hasRadius) {
2207 TrErrorAt (tr, line, col, "Redefinition of 'radius'.\n");
2208 return (0);
2209 }
2210 if (! TrReadOperator (tr, "="))
2211 return (0);
2212 if (! TrReadFloat (tr, MIN_RADIUS, MAX_RADIUS, & fpVal))
2213 return (0);
2214 hData -> mRadius = fpVal;
2215 hasRadius = 1;
2216 } else if (strcasecmp (ident, "distance") == 0) {
2217 if (hasDistance) {
2218 TrErrorAt (tr, line, col, "Redefinition of 'distance'.\n");
2219 return (0);
2220 }
2221 if (! TrReadOperator (tr, "="))
2222 return (0);
2223 if (! TrReadFloat (tr, MIN_DISTANCE, MAX_DISTANCE, & fpVal))
2224 return (0);
2225 hData -> mDistance = fpVal;
2226 hasDistance = 1;
2227 } else {
2228 TrErrorAt (tr, line, col, "Expected a metric name.\n");
2229 return (0);
2230 }
2231 TrSkipWhitespace (tr);
2232 }
2233 return (1);
2234 }
2235
2236 // Parse an index pair from the data set definition.
2237 static int ReadIndexPair (TokenReaderT * tr, const HrirDataT * hData, uint * ei, uint * ai) {
2238 int intVal;
2239
2240 if (! TrReadInt (tr, 0, (int) hData -> mEvCount, & intVal))
2241 return (0);
2242 (* ei) = (uint) intVal;
2243 if (! TrReadOperator (tr, ","))
2244 return (0);
2245 if (! TrReadInt (tr, 0, (int) hData -> mAzCount [(* ei)], & intVal))
2246 return (0);
2247 (* ai) = (uint) intVal;
2248 return (1);
2249 }
2250
2251 // Match the source format from a given identifier.
2252 static SourceFormatT MatchSourceFormat (const char * ident) {
2253 if (strcasecmp (ident, "wave") == 0)
2254 return (SF_WAVE);
2255 else if (strcasecmp (ident, "bin_le") == 0)
2256 return (SF_BIN_LE);
2257 else if (strcasecmp (ident, "bin_be") == 0)
2258 return (SF_BIN_BE);
2259 else if (strcasecmp (ident, "ascii") == 0)
2260 return (SF_ASCII);
2261 return (SF_NONE);
2262 }
2263
2264 // Match the source element type from a given identifier.
2265 static ElementTypeT MatchElementType (const char * ident) {
2266 if (strcasecmp (ident, "int") == 0)
2267 return (ET_INT);
2268 else if (strcasecmp (ident, "fp") == 0)
2269 return (ET_FP);
2270 return (ET_NONE);
2271 }
2272
2273 // Parse and validate a source reference from the data set definition.
2274 static int ReadSourceRef (TokenReaderT * tr, SourceRefT * src) {
2275 uint line, col;
2276 char ident [MAX_IDENT_LEN + 1];
2277 int intVal;
2278
2279 TrIndication (tr, & line, & col);
2280 if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
2281 return (0);
2282 src -> mFormat = MatchSourceFormat (ident);
2283 if (src -> mFormat == SF_NONE) {
2284 TrErrorAt (tr, line, col, "Expected a source format.\n");
2285 return (0);
2286 }
2287 if (! TrReadOperator (tr, "("))
2288 return (0);
2289 if (src -> mFormat == SF_WAVE) {
2290 if (! TrReadInt (tr, 0, MAX_WAVE_CHANNELS, & intVal))
2291 return (0);
2292 src -> mType = ET_NONE;
2293 src -> mSize = 0;
2294 src -> mBits = 0;
2295 src -> mChannel = (uint) intVal;
2296 src -> mSkip = 0;
2297 } else {
2298 TrIndication (tr, & line, & col);
2299 if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
2300 return (0);
2301 src -> mType = MatchElementType (ident);
2302 if (src -> mType == ET_NONE) {
2303 TrErrorAt (tr, line, col, "Expected a source element type.\n");
2304 return (0);
2305 }
2306 if ((src -> mFormat == SF_BIN_LE) || (src -> mFormat == SF_BIN_BE)) {
2307 if (! TrReadOperator (tr, ","))
2308 return (0);
2309 if (src -> mType == ET_INT) {
2310 if (! TrReadInt (tr, MIN_BIN_SIZE, MAX_BIN_SIZE, & intVal))
2311 return (0);
2312 src -> mSize = (uint) intVal;
2313 if (TrIsOperator (tr, ",")) {
2314 TrReadOperator (tr, ",");
2315 TrIndication (tr, & line, & col);
2316 if (! TrReadInt (tr, -2147483647 - 1, 2147483647, & intVal))
2317 return (0);
2318 if ((abs (intVal) < MIN_BIN_BITS) || (((uint) abs (intVal)) > (8 * src -> mSize))) {
2319 TrErrorAt (tr, line, col, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS, 8 * src -> mSize);
2320 return (0);
2321 }
2322 src -> mBits = intVal;
2323 } else {
2324 src -> mBits = (int) (8 * src -> mSize);
2325 }
2326 } else {
2327 TrIndication (tr, & line, & col);
2328 if (! TrReadInt (tr, -2147483647 - 1, 2147483647, & intVal))
2329 return (0);
2330 if ((intVal != 4) && (intVal != 8)) {
2331 TrErrorAt (tr, line, col, "Expected a value of 4 or 8.\n");
2332 return (0);
2333 }
2334 src -> mSize = (uint) intVal;
2335 src -> mBits = 0;
2336 }
2337 } else if ((src -> mFormat == SF_ASCII) && (src -> mType == ET_INT)) {
2338 if (! TrReadOperator (tr, ","))
2339 return (0);
2340 if (! TrReadInt (tr, MIN_ASCII_BITS, MAX_ASCII_BITS, & intVal))
2341 return (0);
2342 src -> mSize = 0;
2343 src -> mBits = intVal;
2344 } else {
2345 src -> mSize = 0;
2346 src -> mBits = 0;
2347 }
2348 if (TrIsOperator (tr, ";")) {
2349 TrReadOperator (tr, ";");
2350 if (! TrReadInt (tr, 0, 0x7FFFFFFF, & intVal))
2351 return (0);
2352 src -> mSkip = (uint) intVal;
2353 } else {
2354 src -> mSkip = 0;
2355 }
2356 }
2357 if (! TrReadOperator (tr, ")"))
2358 return (0);
2359 if (TrIsOperator (tr, "@")) {
2360 TrReadOperator (tr, "@");
2361 if (! TrReadInt (tr, 0, 0x7FFFFFFF, & intVal))
2362 return (0);
2363 src -> mOffset = (uint) intVal;
2364 } else {
2365 src -> mOffset = 0;
2366 }
2367 if (! TrReadOperator (tr, ":"))
2368 return (0);
2369 if (! TrReadString (tr, MAX_PATH_LEN, src -> mPath))
2370 return (0);
2371 return (1);
2372 }
2373
2374 // Process the list of sources in the data set definition.
2375 static int ProcessSources (TokenReaderT * tr, HrirDataT * hData) {
2376 uint * setCount = NULL, * setFlag = NULL;
2377 double * hrir = NULL;
2378 uint line, col, ei, ai;
2379 SourceRefT src;
2380 double factor;
2381
2382 setCount = (uint *) calloc (hData -> mEvCount, sizeof (uint));
2383 setFlag = (uint *) calloc (hData -> mIrCount, sizeof (uint));
2384 hrir = CreateArray (hData -> mIrPoints);
2385 while (TrIsOperator (tr, "[")) {
2386 TrIndication (tr, & line, & col);
2387 TrReadOperator (tr, "[");
2388 if (ReadIndexPair (tr, hData, & ei, & ai)) {
2389 if (TrReadOperator (tr, "]")) {
2390 if (! setFlag [hData -> mEvOffset [ei] + ai]) {
2391 if (TrReadOperator (tr, "=")) {
2392 factor = 1.0;
2393 for (;;) {
2394 if (ReadSourceRef (tr, & src)) {
2395 if (LoadSource (& src, hData -> mIrRate, hData -> mIrPoints, hrir)) {
2396 AverageHrirMagnitude (hrir, 1.0 / factor, ei, ai, hData);
2397 factor += 1.0;
2398 if (! TrIsOperator (tr, "+"))
2399 break;
2400 TrReadOperator (tr, "+");
2401 continue;
2402 }
2403 }
2404 DestroyArray (hrir);
2405 free (setFlag);
2406 free (setCount);
2407 return (0);
2408 }
2409 setFlag [hData -> mEvOffset [ei] + ai] = 1;
2410 setCount [ei] ++;
2411 continue;
2412 }
2413 } else {
2414 TrErrorAt (tr, line, col, "Redefinition of source.\n");
2415 }
2416 }
2417 }
2418 DestroyArray (hrir);
2419 free (setFlag);
2420 free (setCount);
2421 return (0);
2422 }
2423 ei = 0;
2424 while ((ei < hData -> mEvCount) && (setCount [ei] < 1))
2425 ei ++;
2426 if (ei < hData -> mEvCount) {
2427 hData -> mEvStart = ei;
2428 while ((ei < hData -> mEvCount) && (setCount [ei] == hData -> mAzCount [ei]))
2429 ei ++;
2430 if (ei >= hData -> mEvCount) {
2431 if (! TrLoad (tr)) {
2432 DestroyArray (hrir);
2433 free (setFlag);
2434 free (setCount);
2435 return (1);
2436 } else {
2437 TrError (tr, "Errant data at end of source list.\n");
2438 }
2439 } else {
2440 TrError (tr, "Missing sources for elevation index %d.\n", ei);
2441 }
2442 } else {
2443 TrError (tr, "Missing source references.\n");
2444 }
2445 DestroyArray (hrir);
2446 free (setFlag);
2447 free (setCount);
2448 return (0);
2449 }
2450
2451 /* Parse the data set definition and process the source data, storing the
2452 * resulting data set as desired. If the input name is NULL it will read
2453 * from standard input.
2454 */
2455 static int ProcessDefinition (const char * inName, const uint outRate, const uint fftSize, const int equalize, const int surface, const double limit, const uint truncSize, const OutputFormatT outFormat, const char * outName) {
2456 FILE * fp = NULL;
2457 TokenReaderT tr;
2458 HrirDataT hData;
2459 double * dfa = NULL;
2460 char rateStr [8 + 1], expName [MAX_PATH_LEN];
2461
2462 hData . mIrRate = 0;
2463 hData . mIrPoints = 0;
2464 hData . mFftSize = 0;
2465 hData . mIrSize = 0;
2466 hData . mIrCount = 0;
2467 hData . mEvCount = 0;
2468 hData . mRadius = 0;
2469 hData . mDistance = 0;
2470 fprintf (stdout, "Reading HRIR definition...\n");
2471 if (inName != NULL) {
2472 fp = fopen (inName, "r");
2473 if (fp == NULL) {
2474 fprintf (stderr, "Error: Could not open definition file '%s'\n", inName);
2475 return (0);
2476 }
2477 TrSetup (fp, inName, & tr);
2478 } else {
2479 fp = stdin;
2480 TrSetup (fp, "<stdin>", & tr);
2481 }
2482 if (! ProcessMetrics (& tr, fftSize, truncSize, & hData)) {
2483 if (inName != NULL)
2484 fclose (fp);
2485 return (0);
2486 }
2487 hData . mHrirs = CreateArray (hData . mIrCount * hData . mIrSize);
2488 hData . mHrtds = CreateArray (hData . mIrCount);
2489 if (! ProcessSources (& tr, & hData)) {
2490 DestroyArray (hData . mHrtds);
2491 DestroyArray (hData . mHrirs);
2492 if (inName != NULL)
2493 fclose (fp);
2494 return (0);
2495 }
2496 if (inName != NULL)
2497 fclose (fp);
2498 if (equalize) {
2499 dfa = CreateArray (1 + (hData . mFftSize / 2));
2500 fprintf (stdout, "Calculating diffuse-field average...\n");
2501 CalculateDiffuseFieldAverage (& hData, surface, limit, dfa);
2502 fprintf (stdout, "Performing diffuse-field equalization...\n");
2503 DiffuseFieldEqualize (dfa, & hData);
2504 DestroyArray (dfa);
2505 }
2506 fprintf (stdout, "Performing minimum phase reconstruction...\n");
2507 ReconstructHrirs (& hData);
2508 if ((outRate != 0) && (outRate != hData . mIrRate)) {
2509 fprintf (stdout, "Resampling HRIRs...\n");
2510 ResampleHrirs (outRate, & hData);
2511 }
2512 fprintf (stdout, "Truncating minimum-phase HRIRs...\n");
2513 hData . mIrPoints = truncSize;
2514 fprintf (stdout, "Synthesizing missing elevations...\n");
2515 SynthesizeHrirs (& hData);
2516 fprintf (stdout, "Normalizing final HRIRs...\n");
2517 NormalizeHrirs (& hData);
2518 fprintf (stdout, "Calculating impulse delays...\n");
2519 CalculateHrtds (& hData);
2520 snprintf (rateStr, 8, "%u", hData . mIrRate);
2521 StrSubst (outName, "%r", rateStr, MAX_PATH_LEN, expName);
2522 switch (outFormat) {
2523 case OF_MHR :
2524 fprintf (stdout, "Creating MHR data set file...\n");
2525 if (! StoreMhr (& hData, expName))
2526 return (0);
2527 break;
2528 case OF_TABLE :
2529 fprintf (stderr, "Creating OpenAL Soft table file...\n");
2530 if (! StoreTable (& hData, expName))
2531 return (0);
2532 break;
2533 default :
2534 break;
2535 }
2536 DestroyArray (hData . mHrtds);
2537 DestroyArray (hData . mHrirs);
2538 return (1);
2539 }
2540
2541 // Standard command line dispatch.
2542 int main (const int argc, const char * argv []) {
2543 const char * inName = NULL, * outName = NULL;
2544 OutputFormatT outFormat;
2545 int argi;
2546 uint outRate, fftSize;
2547 int equalize, surface;
2548 double limit;
2549 uint truncSize;
2550 char * end = NULL;
2551
2552 if (argc < 2) {
2553 fprintf (stderr, "Error: No command specified. See '%s -h' for help.\n", argv [0]);
2554 return (-1);
2555 }
2556 if ((strcmp (argv [1], "--help") == 0) || (strcmp (argv [1], "-h") == 0)) {
2557 fprintf (stdout, "HRTF Processing and Composition Utility\n\n");
2558 fprintf (stdout, "Usage: %s <command> [<option>...]\n\n", argv [0]);
2559 fprintf (stdout, "Commands:\n");
2560 fprintf (stdout, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n");
2561 fprintf (stdout, " Defaults output to: ./oalsoft_hrtf_%%r.mhr\n");
2562 fprintf (stdout, " -t, --make-tab Makes the built-in table used when compiling OpenAL Soft.\n");
2563 fprintf (stdout, " Defaults output to: ./hrtf_tables.inc\n");
2564 fprintf (stdout, " -h, --help Displays this help information.\n\n");
2565 fprintf (stdout, "Options:\n");
2566 fprintf (stdout, " -r=<rate> Change the data set sample rate to the specified value and\n");
2567 fprintf (stdout, " resample the HRIRs accordingly.\n");
2568 fprintf (stdout, " -f=<points> Override the FFT window size (defaults to the first power-\n");
2569 fprintf (stdout, " of-two that fits four times the number of HRIR points).\n");
2570 fprintf (stdout, " -e={on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE ? "on" : "off"));
2571 fprintf (stdout, " -s={on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE ? "on" : "off"));
2572 fprintf (stdout, " -l={<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
2573 fprintf (stdout, " average (default: %.2f).\n", DEFAULT_LIMIT);
2574 fprintf (stdout, " -w=<points> Specify the size of the truncation window that's applied\n");
2575 fprintf (stdout, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE);
2576 fprintf (stdout, " -i=<filename> Specify an HRIR definition file to use (defaults to stdin).\n");
2577 fprintf (stdout, " -o=<filename> Specify an output file. Overrides command-selected default.\n");
2578 fprintf (stdout, " Use of '%%r' will be substituted with the data set sample rate.\n");
2579 return (0);
2580 }
2581 if ((strcmp (argv [1], "--make-mhr") == 0) || (strcmp (argv [1], "-m") == 0)) {
2582 if (argc > 3)
2583 outName = argv [3];
2584 else
2585 outName = "./oalsoft_hrtf_%r.mhr";
2586 outFormat = OF_MHR;
2587 } else if ((strcmp (argv [1], "--make-tab") == 0) || (strcmp (argv [1], "-t") == 0)) {
2588 if (argc > 3)
2589 outName = argv [3];
2590 else
2591 outName = "./hrtf_tables.inc";
2592 outFormat = OF_TABLE;
2593 } else {
2594 fprintf (stderr, "Error: Invalid command '%s'.\n", argv [1]);
2595 return (-1);
2596 }
2597 argi = 2;
2598 outRate = 0;
2599 fftSize = 0;
2600 equalize = DEFAULT_EQUALIZE;
2601 surface = DEFAULT_SURFACE;
2602 limit = DEFAULT_LIMIT;
2603 truncSize = DEFAULT_TRUNCSIZE;
2604 while (argi < argc) {
2605 if (strncmp (argv [argi], "-r=", 3) == 0) {
2606 outRate = strtoul (& argv [argi] [3], & end, 10);
2607 if ((end [0] != '\0') || (outRate < MIN_RATE) || (outRate > MAX_RATE)) {
2608 fprintf (stderr, "Error: Expected a value from %u to %u for '-r'.\n", MIN_RATE, MAX_RATE);
2609 return (-1);
2610 }
2611 } else if (strncmp (argv [argi], "-f=", 3) == 0) {
2612 fftSize = strtoul (& argv [argi] [3], & end, 10);
2613 if ((end [0] != '\0') || (fftSize & (fftSize - 1)) || (fftSize < MIN_FFTSIZE) || (fftSize > MAX_FFTSIZE)) {
2614 fprintf (stderr, "Error: Expected a power-of-two value from %u to %u for '-f'.\n", MIN_FFTSIZE, MAX_FFTSIZE);
2615 return (-1);
2616 }
2617 } else if (strncmp (argv [argi], "-e=", 3) == 0) {
2618 if (strcmp (& argv [argi] [3], "on") == 0) {
2619 equalize = 1;
2620 } else if (strcmp (& argv [argi] [3], "off") == 0) {
2621 equalize = 0;
2622 } else {
2623 fprintf (stderr, "Error: Expected 'on' or 'off' for '-e'.\n");
2624 return (-1);
2625 }
2626 } else if (strncmp (argv [argi], "-s=", 3) == 0) {
2627 if (strcmp (& argv [argi] [3], "on") == 0) {
2628 surface = 1;
2629 } else if (strcmp (& argv [argi] [3], "off") == 0) {
2630 surface = 0;
2631 } else {
2632 fprintf (stderr, "Error: Expected 'on' or 'off' for '-s'.\n");
2633 return (-1);
2634 }
2635 } else if (strncmp (argv [argi], "-l=", 3) == 0) {
2636 if (strcmp (& argv [argi] [3], "none") == 0) {
2637 limit = 0.0;
2638 } else {
2639 limit = strtod (& argv [argi] [3], & end);
2640 if ((end [0] != '\0') || (limit < MIN_LIMIT) || (limit > MAX_LIMIT)) {
2641 fprintf (stderr, "Error: Expected 'none' or a value from %.2f to %.2f for '-l'.\n", MIN_LIMIT, MAX_LIMIT);
2642 return (-1);
2643 }
2644 }
2645 } else if (strncmp (argv [argi], "-w=", 3) == 0) {
2646 truncSize = strtoul (& argv [argi] [3], & end, 10);
2647 if ((end [0] != '\0') || (truncSize < MIN_TRUNCSIZE) || (truncSize > MAX_TRUNCSIZE) || (truncSize % MOD_TRUNCSIZE)) {
2648 fprintf (stderr, "Error: Expected a value from %u to %u in multiples of %u for '-w'.\n", MIN_TRUNCSIZE, MAX_TRUNCSIZE, MOD_TRUNCSIZE);
2649 return (-1);
2650 }
2651 } else if (strncmp (argv [argi], "-i=", 3) == 0) {
2652 inName = & argv [argi] [3];
2653 } else if (strncmp (argv [argi], "-o=", 3) == 0) {
2654 outName = & argv [argi] [3];
2655 } else {
2656 fprintf (stderr, "Error: Invalid option '%s'.\n", argv [argi]);
2657 return (-1);
2658 }
2659 argi ++;
2660 }
2661 if (! ProcessDefinition (inName, outRate, fftSize, equalize, surface, limit, truncSize, outFormat, outName))
2662 return (-1);
2663 fprintf (stdout, "Operation completed.\n");
2664 return (0);
2665 }
2666
2667
2668
2669
StoreTable
static int StoreTable(const HrirDataT *hData, const char *filename)
Definition: makehrtf.c:2039
MIN_POINTS
#define MIN_POINTS
Definition: makehrtf.c:106
ReadWaveData
static int ReadWaveData(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
Definition: makehrtf.c:1488
pi
static const double pi
Definition: e_atan2.c:47
ALubyte
unsigned char ALubyte
Definition: al.h:47
OutputFormatT
OutputFormatT
Definition: makehrtf.c:212
MIN_DISTANCE
#define MIN_DISTANCE
Definition: makehrtf.c:123
TrSkipWhitespace
static int TrSkipWhitespace(TokenReaderT *tr)
Definition: makehrtf.c:399
TrSetup
static void TrSetup(FILE *fp, const char *filename, TokenReaderT *tr)
Definition: makehrtf.c:308
val
GLuint const GLfloat * val
Definition: glew.h:2715
AverageHrirMagnitude
static void AverageHrirMagnitude(const double *hrir, const double f, const uint ei, const uint ai, const HrirDataT *hData)
Definition: makehrtf.c:1688
int4
ALint int4
Definition: makehrtf.c:224
s
GLdouble s
Definition: glew.h:1376
MIN_RADIUS
#define MIN_RADIUS
Definition: makehrtf.c:118
FOURCC_WAVL
#define FOURCC_WAVL
Definition: makehrtf.c:173
MAX_AZ_COUNT
#define MAX_AZ_COUNT
Definition: makehrtf.c:115
SynthesizeHrirs
static void SynthesizeHrirs(HrirDataT *hData)
Definition: makehrtf.c:1875
FftForward
static void FftForward(const uint n, const double *inR, const double *inI, double *outR, double *outI)
Definition: makehrtf.c:888
TrIndication
static void TrIndication(TokenReaderT *tr, uint *line, uint *column)
Definition: makehrtf.c:422
type
GLint GLenum GLsizei GLsizei GLsizei GLint GLenum GLenum type
Definition: gl2ext.h:845
DEFAULT_TRUNCSIZE
#define DEFAULT_TRUNCSIZE
Definition: makehrtf.c:164
CalcLTD
static double CalcLTD(const double ev, const double az, const double rad, const double dist)
Definition: makehrtf.c:1957
main
int main(int argc, char **argv)
Definition: bootstrap.cpp:102
MAX_DISTANCE
#define MAX_DISTANCE
Definition: makehrtf.c:124
ComplexAbs
static double ComplexAbs(const double r, const double i)
Definition: makehrtf.c:787
NULL
#define NULL
Definition: ftobjs.h:61
DiffuseFieldEqualize
static void DiffuseFieldEqualize(const double *dfa, const HrirDataT *hData)
Definition: makehrtf.c:1803
start
GLuint start
Definition: glew.h:1239
MIN_BIN_SIZE
#define MIN_BIN_SIZE
Definition: makehrtf.c:132
f
GLclampf f
Definition: glew.h:3390
MIN_TRUNCSIZE
#define MIN_TRUNCSIZE
Definition: makehrtf.c:153
MAX_ASCII_BITS
#define MAX_ASCII_BITS
Definition: makehrtf.c:142
MagnitudeResponse
static void MagnitudeResponse(const uint n, const double *inR, const double *inI, double *out)
Definition: makehrtf.c:949
LimitMagnitudeResponse
static void LimitMagnitudeResponse(const uint n, const double limit, const double *in, double *out)
Definition: makehrtf.c:961
k
int32_t k
Definition: e_log.c:102
CalcAzIndices
static void CalcAzIndices(const HrirDataT *hData, const uint ei, const double az, uint *j0, uint *j1, double *jf)
Definition: makehrtf.c:1858
n
GLclampd n
Definition: glew.h:7287
FOURCC_RIFX
#define FOURCC_RIFX
Definition: makehrtf.c:168
MIN_RATE
#define MIN_RATE
Definition: makehrtf.c:102
x
EGLSurface EGLint x
Definition: eglext.h:293
free
SDL_EventEntry * free
Definition: SDL_events.c:80
surface
EGLSurface surface
Definition: eglext.h:74
e
return Display return Display Bool Bool int int e
Definition: SDL_x11sym.h:30
uint64_t
unsigned long long uint64_t
Definition: SDL_config_windows.h:64
t
GLdouble GLdouble t
Definition: glew.h:1384
in
GLuint in
Definition: glew.h:10672
ALint
int ALint
Definition: al.h:56
j
int32_t j
Definition: e_log.c:102
Kaiser
static double Kaiser(const double b, const double k)
Definition: makehrtf.c:1069
ReadAsciiAsDouble
static int ReadAsciiAsDouble(TokenReaderT *tr, const char *filename, const ElementTypeT type, const uint bits, double *out)
Definition: makehrtf.c:1379
SincFilter
static double SincFilter(const int l, const double b, const double gain, const double cutoff, const int i)
Definition: makehrtf.c:1132
a
GLboolean GLboolean GLboolean GLboolean a
Definition: glew.h:8736
BesselI_0
static double BesselI_0(const double x)
Definition: makehrtf.c:1034
FftInverse
static void FftInverse(const uint n, const double *inR, const double *inI, double *outR, double *outI)
Definition: makehrtf.c:894
name
EGLImageKHR EGLint * name
Definition: eglext.h:284
MatchElementType
static ElementTypeT MatchElementType(const char *ident)
Definition: makehrtf.c:2265
CalculateHrtds
static void CalculateHrtds(HrirDataT *hData)
Definition: makehrtf.c:1971
SourceFormatT
SourceFormatT
Definition: makehrtf.c:196
ProcessSources
static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
Definition: makehrtf.c:2375
len
GLenum GLsizei len
Definition: glew.h:7035
StrSubst
static int StrSubst(const char *in, const char *pat, const char *rep, const size_t maxLen, char *out)
Definition: makehrtf.c:694
weight
GLuint GLuint GLfloat weight
Definition: glew.h:12401
bits
#define bits
Definition: infblock.c:15
MAX_IDENT_LEN
#define MAX_IDENT_LEN
Definition: makehrtf.c:95
MHR_FORMAT
#define MHR_FORMAT
Definition: makehrtf.c:186
s0
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s0
Definition: glew.h:11582
ResamplerSetup
static void ResamplerSetup(ResamplerT *rs, const uint srcRate, const uint dstRate)
Definition: makehrtf.c:1164
MIN_FFTSIZE
#define MIN_FFTSIZE
Definition: makehrtf.c:145
calloc
#define calloc
Definition: SDL_malloc.c:636
TrErrorAt
static void TrErrorAt(const TokenReaderT *tr, uint line, uint column, const char *format,...)
Definition: makehrtf.c:365
LoadWaveSource
static int LoadWaveSource(FILE *fp, SourceRefT *src, const uint hrirRate, const uint n, double *hrir)
Definition: makehrtf.c:1596
BO_LITTLE
Definition: makehrtf.c:191
SF_BIN_BE
Definition: makehrtf.c:200
Lerp
static double Lerp(const double a, const double b, const double f)
Definition: makehrtf.c:750
ProcessDefinition
static int ProcessDefinition(const char *inName, const uint outRate, const uint fftSize, const int equalize, const int surface, const double limit, const uint truncSize, const OutputFormatT outFormat, const char *outName)
Definition: makehrtf.c:2455
ResamplerRun
static void ResamplerRun(ResamplerT *rs, const uint inN, const double *in, const uint outN, double *out)
Definition: makehrtf.c:1201
OF_TABLE
Definition: makehrtf.c:215
TrReadInt
static int TrReadInt(TokenReaderT *tr, const int loBound, const int hiBound, int *value)
Definition: makehrtf.c:491
MAX_POINTS
#define MAX_POINTS
Definition: makehrtf.c:107
FOURCC_DATA
#define FOURCC_DATA
Definition: makehrtf.c:171
MAX_HRTD
#define MAX_HRTD
Definition: makehrtf.c:182
MinimumPhase
static void MinimumPhase(const uint n, const double *in, double *outR, double *outI)
Definition: makehrtf.c:991
ComplexMul
static void ComplexMul(const double aR, const double aI, const double bR, const double bI, double *outR, double *outI)
Definition: makehrtf.c:792
SEEK_CUR
#define SEEK_CUR
Definition: zconf.h:250
FOURCC_WAVE
#define FOURCC_WAVE
Definition: makehrtf.c:169
TokenReaderT
struct TokenReaderT TokenReaderT
Definition: makehrtf.c:244
WAVE_FORMAT_EXTENSIBLE
#define WAVE_FORMAT_EXTENSIBLE
Definition: makehrtf.c:179
CalculateDfWeights
static void CalculateDfWeights(const HrirDataT *hData, double *weights)
Definition: makehrtf.c:1717
FftArrange
static void FftArrange(const uint n, const double *inR, const double *inI, double *outR, double *outI)
Definition: makehrtf.c:812
ResampleHrirs
static void ResampleHrirs(const uint rate, HrirDataT *hData)
Definition: makehrtf.c:1839
uint
unsigned int uint
Definition: makehrtf.c:219
FOURCC_LIST
#define FOURCC_LIST
Definition: makehrtf.c:172
ReadWaveList
static int ReadWaveList(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
Definition: makehrtf.c:1510
ReadSourceRef
static int ReadSourceRef(TokenReaderT *tr, SourceRefT *src)
Definition: makehrtf.c:2274
MatchSourceFormat
static SourceFormatT MatchSourceFormat(const char *ident)
Definition: makehrtf.c:2252
uint4
ALuint uint4
Definition: makehrtf.c:225
result
GLuint64EXT * result
Definition: glew.h:12708
v
const GLdouble * v
Definition: glew.h:1377
int
int
Definition: SDL_systhread.c:37
Gcd
static uint Gcd(const uint a, const uint b)
Definition: makehrtf.c:1080
TrReadIdent
static int TrReadIdent(TokenReaderT *tr, const uint maxLen, char *ident)
Definition: makehrtf.c:458
width
EGLSurface EGLint EGLint EGLint width
Definition: eglext.h:293
sort_controllers.filename
string filename
Definition: sort_controllers.py:8
MAX_FFTSIZE
#define MAX_FFTSIZE
Definition: makehrtf.c:146
FOURCC_FMT
#define FOURCC_FMT
Definition: makehrtf.c:170
TrIsOperator
static int TrIsOperator(TokenReaderT *tr, const char *op)
Definition: makehrtf.c:432
CreateArray
static double * CreateArray(const size_t n)
Definition: makehrtf.c:768
EPSILON
#define EPSILON
Definition: makehrtf.c:83
MIN_EV_COUNT
#define MIN_EV_COUNT
Definition: makehrtf.c:110
count
GLint GLsizei count
Definition: gl2ext.h:1011
p
GLfloat GLfloat p
Definition: glew.h:14938
WAVE_FORMAT_IEEE_FLOAT
#define WAVE_FORMAT_IEEE_FLOAT
Definition: makehrtf.c:178
MAX_WAVE_CHANNELS
#define MAX_WAVE_CHANNELS
Definition: makehrtf.c:128
j1
const GLubyte GLuint GLuint GLuint GLuint alpha GLboolean GLboolean GLboolean GLboolean alpha GLint GLint GLsizei GLsizei GLenum type GLenum GLint GLenum GLint GLint GLsizei GLsizei GLint border GLenum GLint GLint GLint GLint GLint GLsizei GLsizei height GLsizei GLsizei GLenum GLenum const GLvoid *pixels GLenum GLint GLint GLint j1
Definition: SDL_glfuncs.h:118
MAX_EV_COUNT
#define MAX_EV_COUNT
Definition: makehrtf.c:111
BO_NONE
Definition: makehrtf.c:190
ReadBin8
static int ReadBin8(FILE *fp, const char *filename, const ByteOrderT order, uint8 *out)
Definition: makehrtf.c:1269
column
GLenum GLenum GLvoid GLvoid * column
Definition: glew.h:4447
ResamplerClear
static void ResamplerClear(ResamplerT *rs)
Definition: makehrtf.c:1194
bits
GLenum GLint GLenum GLsizei GLsizei GLsizei GLint GLsizei const GLvoid * bits
Definition: SDL_opengl.h:10449
NormalizeHrirs
static void NormalizeHrirs(const HrirDataT *hData)
Definition: makehrtf.c:1937
SF_ASCII
Definition: makehrtf.c:201
format
GLint GLenum GLsizei GLsizei GLsizei GLint GLenum format
Definition: gl2ext.h:845
limit
GLint limit
Definition: glew.h:11829
WriteBin4
static int WriteBin4(const ByteOrderT order, const uint bytes, const uint4 in, FILE *fp, const char *filename)
Definition: makehrtf.c:1310
MIN_AZ_COUNT
#define MIN_AZ_COUNT
Definition: makehrtf.c:114
ALuint
unsigned int ALuint
Definition: al.h:59
MAX_RATE
#define MAX_RATE
Definition: makehrtf.c:103
floor
double floor(double x)
Definition: s_floor.c:40
ReconstructHrirs
static void ReconstructHrirs(const HrirDataT *hData)
Definition: makehrtf.c:1818
s1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
Definition: glew.h:11582
Clamp
static double Clamp(const double val, const double lower, const double upper)
Definition: makehrtf.c:745
ReadBinAsDouble
static int ReadBinAsDouble(FILE *fp, const char *filename, const ByteOrderT order, const ElementTypeT type, const uint bytes, const int bits, double *out)
Definition: makehrtf.c:1339
MAX_BIN_SIZE
#define MAX_BIN_SIZE
Definition: makehrtf.c:133
ComplexExp
static void ComplexExp(const double inR, const double inI, double *outR, double *outI)
Definition: makehrtf.c:798
TrReadFloat
static int TrReadFloat(TokenReaderT *tr, const double loBound, const double hiBound, double *value)
Definition: makehrtf.c:536
WAVE_FORMAT_PCM
#define WAVE_FORMAT_PCM
Definition: makehrtf.c:177
LoadSource
static int LoadSource(SourceRefT *src, const uint hrirRate, const uint n, double *hrir)
Definition: makehrtf.c:1662
SEEK_SET
#define SEEK_SET
Definition: zconf.h:249
sin
double sin(double x)
Definition: s_sin.c:56
FOURCC_SLNT
#define FOURCC_SLNT
Definition: makehrtf.c:174
OF_MHR
Definition: makehrtf.c:214
TrSkipLine
static void TrSkipLine(TokenReaderT *tr)
Definition: makehrtf.c:383
LoadAsciiSource
static int LoadAsciiSource(FILE *fp, const SourceRefT *src, const uint n, double *hrir)
Definition: makehrtf.c:1640
SourceRefT
struct SourceRefT SourceRefT
Definition: makehrtf.c:245
TrReadOperator
static int TrReadOperator(TokenReaderT *tr, const char *op)
Definition: makehrtf.c:667
y
EGLSurface EGLint EGLint y
Definition: eglext.h:293
MAX_PATH_LEN
#define MAX_PATH_LEN
Definition: makehrtf.c:98
TR_RING_SIZE
#define TR_RING_SIZE
Definition: makehrtf.c:87
l
GLdouble l
Definition: glew.h:8383
M_PI
#define M_PI
Definition: makehrtf.c:75
StoreMhr
static int StoreMhr(const HrirDataT *hData, const char *filename)
Definition: makehrtf.c:1996
Sinc
static double Sinc(const double x)
Definition: makehrtf.c:1022
ByteOrderT
ByteOrderT
Definition: makehrtf.c:189
value
EGLSurface EGLint void ** value
Definition: eglext.h:301
ET_NONE
Definition: makehrtf.c:206
ElementTypeT
ElementTypeT
Definition: makehrtf.c:205
offset
GLintptr offset
Definition: glew.h:1668
OF_NONE
Definition: makehrtf.c:213
ProcessMetrics
static int ProcessMetrics(TokenReaderT *tr, const uint fftSize, const uint truncSize, HrirDataT *hData)
Definition: makehrtf.c:2117
ET_FP
Definition: makehrtf.c:208
DEFAULT_EQUALIZE
#define DEFAULT_EQUALIZE
Definition: makehrtf.c:161
ReadBin4
static int ReadBin4(FILE *fp, const char *filename, const ByteOrderT order, const uint bytes, uint4 *out)
Definition: makehrtf.c:1241
MAX_RADIUS
#define MAX_RADIUS
Definition: makehrtf.c:119
MOD_TRUNCSIZE
#define MOD_TRUNCSIZE
Definition: makehrtf.c:158
q
GLdouble GLdouble GLdouble GLdouble q
Definition: glew.h:1400
SF_BIN_LE
Definition: makehrtf.c:199
CalcKaiserOrder
static uint CalcKaiserOrder(const double rejection, const double transition)
Definition: makehrtf.c:1100
TR_LOAD_SIZE
#define TR_LOAD_SIZE
Definition: makehrtf.c:91
Hilbert
static void Hilbert(const uint n, const double *in, double *outR, double *outI)
Definition: makehrtf.c:912
SF_WAVE
Definition: makehrtf.c:198
ReadIndexPair
static int ReadIndexPair(TokenReaderT *tr, const HrirDataT *hData, uint *ei, uint *ai)
Definition: makehrtf.c:2237
x2
GLfixed GLfixed x2
Definition: glext.h:4559
String.h
r
GLdouble GLdouble GLdouble r
Definition: glew.h:1392
SF_NONE
Definition: makehrtf.c:197
TrError
static void TrError(const TokenReaderT *tr, const char *format,...)
Definition: makehrtf.c:374
b
GLdouble GLdouble GLdouble b
Definition: glew.h:8383
end
GLuint GLuint end
Definition: glew.h:1239
z
GLint GLint GLint GLint z
Definition: gl2ext.h:1214
TrErrorVA
static void TrErrorVA(const TokenReaderT *tr, uint line, uint column, const char *format, va_list argPtr)
Definition: makehrtf.c:357
MIN_LIMIT
#define MIN_LIMIT
Definition: makehrtf.c:149
order
GLuint GLdouble GLdouble GLint GLint order
Definition: glew.h:3337
FOURCC_RIFF
#define FOURCC_RIFF
Definition: makehrtf.c:167
skip
png_uint_32 skip
Definition: pngrutil.c:1241
uint1
ALubyte uint1
Definition: makehrtf.c:222
BO_BIG
Definition: makehrtf.c:192
points
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
Definition: glew.h:3337
MAX_LIMIT
#define MAX_LIMIT
Definition: makehrtf.c:150
MIN_ASCII_BITS
#define MIN_ASCII_BITS
Definition: makehrtf.c:141
src
GLenum src
Definition: glew.h:2396
ET_INT
Definition: makehrtf.c:207
DestroyArray
static void DestroyArray(const double *a)
Definition: makehrtf.c:780
DEFAULT_SURFACE
#define DEFAULT_SURFACE
Definition: makehrtf.c:162
ResamplerT
struct ResamplerT ResamplerT
Definition: makehrtf.c:247
i
int i
Definition: pngrutil.c:1377
FftSummation
static void FftSummation(const uint n, const double s, double *re, double *im)
Definition: makehrtf.c:848
ReadWaveFormat
static int ReadWaveFormat(FILE *fp, const ByteOrderT order, const uint hrirRate, SourceRefT *src)
Definition: makehrtf.c:1407
MAX_TRUNCSIZE
#define MAX_TRUNCSIZE
Definition: makehrtf.c:154
cos
double cos(double x)
Definition: s_cos.c:56
WriteAscii
static int WriteAscii(const char *out, FILE *fp, const char *filename)
Definition: makehrtf.c:1296
CalculateDiffuseFieldAverage
static void CalculateDiffuseFieldAverage(const HrirDataT *hData, const int weighted, const double limit, double *dfa)
Definition: makehrtf.c:1752
HpTpdfDither
static int HpTpdfDither(const double in, int *hpHist)
Definition: makehrtf.c:756
weights
GLbyte * weights
Definition: glew.h:6834
TrReadString
static int TrReadString(TokenReaderT *tr, const uint maxLen, char *text)
Definition: makehrtf.c:624
HrirDataT
struct HrirDataT HrirDataT
Definition: makehrtf.c:246
fabs
double fabs(double x)
Definition: s_fabs.c:29
m
#define m(i, j)
TrLoad
static int TrLoad(TokenReaderT *tr)
Definition: makehrtf.c:331
TR_RING_MASK
#define TR_RING_MASK
Definition: makehrtf.c:88
MIN_BIN_BITS
#define MIN_BIN_BITS
Definition: makehrtf.c:137
CalcKaiserBeta
static double CalcKaiserBeta(const double rejection)
Definition: makehrtf.c:1110
DEFAULT_LIMIT
#define DEFAULT_LIMIT
Definition: makehrtf.c:163
term
#define term(a, b, c, d)
LoadBinarySource
static int LoadBinarySource(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
Definition: makehrtf.c:1625
size
GLsizei size
Definition: gl2ext.h:1467
m
const GLdouble * m
Definition: glew.h:8385