2 * HRTF utility for producing and demonstrating the process of creating an
3 * OpenAL Soft compatible HRIR data set.
5 * Copyright (C) 2011-2012 Christopher Fitzgerald
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.
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.
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.
21 * Or visit: http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
23 * --------------------------------------------------------------------------
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.
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:
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
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
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.
44 * The spherical head algorithm for calculating propagation delay was adapted
47 * Modeling Interaural Time Difference Assuming a Spherical Head
49 * Music 150, Musical Acoustics, Stanford University
61 // Needed for 64-bit unsigned integer.
64 // Rely (if naively) on OpenAL's header for the types used for serialization.
68 #define M_PI 3.14159265358979323846
71 // The epsilon used to maintain signal stability.
72 #define EPSILON (1e-15)
74 // Constants for accessing the token reader's ring buffer.
75 #define TR_RING_BITS (16)
76 #define TR_RING_SIZE (1 << TR_RING_BITS)
77 #define TR_RING_MASK (TR_RING_SIZE - 1)
79 // The token reader's load interval in bytes.
80 #define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
82 // The maximum identifier length used when processing the data set
84 #define MAX_IDENT_LEN (16)
86 // The maximum path length used when processing filenames.
87 #define MAX_PATH_LEN (256)
89 // The limits for the sample 'rate' metric in the data set definition.
90 #define MIN_RATE (32000)
91 #define MAX_RATE (96000)
93 // The limits for the HRIR 'points' metric in the data set definition.
94 #define MIN_POINTS (16)
95 #define MAX_POINTS (8192)
97 // The limits to the number of 'azimuths' listed in the data set definition.
98 #define MIN_EV_COUNT (5)
99 #define MAX_EV_COUNT (128)
101 // The limits for each of the 'azimuths' listed in the data set definition.
102 #define MIN_AZ_COUNT (1)
103 #define MAX_AZ_COUNT (128)
105 // The limits for the listener's head 'radius' in the data set definition.
106 #define MIN_RADIUS (0.05)
107 #define MAX_RADIUS (0.15)
109 // The limits for the 'distance' from source to listener in the definition
111 #define MIN_DISTANCE (0.5)
112 #define MAX_DISTANCE (2.5)
114 // The maximum number of channels that can be addressed for a WAVE file
115 // source listed in the data set definition.
116 #define MAX_WAVE_CHANNELS (65535)
118 // The limits to the byte size for a binary source listed in the definition
120 #define MIN_BIN_SIZE (2)
121 #define MAX_BIN_SIZE (4)
123 // The minimum number of significant bits for binary sources listed in the
124 // data set definition. The maximum is calculated from the byte size.
125 #define MIN_BIN_BITS (16)
127 // The limits to the number of significant bits for an ASCII source listed in
128 // the data set definition.
129 #define MIN_ASCII_BITS (16)
130 #define MAX_ASCII_BITS (32)
132 // The limits to the FFT window size override on the command line.
133 #define MIN_FFTSIZE (512)
134 #define MAX_FFTSIZE (16384)
136 // The limits to the equalization range limit on the command line.
137 #define MIN_LIMIT (2.0)
138 #define MAX_LIMIT (120.0)
140 // The limits to the truncation window size on the command line.
141 #define MIN_TRUNCSIZE (8)
142 #define MAX_TRUNCSIZE (128)
144 // The truncation window size must be a multiple of the below value to allow
145 // for vectorized convolution.
146 #define MOD_TRUNCSIZE (8)
148 // The defaults for the command line options.
149 #define DEFAULT_EQUALIZE (1)
150 #define DEFAULT_SURFACE (1)
151 #define DEFAULT_LIMIT (24.0)
152 #define DEFAULT_TRUNCSIZE (32)
154 // The four-character-codes for RIFF/RIFX WAVE file chunks.
155 #define FOURCC_RIFF (0x46464952) // 'RIFF'
156 #define FOURCC_RIFX (0x58464952) // 'RIFX'
157 #define FOURCC_WAVE (0x45564157) // 'WAVE'
158 #define FOURCC_FMT (0x20746D66) // 'fmt '
159 #define FOURCC_DATA (0x61746164) // 'data'
160 #define FOURCC_LIST (0x5453494C) // 'LIST'
161 #define FOURCC_WAVL (0x6C766177) // 'wavl'
162 #define FOURCC_SLNT (0x746E6C73) // 'slnt'
164 // The supported wave formats.
165 #define WAVE_FORMAT_PCM (0x0001)
166 #define WAVE_FORMAT_IEEE_FLOAT (0x0003)
167 #define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
169 // The maximum propagation delay value supported by OpenAL Soft.
170 #define MAX_HRTD (63.0)
172 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
173 // response protocol 01.
174 #define MHR_FORMAT ("MinPHR01")
176 // Byte order for the serialization routines.
183 // Source format for the references listed in the data set definition.
186 SF_WAVE
, // RIFF/RIFX WAVE file.
187 SF_BIN_LE
, // Little-endian binary file.
188 SF_BIN_BE
, // Big-endian binary file.
189 SF_ASCII
// ASCII text file.
192 // Element types for the references listed in the data set definition.
195 ET_INT
, // Integer elements.
196 ET_FP
// Floating-point elements.
199 // Desired output format from the command line.
202 OF_MHR
, // OpenAL Soft MHR data set file.
203 OF_TABLE
// OpenAL Soft built-in table file (used when compiling).
206 // Unsigned integer type.
207 typedef unsigned int uint
;
209 // Serialization types. The trailing digit indicates the number of bytes.
210 typedef ALubyte uint1
;
213 typedef ALuint uint4
;
215 #if defined (HAVE_STDINT_H)
218 typedef uint64_t uint8
;
219 #elif defined (HAVE___INT64)
220 typedef unsigned __int64 uint8
;
221 #elif (SIZEOF_LONG == 8)
222 typedef unsigned long uint8
;
223 #elif (SIZEOF_LONG_LONG == 8)
224 typedef unsigned long long uint8
;
227 typedef enum ByteOrderT ByteOrderT
;
228 typedef enum SourceFormatT SourceFormatT
;
229 typedef enum ElementTypeT ElementTypeT
;
230 typedef enum OutputFormatT OutputFormatT
;
232 typedef struct TokenReaderT TokenReaderT
;
233 typedef struct SourceRefT SourceRefT
;
234 typedef struct HrirDataT HrirDataT
;
236 // Token reader state for parsing the data set definition.
237 struct TokenReaderT
{
242 char mRing
[TR_RING_SIZE
];
247 // Source reference state used when loading sources.
249 SourceFormatT mFormat
;
256 char mPath
[MAX_PATH_LEN
+ 1];
259 // The HRIR metrics and data set used when loading, processing, and storing
260 // the resulting HRTF.
269 mAzCount
[MAX_EV_COUNT
],
270 mEvOffset
[MAX_EV_COUNT
];
278 /* Token reader routines for parsing text files. Whitespace is not
279 * significant. It can process tokens as identifiers, numbers (integer and
280 * floating-point), strings, and operators. Strings must be encapsulated by
281 * double-quotes and cannot span multiple lines.
284 // Setup the reader on the given file. The filename can be NULL if no error
285 // output is desired.
286 static void TrSetup (FILE * fp
, const char * filename
, TokenReaderT
* tr
) {
287 const char * name
= NULL
;
292 // If a filename was given, store a pointer to the base name.
293 if (filename
!= NULL
) {
294 while ((ch
= (* filename
)) != '\0') {
295 if ((ch
== '/') || (ch
== '\\'))
307 // Prime the reader's ring buffer, and return a result indicating that there
308 // is text to process.
309 static int TrLoad (TokenReaderT
* tr
) {
310 size_t toLoad
, in
, count
;
312 toLoad
= TR_RING_SIZE
- (tr
-> mIn
- tr
-> mOut
);
313 if ((toLoad
>= TR_LOAD_SIZE
) && (! feof (tr
-> mFile
))) {
314 // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
315 toLoad
= TR_LOAD_SIZE
;
316 in
= tr
-> mIn
& TR_RING_MASK
;
317 count
= TR_RING_SIZE
- in
;
318 if (count
< toLoad
) {
319 tr
-> mIn
+= fread (& tr
-> mRing
[in
], 1, count
, tr
-> mFile
);
320 tr
-> mIn
+= fread (& tr
-> mRing
[0], 1, toLoad
- count
, tr
-> mFile
);
322 tr
-> mIn
+= fread (& tr
-> mRing
[in
], 1, toLoad
, tr
-> mFile
);
324 if (tr
-> mOut
>= TR_RING_SIZE
) {
325 tr
-> mOut
-= TR_RING_SIZE
;
326 tr
-> mIn
-= TR_RING_SIZE
;
329 if (tr
-> mIn
> tr
-> mOut
)
334 // Error display routine. Only displays when the base name is not NULL.
335 static void TrErrorVA (const TokenReaderT
* tr
, uint line
, uint column
, const char * format
, va_list argPtr
) {
336 if (tr
-> mName
!= NULL
) {
337 fprintf (stderr
, "Error (%s:%d:%d): ", tr
-> mName
, line
, column
);
338 vfprintf (stderr
, format
, argPtr
);
342 // Used to display an error at a saved line/column.
343 static void TrErrorAt (const TokenReaderT
* tr
, uint line
, uint column
, const char * format
, ...) {
346 va_start (argPtr
, format
);
347 TrErrorVA (tr
, line
, column
, format
, argPtr
);
351 // Used to display an error at the current line/column.
352 static void TrError (const TokenReaderT
* tr
, const char * format
, ...) {
355 va_start (argPtr
, format
);
356 TrErrorVA (tr
, tr
-> mLine
, tr
-> mColumn
, format
, argPtr
);
360 // Skips to the next line.
361 static void TrSkipLine (TokenReaderT
* tr
) {
364 while (TrLoad (tr
)) {
365 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
376 // Skips to the next token.
377 static int TrSkipWhitespace (TokenReaderT
* tr
) {
380 while (TrLoad (tr
)) {
381 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
390 } else if (ch
== '#') {
399 // Get the line and/or column of the next token (or the end of input).
400 static void TrIndication (TokenReaderT
* tr
, uint
* line
, uint
* column
) {
401 TrSkipWhitespace (tr
);
403 (* line
) = tr
-> mLine
;
405 (* column
) = tr
-> mColumn
;
408 // Checks to see if a token is the given operator. It does not display any
409 // errors and will not proceed to the next token.
410 static int TrIsOperator (TokenReaderT
* tr
, const char * op
) {
414 if (! TrSkipWhitespace (tr
))
418 while ((op
[len
] != '\0') && (out
< tr
-> mIn
)) {
419 ch
= tr
-> mRing
[out
& TR_RING_MASK
];
425 if (op
[len
] == '\0')
430 /* The TrRead*() routines obtain the value of a matching token type. They
431 * display type, form, and boundary errors and will proceed to the next
435 // Reads and validates an identifier token.
436 static int TrReadIdent (TokenReaderT
* tr
, const size_t maxLen
, char * ident
) {
442 if (TrSkipWhitespace (tr
)) {
444 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
445 if ((ch
== '_') || isalpha (ch
)) {
454 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
455 } while ((ch
== '_') || isdigit (ch
) || isalpha (ch
));
456 tr
-> mColumn
+= len
;
458 TrErrorAt (tr
, tr
-> mLine
, col
, "Identifier is too long.\n");
465 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected an identifier.\n");
469 // Reads and validates (including bounds) an integer token.
470 static int TrReadInt (TokenReaderT
* tr
, const int loBound
, const int hiBound
, int * value
) {
473 char ch
, temp
[64 + 1];
476 if (TrSkipWhitespace (tr
)) {
479 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
480 if ((ch
== '+') || (ch
== '-')) {
486 while (TrLoad (tr
)) {
487 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
496 tr
-> mColumn
+= len
;
497 if ((digis
> 0) && (ch
!= '.') && (! isalpha (ch
))) {
499 TrErrorAt (tr
, tr
-> mLine
, col
, "Integer is too long.");
503 (* value
) = strtol (temp
, NULL
, 10);
504 if (((* value
) < loBound
) || ((* value
) > hiBound
)) {
505 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a value from %d to %d.\n", loBound
, hiBound
);
511 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected an integer.\n");
515 // Reads and validates (including bounds) a float token.
516 static int TrReadFloat (TokenReaderT
* tr
, const double loBound
, const double hiBound
, double * value
) {
519 char ch
, temp
[64 + 1];
522 if (TrSkipWhitespace (tr
)) {
525 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
526 if ((ch
== '+') || (ch
== '-')) {
532 while (TrLoad (tr
)) {
533 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
548 while (TrLoad (tr
)) {
549 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
559 if ((ch
== 'E') || (ch
== 'e')) {
565 if ((ch
== '+') || (ch
== '-')) {
571 while (TrLoad (tr
)) {
572 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
582 tr
-> mColumn
+= len
;
583 if ((digis
> 0) && (ch
!= '.') && (! isalpha (ch
))) {
585 TrErrorAt (tr
, tr
-> mLine
, col
, "Float is too long.");
589 (* value
) = strtod (temp
, NULL
);
590 if (((* value
) < loBound
) || ((* value
) > hiBound
)) {
591 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a value from %f to %f.\n", loBound
, hiBound
);
597 tr
-> mColumn
+= len
;
600 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a float.\n");
604 // Reads and validates a string token.
605 static int TrReadString (TokenReaderT
* tr
, const size_t maxLen
, char * text
) {
611 if (TrSkipWhitespace (tr
)) {
613 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
617 while (TrLoad (tr
)) {
618 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
623 TrErrorAt (tr
, tr
-> mLine
, col
, "Unterminated string at end of line.\n");
631 tr
-> mColumn
+= 1 + len
;
632 TrErrorAt (tr
, tr
-> mLine
, col
, "Unterminated string at end of input.\n");
635 tr
-> mColumn
+= 2 + len
;
637 TrErrorAt (tr
, tr
-> mLine
, col
, "String is too long.\n");
644 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a string.\n");
648 // Reads and validates the given operator.
649 static int TrReadOperator (TokenReaderT
* tr
, const char * op
) {
655 if (TrSkipWhitespace (tr
)) {
658 while ((op
[len
] != '\0') && TrLoad (tr
)) {
659 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
665 tr
-> mColumn
+= len
;
666 if (op
[len
] == '\0')
669 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected '%s' operator.\n", op
);
673 /* Performs a string substitution. Any case-insensitive occurrences of the
674 * pattern string are replaced with the replacement string. The result is
675 * truncated if necessary.
677 static int StrSubst (const char * in
, const char * pat
, const char * rep
, const size_t maxLen
, char * out
) {
678 size_t inLen
, patLen
, repLen
;
683 patLen
= strlen (pat
);
684 repLen
= strlen (rep
);
688 while ((si
< inLen
) && (di
< maxLen
)) {
689 if (patLen
<= (inLen
- si
)) {
690 if (strncasecmp (& in
[si
], pat
, patLen
) == 0) {
691 if (repLen
> (maxLen
- di
)) {
692 repLen
= maxLen
- di
;
695 strncpy (& out
[di
], rep
, repLen
);
707 return (! truncated
);
710 // Provide missing math routines for MSVC.
712 static double round (double val
) {
714 return (ceil (val
- 0.5));
715 return (floor (val
+ 0.5));
718 static double fmin (double a
, double b
) {
719 return ((a
< b
) ? a
: b
);
722 static double fmax (double a
, double b
) {
723 return ((a
> b
) ? a
: b
);
727 // Simple clamp routine.
728 static double Clamp (const double val
, const double lower
, const double upper
) {
729 return (fmin (fmax (val
, lower
), upper
));
732 // Performs linear interpolation.
733 static double Lerp (const double a
, const double b
, const double f
) {
734 return (a
+ (f
* (b
- a
)));
737 // Performs a high-passed triangular probability density function dither from
738 // a double to an integer. It assumes the input sample is already scaled.
739 static int HpTpdfDither (const double in
, int * hpHist
) {
740 const double PRNG_SCALE
= 1.0 / RAND_MAX
;
745 out
= round (in
+ (PRNG_SCALE
* (prn
- (* hpHist
))));
750 // Allocates an array of doubles.
751 static double * CreateArray (const size_t n
) {
754 a
= (double *) calloc (n
, sizeof (double));
756 fprintf (stderr
, "Error: Out of memory.\n");
762 // Frees an array of doubles.
763 static void DestroyArray (const double * a
) {
767 // Complex number routines. All outputs must be non-NULL.
769 // Magnitude/absolute value.
770 static double ComplexAbs (const double r
, const double i
) {
771 return (sqrt ((r
* r
) + (i
* i
)));
775 static void ComplexMul (const double aR
, const double aI
, const double bR
, const double bI
, double * outR
, double * outI
) {
776 (* outR
) = (aR
* bR
) - (aI
* bI
);
777 (* outI
) = (aI
* bR
) + (aR
* bI
);
781 static void ComplexExp (const double inR
, const double inI
, double * outR
, double * outI
) {
785 (* outR
) = e
* cos (inI
);
786 (* outI
) = e
* sin (inI
);
789 /* Fast Fourier transform routines. The number of points must be a power of
790 * two. In-place operation is possible only if both the real and imaginary
791 * parts are in-place together.
794 // Performs bit-reversal ordering.
795 static void FftArrange (const size_t n
, const double * inR
, const double * inI
, double * outR
, double * outI
) {
799 if ((inR
== outR
) && (inI
== outI
)) {
800 // Handle in-place arrangement.
802 for (k
= 0; k
< n
; k
++) {
812 while (rk
& (m
>>= 1))
817 // Handle copy arrangement.
819 for (k
= 0; k
< n
; k
++) {
823 while (rk
& (m
>>= 1))
830 // Performs the summation.
831 static void FftSummation (const size_t n
, const double s
, double * re
, double * im
) {
834 double vR
, vI
, wR
, wI
;
839 for (m
= 1, m2
= 2; m
< n
; m
<<= 1, m2
<<= 1) {
840 // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
841 vR
= sin (0.5 * pi
/ m
);
844 // w = Complex (1.0, 0.0)
847 for (i
= 0; i
< m
; i
++) {
848 for (k
= i
; k
< n
; k
+= m2
) {
850 // t = ComplexMul (w, out [km2])
851 tR
= (wR
* re
[mk
]) - (wI
* im
[mk
]);
852 tI
= (wR
* im
[mk
]) + (wI
* re
[mk
]);
853 // out [mk] = ComplexSub (out [k], t)
854 re
[mk
] = re
[k
] - tR
;
855 im
[mk
] = im
[k
] - tI
;
856 // out [k] = ComplexAdd (out [k], t)
860 // t = ComplexMul (v, w)
861 tR
= (vR
* wR
) - (vI
* wI
);
862 tI
= (vR
* wI
) + (vI
* wR
);
863 // w = ComplexAdd (w, t)
870 // Performs a forward FFT.
871 static void FftForward (const size_t n
, const double * inR
, const double * inI
, double * outR
, double * outI
) {
872 FftArrange (n
, inR
, inI
, outR
, outI
);
873 FftSummation (n
, 1.0, outR
, outI
);
876 // Performs an inverse FFT.
877 static void FftInverse (const size_t n
, const double * inR
, const double * inI
, double * outR
, double * outI
) {
881 FftArrange (n
, inR
, inI
, outR
, outI
);
882 FftSummation (n
, -1.0, outR
, outI
);
884 for (i
= 0; i
< n
; i
++) {
890 /* Calculate the complex helical sequence (or discrete-time analytical
891 * signal) of the given input using the Hilbert transform. Given the
892 * negative natural logarithm of a signal's magnitude response, the imaginary
893 * components can be used as the angles for minimum-phase reconstruction.
895 static void Hilbert (const size_t n
, const double * in
, double * outR
, double * outI
) {
899 // Handle in-place operation.
900 for (i
= 0; i
< n
; i
++)
903 // Handle copy operation.
904 for (i
= 0; i
< n
; i
++) {
909 FftForward (n
, outR
, outI
, outR
, outI
);
910 /* Currently the Fourier routines operate only on point counts that are
911 * powers of two. If that changes and n is odd, the following conditional
912 * should be: i < (n + 1) / 2.
914 for (i
= 1; i
< (n
/ 2); i
++) {
918 // If n is odd, the following increment should be skipped.
920 for (; i
< n
; i
++) {
924 FftInverse (n
, outR
, outI
, outR
, outI
);
927 /* Calculate the magnitude response of the given input. This is used in
928 * place of phase decomposition, since the phase residuals are discarded for
929 * minimum phase reconstruction. The mirrored half of the response is also
932 static void MagnitudeResponse (const size_t n
, const double * inR
, const double * inI
, double * out
) {
933 const size_t m
= 1 + (n
/ 2);
936 for (i
= 0; i
< m
; i
++)
937 out
[i
] = fmax (ComplexAbs (inR
[i
], inI
[i
]), EPSILON
);
940 /* Apply a range limit (in dB) to the given magnitude response. This is used
941 * to adjust the effects of the diffuse-field average on the equalization
944 static void LimitMagnitudeResponse (const size_t n
, const double limit
, const double * in
, double * out
) {
945 const size_t m
= 1 + (n
/ 2);
947 size_t i
, lower
, upper
;
950 halfLim
= limit
/ 2.0;
951 // Convert the response to dB.
952 for (i
= 0; i
< m
; i
++)
953 out
[i
] = 20.0 * log10 (in
[i
]);
954 // Use six octaves to calculate the average magnitude of the signal.
955 lower
= ((size_t) ceil (n
/ pow (2.0, 8.0))) - 1;
956 upper
= ((size_t) floor (n
/ pow (2.0, 2.0))) - 1;
958 for (i
= lower
; i
<= upper
; i
++)
960 ave
/= upper
- lower
+ 1;
961 // Keep the response within range of the average magnitude.
962 for (i
= 0; i
< m
; i
++)
963 out
[i
] = Clamp (out
[i
], ave
- halfLim
, ave
+ halfLim
);
964 // Convert the response back to linear magnitude.
965 for (i
= 0; i
< m
; i
++)
966 out
[i
] = pow (10.0, out
[i
] / 20.0);
969 /* Reconstructs the minimum-phase component for the given magnitude response
970 * of a signal. This is equivalent to phase recomposition, sans the missing
971 * residuals (which were discarded). The mirrored half of the response is
974 static void MinimumPhase (const size_t n
, const double * in
, double * outR
, double * outI
) {
975 const size_t m
= 1 + (n
/ 2);
976 double * mags
= NULL
;
980 mags
= CreateArray (n
);
981 for (i
= 0; i
< m
; i
++) {
982 mags
[i
] = fmax (in
[i
], EPSILON
);
983 outR
[i
] = -log (mags
[i
]);
985 for (; i
< n
; i
++) {
986 mags
[i
] = mags
[n
- i
];
987 outR
[i
] = outR
[n
- i
];
989 Hilbert (n
, outR
, outR
, outI
);
990 // Remove any DC offset the filter has.
993 for (i
= 1; i
< n
; i
++) {
994 ComplexExp (0.0, outI
[i
], & aR
, & aI
);
995 ComplexMul (mags
[i
], 0.0, aR
, aI
, & outR
[i
], & outI
[i
]);
1000 // Read a binary value of the specified byte order and byte size from a file,
1001 // storing it as a 32-bit unsigned integer.
1002 static int ReadBin4 (FILE * fp
, const char * filename
, const ByteOrderT order
, const uint bytes
, uint4
* out
) {
1007 if (fread (in
, 1, bytes
, fp
) != bytes
) {
1008 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1014 for (i
= 0; i
< bytes
; i
++)
1015 accum
= (accum
<< 8) | in
[bytes
- i
- 1];
1018 for (i
= 0; i
< bytes
; i
++)
1019 accum
= (accum
<< 8) | in
[i
];
1028 // Read a binary value of the specified byte order from a file, storing it as
1029 // a 64-bit unsigned integer.
1030 static int ReadBin8 (FILE * fp
, const char * filename
, const ByteOrderT order
, uint8
* out
) {
1035 if (fread (in
, 1, 8, fp
) != 8) {
1036 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1042 for (i
= 0; i
< 8; i
++)
1043 accum
= (accum
<< 8) | in
[8 - i
- 1];
1046 for (i
= 0; i
< 8; i
++)
1047 accum
= (accum
<< 8) | in
[i
];
1056 // Write an ASCII string to a file.
1057 static int WriteAscii (const char * out
, FILE * fp
, const char * filename
) {
1061 if (fwrite (out
, 1, len
, fp
) != len
) {
1063 fprintf (stderr
, "Error: Bad write to file '%s'.\n", filename
);
1069 // Write a binary value of the given byte order and byte size to a file,
1070 // loading it from a 32-bit unsigned integer.
1071 static int WriteBin4 (const ByteOrderT order
, const uint bytes
, const uint4 in
, FILE * fp
, const char * filename
) {
1077 for (i
= 0; i
< bytes
; i
++)
1078 out
[i
] = (in
>> (i
* 8)) & 0x000000FF;
1081 for (i
= 0; i
< bytes
; i
++)
1082 out
[bytes
- i
- 1] = (in
>> (i
* 8)) & 0x000000FF;
1087 if (fwrite (out
, 1, bytes
, fp
) != bytes
) {
1088 fprintf (stderr
, "Error: Bad write to file '%s'.\n", filename
);
1094 /* Read a binary value of the specified type, byte order, and byte size from
1095 * a file, converting it to a double. For integer types, the significant
1096 * bits are used to normalize the result. The sign of bits determines
1097 * whether they are padded toward the MSB (negative) or LSB (positive).
1098 * Floating-point types are not normalized.
1100 static int ReadBinAsDouble (FILE * fp
, const char * filename
, const ByteOrderT order
, const ElementTypeT type
, const uint bytes
, const int bits
, double * out
) {
1113 if (! ReadBin8 (fp
, filename
, order
, & v8
. ui
))
1118 if (! ReadBin4 (fp
, filename
, order
, bytes
, & v4
. ui
))
1120 if (type
== ET_FP
) {
1121 (* out
) = (double) v4
. f
;
1124 v4
. ui
>>= (8 * bytes
) - bits
;
1126 v4
. ui
&= (0xFFFFFFFF >> (32 + bits
));
1127 if (v4
. ui
& (1 << (abs (bits
) - 1)))
1128 v4
. ui
|= (0xFFFFFFFF << abs (bits
));
1129 (* out
) = v4
. i
/ ((double) (1 << (abs (bits
) - 1)));
1135 /* Read an ascii value of the specified type from a file, converting it to a
1136 * double. For integer types, the significant bits are used to normalize the
1137 * result. The sign of the bits should always be positive. This also skips
1138 * up to one separator character before the element itself.
1140 static int ReadAsciiAsDouble (TokenReaderT
* tr
, const char * filename
, const ElementTypeT type
, const uint bits
, double * out
) {
1143 if (TrIsOperator (tr
, ","))
1144 TrReadOperator (tr
, ",");
1145 else if (TrIsOperator (tr
, ":"))
1146 TrReadOperator (tr
, ":");
1147 else if (TrIsOperator (tr
, ";"))
1148 TrReadOperator (tr
, ";");
1149 else if (TrIsOperator (tr
, "|"))
1150 TrReadOperator (tr
, "|");
1151 if (type
== ET_FP
) {
1152 if (! TrReadFloat (tr
, -1.0 / 0.0, 1.0 / 0.0, out
)) {
1153 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1157 if (! TrReadInt (tr
, -(1 << (bits
- 1)), (1 << (bits
- 1)) - 1, & v
)) {
1158 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1161 (* out
) = v
/ ((double) ((1 << (bits
- 1)) - 1));
1166 // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
1167 // the source parameters and data set metrics.
1168 static int ReadWaveFormat (FILE * fp
, const ByteOrderT order
, const uint hrirRate
, SourceRefT
* src
) {
1169 uint4 fourCC
, chunkSize
;
1170 uint4 format
, channels
, rate
, dummy
, block
, size
, bits
;
1175 fseek (fp
, chunkSize
, SEEK_CUR
);
1176 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1177 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & chunkSize
)))
1179 } while (fourCC
!= FOURCC_FMT
);
1180 if ((! ReadBin4 (fp
, src
-> mPath
, order
, 2, & format
)) ||
1181 (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & channels
)) ||
1182 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & rate
)) ||
1183 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & dummy
)) ||
1184 (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & block
)))
1187 if (chunkSize
> 14) {
1188 if (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & size
))
1196 if (format
== WAVE_FORMAT_EXTENSIBLE
) {
1197 fseek (fp
, 2, SEEK_CUR
);
1198 if (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & bits
))
1202 fseek (fp
, 4, SEEK_CUR
);
1203 if (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & format
))
1205 fseek (fp
, chunkSize
- 26, SEEK_CUR
);
1209 fseek (fp
, chunkSize
- 16, SEEK_CUR
);
1211 fseek (fp
, chunkSize
- 14, SEEK_CUR
);
1213 if ((format
!= WAVE_FORMAT_PCM
) && (format
!= WAVE_FORMAT_IEEE_FLOAT
)) {
1214 fprintf (stderr
, "Error: Unsupported WAVE format in file '%s'.\n", src
-> mPath
);
1217 if (src
-> mChannel
>= channels
) {
1218 fprintf (stderr
, "Error: Missing source channel in WAVE file '%s'.\n", src
-> mPath
);
1221 if (rate
!= hrirRate
) {
1222 fprintf (stderr
, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src
-> mPath
);
1225 if (format
== WAVE_FORMAT_PCM
) {
1226 if ((size
< 2) || (size
> 4)) {
1227 fprintf (stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
-> mPath
);
1230 if ((bits
< 16) || (bits
> (8 * size
))) {
1231 fprintf (stderr
, "Error: Bad significant bits in WAVE file '%s'.\n", src
-> mPath
);
1234 src
-> mType
= ET_INT
;
1236 if ((size
!= 4) && (size
!= 8)) {
1237 fprintf (stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
-> mPath
);
1240 src
-> mType
= ET_FP
;
1242 src
-> mSize
= size
;
1243 src
-> mBits
= (int) bits
;
1244 src
-> mSkip
= channels
;
1248 // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
1249 static int ReadWaveData (FILE * fp
, const SourceRefT
* src
, const ByteOrderT order
, const size_t n
, double * hrir
) {
1250 int pre
, post
, skip
;
1253 pre
= (int) (src
-> mSize
* src
-> mChannel
);
1254 post
= (int) (src
-> mSize
* (src
-> mSkip
- src
-> mChannel
- 1));
1256 for (i
= 0; i
< n
; i
++) {
1259 fseek (fp
, skip
, SEEK_CUR
);
1260 if (! ReadBinAsDouble (fp
, src
-> mPath
, order
, src
-> mType
, src
-> mSize
, src
-> mBits
, & hrir
[i
]))
1265 fseek (fp
, skip
, SEEK_CUR
);
1269 // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
1271 static int ReadWaveList (FILE * fp
, const SourceRefT
* src
, const ByteOrderT order
, const size_t n
, double * hrir
) {
1272 uint4 fourCC
, chunkSize
, listSize
;
1273 size_t block
, skip
, count
, offset
, i
;
1277 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1278 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & chunkSize
)))
1280 if (fourCC
== FOURCC_DATA
) {
1281 block
= src
-> mSize
* src
-> mSkip
;
1282 count
= chunkSize
/ block
;
1283 if (count
< (src
-> mOffset
+ n
)) {
1284 fprintf (stderr
, "Error: Bad read from file '%s'.\n", src
-> mPath
);
1287 fseek (fp
, src
-> mOffset
* block
, SEEK_CUR
);
1288 if (! ReadWaveData (fp
, src
, order
, n
, & hrir
[0]))
1291 } else if (fourCC
== FOURCC_LIST
) {
1292 if (! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
))
1295 if (fourCC
== FOURCC_WAVL
)
1299 fseek (fp
, chunkSize
, SEEK_CUR
);
1301 listSize
= chunkSize
;
1302 block
= src
-> mSize
* src
-> mSkip
;
1303 skip
= src
-> mOffset
;
1306 while ((offset
< n
) && (listSize
> 8)) {
1307 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1308 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & chunkSize
)))
1310 listSize
-= 8 + chunkSize
;
1311 if (fourCC
== FOURCC_DATA
) {
1312 count
= chunkSize
/ block
;
1314 fseek (fp
, skip
* block
, SEEK_CUR
);
1315 chunkSize
-= skip
* block
;
1318 if (count
> (n
- offset
))
1320 if (! ReadWaveData (fp
, src
, order
, count
, & hrir
[offset
]))
1322 chunkSize
-= count
* block
;
1324 lastSample
= hrir
[offset
- 1];
1329 } else if (fourCC
== FOURCC_SLNT
) {
1330 if (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & count
))
1336 if (count
> (n
- offset
))
1338 for (i
= 0; i
< count
; i
++)
1339 hrir
[offset
+ i
] = lastSample
;
1347 fseek (fp
, chunkSize
, SEEK_CUR
);
1350 fprintf (stderr
, "Error: Bad read from file '%s'.\n", src
-> mPath
);
1356 // Load a source HRIR from a RIFF/RIFX WAVE file.
1357 static int LoadWaveSource (FILE * fp
, SourceRefT
* src
, const uint hrirRate
, const size_t n
, double * hrir
) {
1358 uint4 fourCC
, dummy
;
1361 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1362 (! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & dummy
)))
1364 if (fourCC
== FOURCC_RIFF
) {
1366 } else if (fourCC
== FOURCC_RIFX
) {
1369 fprintf (stderr
, "Error: No RIFF/RIFX chunk in file '%s'.\n", src
-> mPath
);
1372 if (! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
))
1374 if (fourCC
!= FOURCC_WAVE
) {
1375 fprintf (stderr
, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src
-> mPath
);
1378 if (! ReadWaveFormat (fp
, order
, hrirRate
, src
))
1380 if (! ReadWaveList (fp
, src
, order
, n
, hrir
))
1385 // Load a source HRIR from a binary file.
1386 static int LoadBinarySource (FILE * fp
, const SourceRefT
* src
, const ByteOrderT order
, const size_t n
, double * hrir
) {
1389 fseek (fp
, src
-> mOffset
, SEEK_SET
);
1390 for (i
= 0; i
< n
; i
++) {
1391 if (! ReadBinAsDouble (fp
, src
-> mPath
, order
, src
-> mType
, src
-> mSize
, src
-> mBits
, & hrir
[i
]))
1393 if (src
-> mSkip
> 0)
1394 fseek (fp
, src
-> mSkip
, SEEK_CUR
);
1399 // Load a source HRIR from an ASCII text file containing a list of elements
1400 // separated by whitespace or common list operators (',', ';', ':', '|').
1401 static int LoadAsciiSource (FILE * fp
, const SourceRefT
* src
, const size_t n
, double * hrir
) {
1406 TrSetup (fp
, NULL
, & tr
);
1407 for (i
= 0; i
< src
-> mOffset
; i
++) {
1408 if (! ReadAsciiAsDouble (& tr
, src
-> mPath
, src
-> mType
, src
-> mBits
, & dummy
))
1411 for (i
= 0; i
< n
; i
++) {
1412 if (! ReadAsciiAsDouble (& tr
, src
-> mPath
, src
-> mType
, src
-> mBits
, & hrir
[i
]))
1414 for (j
= 0; j
< src
-> mSkip
; j
++) {
1415 if (! ReadAsciiAsDouble (& tr
, src
-> mPath
, src
-> mType
, src
-> mBits
, & dummy
))
1422 // Load a source HRIR from a supported file type.
1423 static int LoadSource (SourceRefT
* src
, const uint hrirRate
, const size_t n
, double * hrir
) {
1427 if (src
-> mFormat
== SF_ASCII
)
1428 fp
= fopen (src
-> mPath
, "r");
1430 fp
= fopen (src
-> mPath
, "rb");
1432 fprintf (stderr
, "Error: Could not open source file '%s'.\n", src
-> mPath
);
1435 if (src
-> mFormat
== SF_WAVE
)
1436 result
= LoadWaveSource (fp
, src
, hrirRate
, n
, hrir
);
1437 else if (src
-> mFormat
== SF_BIN_LE
)
1438 result
= LoadBinarySource (fp
, src
, BO_LITTLE
, n
, hrir
);
1439 else if (src
-> mFormat
== SF_BIN_BE
)
1440 result
= LoadBinarySource (fp
, src
, BO_BIG
, n
, hrir
);
1442 result
= LoadAsciiSource (fp
, src
, n
, hrir
);
1447 // Calculate the magnitude response of an HRIR and average it with any
1448 // existing responses for its elevation and azimuth.
1449 static void AverageHrirMagnitude (const double * hrir
, const double f
, const uint ei
, const uint ai
, const HrirDataT
* hData
) {
1450 double * re
= NULL
, * im
= NULL
;
1453 n
= hData
-> mFftSize
;
1454 re
= CreateArray (n
);
1455 im
= CreateArray (n
);
1456 for (i
= 0; i
< hData
-> mIrPoints
; i
++) {
1460 for (; i
< n
; i
++) {
1464 FftForward (n
, re
, im
, re
, im
);
1465 MagnitudeResponse (n
, re
, im
, re
);
1467 j
= (hData
-> mEvOffset
[ei
] + ai
) * hData
-> mIrSize
;
1468 for (i
= 0; i
< m
; i
++)
1469 hData
-> mHrirs
[j
+ i
] = Lerp (hData
-> mHrirs
[j
+ i
], re
[i
], f
);
1474 /* Calculate the contribution of each HRIR to the diffuse-field average based
1475 * on the area of its surface patch. All patches are centered at the HRIR
1476 * coordinates on the unit sphere and are measured by solid angle.
1478 static void CalculateDfWeights (const HrirDataT
* hData
, double * weights
) {
1480 double evs
, sum
, ev
, up_ev
, down_ev
, solidAngle
;
1482 evs
= 90.0 / (hData
-> mEvCount
- 1);
1484 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++) {
1485 // For each elevation, calculate the upper and lower limits of the
1487 ev
= -90.0 + (ei
* 2.0 * evs
);
1488 if (ei
< (hData
-> mEvCount
- 1))
1489 up_ev
= (ev
+ evs
) * M_PI
/ 180.0;
1493 down_ev
= (ev
- evs
) * M_PI
/ 180.0;
1495 down_ev
= -M_PI
/ 2.0;
1496 // Calculate the area of the patch band.
1497 solidAngle
= 2.0 * M_PI
* (sin (up_ev
) - sin (down_ev
));
1498 // Each weight is the area of one patch.
1499 weights
[ei
] = solidAngle
/ hData
-> mAzCount
[ei
];
1500 // Sum the total surface area covered by the HRIRs.
1503 // Normalize the weights given the total surface coverage.
1504 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++)
1505 weights
[ei
] /= sum
;
1508 /* Calculate the diffuse-field average from the given magnitude responses of
1509 * the HRIR set. Weighting can be applied to compensate for the varying
1510 * surface area covered by each HRIR. The final average can then be limited
1511 * by the specified magnitude range (in positive dB; 0.0 to skip).
1513 static void CalculateDiffuseFieldAverage (const HrirDataT
* hData
, const int weighted
, const double limit
, double * dfa
) {
1514 double * weights
= NULL
;
1516 size_t count
, step
, start
, end
, m
, j
, i
;
1519 weights
= CreateArray (hData
-> mEvCount
);
1521 // Use coverage weighting to calculate the average.
1522 CalculateDfWeights (hData
, weights
);
1524 // If coverage weighting is not used, the weights still need to be
1525 // averaged by the number of HRIRs.
1527 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++)
1528 count
+= hData
-> mAzCount
[ei
];
1529 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++)
1530 weights
[ei
] = 1.0 / count
;
1532 ei
= hData
-> mEvStart
;
1534 step
= hData
-> mIrSize
;
1535 start
= hData
-> mEvOffset
[ei
] * step
;
1536 end
= hData
-> mIrCount
* step
;
1537 m
= 1 + (hData
-> mFftSize
/ 2);
1538 for (i
= 0; i
< m
; i
++)
1540 for (j
= start
; j
< end
; j
+= step
) {
1541 // Get the weight for this HRIR's contribution.
1542 weight
= weights
[ei
];
1543 // Add this HRIR's weighted power average to the total.
1544 for (i
= 0; i
< m
; i
++)
1545 dfa
[i
] += weight
* hData
-> mHrirs
[j
+ i
] * hData
-> mHrirs
[j
+ i
];
1546 // Determine the next weight to use.
1548 if (ai
>= hData
-> mAzCount
[ei
]) {
1553 // Finish the average calculation and keep it from being too small.
1554 for (i
= 0; i
< m
; i
++)
1555 dfa
[i
] = fmax (sqrt (dfa
[i
]), EPSILON
);
1556 // Apply a limit to the magnitude range of the diffuse-field average if
1559 LimitMagnitudeResponse (hData
-> mFftSize
, limit
, dfa
, dfa
);
1560 DestroyArray (weights
);
1563 // Perform diffuse-field equalization on the magnitude responses of the HRIR
1564 // set using the given average response.
1565 static void DiffuseFieldEqualize (const double * dfa
, const HrirDataT
* hData
) {
1566 size_t step
, start
, end
, m
, j
, i
;
1568 step
= hData
-> mIrSize
;
1569 start
= hData
-> mEvOffset
[hData
-> mEvStart
] * step
;
1570 end
= hData
-> mIrCount
* step
;
1571 m
= 1 + (hData
-> mFftSize
/ 2);
1572 for (j
= start
; j
< end
; j
+= step
) {
1573 for (i
= 0; i
< m
; i
++)
1574 hData
-> mHrirs
[j
+ i
] /= dfa
[i
];
1578 // Perform minimum-phase reconstruction using the magnitude responses of the
1580 static void ReconstructHrirs (const HrirDataT
* hData
) {
1581 double * re
= NULL
, * im
= NULL
;
1582 size_t step
, start
, end
, n
, j
, i
;
1584 step
= hData
-> mIrSize
;
1585 start
= hData
-> mEvOffset
[hData
-> mEvStart
] * step
;
1586 end
= hData
-> mIrCount
* step
;
1587 n
= hData
-> mFftSize
;
1588 re
= CreateArray (n
);
1589 im
= CreateArray (n
);
1590 for (j
= start
; j
< end
; j
+= step
) {
1591 MinimumPhase (n
, & hData
-> mHrirs
[j
], re
, im
);
1592 FftInverse (n
, re
, im
, re
, im
);
1593 for (i
= 0; i
< hData
-> mIrPoints
; i
++)
1594 hData
-> mHrirs
[j
+ i
] = re
[i
];
1600 /* Given an elevation index and an azimuth, calculate the indices of the two
1601 * HRIRs that bound the coordinate along with a factor for calculating the
1602 * continous HRIR using interpolation.
1604 static void CalcAzIndices (const HrirDataT
* hData
, const uint ei
, const double az
, uint
* j0
, uint
* j1
, double * jf
) {
1608 af
= ((2.0 * M_PI
) + az
) * hData
-> mAzCount
[ei
] / (2.0 * M_PI
);
1609 ai
= ((uint
) af
) % hData
-> mAzCount
[ei
];
1611 (* j0
) = hData
-> mEvOffset
[ei
] + ai
;
1612 (* j1
) = hData
-> mEvOffset
[ei
] + ((ai
+ 1) % hData
-> mAzCount
[ei
]);
1616 /* Attempt to synthesize any missing HRIRs at the bottom elevations. Right
1617 * now this just blends the lowest elevation HRIRs together and applies some
1618 * attenuation and high frequency damping. It is a simple, if inaccurate
1621 static void SynthesizeHrirs (HrirDataT
* hData
) {
1623 size_t step
, n
, i
, j
;
1627 double lp
[4], s0
, s1
;
1629 if (hData
-> mEvStart
<= 0)
1631 step
= hData
-> mIrSize
;
1632 oi
= hData
-> mEvStart
;
1633 n
= hData
-> mIrPoints
;
1634 for (i
= 0; i
< n
; i
++)
1635 hData
-> mHrirs
[i
] = 0.0;
1636 for (a
= 0; a
< hData
-> mAzCount
[oi
]; a
++) {
1637 j
= (hData
-> mEvOffset
[oi
] + a
) * step
;
1638 for (i
= 0; i
< n
; i
++)
1639 hData
-> mHrirs
[i
] += hData
-> mHrirs
[j
+ i
] / hData
-> mAzCount
[oi
];
1641 for (e
= 1; e
< hData
-> mEvStart
; e
++) {
1642 of
= ((double) e
) / hData
-> mEvStart
;
1643 for (a
= 0; a
< hData
-> mAzCount
[e
]; a
++) {
1644 j
= (hData
-> mEvOffset
[e
] + a
) * step
;
1645 CalcAzIndices (hData
, oi
, a
* 2.0 * M_PI
/ hData
-> mAzCount
[e
], & j0
, & j1
, & jf
);
1652 for (i
= 0; i
< n
; i
++) {
1653 s0
= hData
-> mHrirs
[i
];
1654 s1
= Lerp (hData
-> mHrirs
[j0
+ i
], hData
-> mHrirs
[j1
+ i
], jf
);
1655 s0
= Lerp (s0
, s1
, of
);
1656 lp
[0] = Lerp (s0
, lp
[0], 0.15 - (0.15 * of
));
1657 lp
[1] = Lerp (lp
[0], lp
[1], 0.15 - (0.15 * of
));
1658 lp
[2] = Lerp (lp
[1], lp
[2], 0.15 - (0.15 * of
));
1659 lp
[3] = Lerp (lp
[2], lp
[3], 0.15 - (0.15 * of
));
1660 hData
-> mHrirs
[j
+ i
] = lp
[3];
1668 for (i
= 0; i
< n
; i
++) {
1669 s0
= hData
-> mHrirs
[i
];
1670 lp
[0] = Lerp (s0
, lp
[0], 0.15);
1671 lp
[1] = Lerp (lp
[0], lp
[1], 0.15);
1672 lp
[2] = Lerp (lp
[1], lp
[2], 0.15);
1673 lp
[3] = Lerp (lp
[2], lp
[3], 0.15);
1674 hData
-> mHrirs
[i
] = lp
[3];
1676 hData
-> mEvStart
= 0;
1679 // The following routines assume a full set of HRIRs for all elevations.
1681 // Normalize the HRIR set and slightly attenuate the result.
1682 static void NormalizeHrirs (const HrirDataT
* hData
) {
1683 size_t step
, end
, n
, j
, i
;
1686 step
= hData
-> mIrSize
;
1687 end
= hData
-> mIrCount
* step
;
1688 n
= hData
-> mIrPoints
;
1690 for (j
= 0; j
< end
; j
+= step
) {
1691 for (i
= 0; i
< n
; i
++)
1692 maxLevel
= fmax (fabs (hData
-> mHrirs
[j
+ i
]), maxLevel
);
1694 maxLevel
= 1.01 * maxLevel
;
1695 for (j
= 0; j
< end
; j
+= step
) {
1696 for (i
= 0; i
< n
; i
++)
1697 hData
-> mHrirs
[j
+ i
] /= maxLevel
;
1701 // Calculate the left-ear time delay using a spherical head model.
1702 static double CalcLTD (const double ev
, const double az
, const double rad
, const double dist
) {
1703 double azp
, dlp
, l
, al
;
1705 azp
= asin (cos (ev
) * sin (az
));
1706 dlp
= sqrt ((dist
* dist
) + (rad
* rad
) + (2.0 * dist
* rad
* sin (azp
)));
1707 l
= sqrt ((dist
* dist
) - (rad
* rad
));
1708 al
= (0.5 * M_PI
) + azp
;
1710 dlp
= l
+ (rad
* (al
- acos (rad
/ dist
)));
1711 return (dlp
/ 343.3);
1714 // Calculate the effective head-related time delays for the each HRIR, now
1715 // that they are minimum-phase.
1716 static void CalculateHrtds (HrirDataT
* hData
) {
1717 double minHrtd
, maxHrtd
;
1723 for (e
= 0; e
< hData
-> mEvCount
; e
++) {
1724 for (a
= 0; a
< hData
-> mAzCount
[e
]; a
++) {
1725 j
= hData
-> mEvOffset
[e
] + a
;
1726 t
= CalcLTD ((-90.0 + (e
* 180.0 / (hData
-> mEvCount
- 1))) * M_PI
/ 180.0,
1727 (a
* 360.0 / hData
-> mAzCount
[e
]) * M_PI
/ 180.0,
1728 hData
-> mRadius
, hData
-> mDistance
);
1729 hData
-> mHrtds
[j
] = t
;
1730 maxHrtd
= fmax (t
, maxHrtd
);
1731 minHrtd
= fmin (t
, minHrtd
);
1735 for (j
= 0; j
< hData
-> mIrCount
; j
++)
1736 hData
-> mHrtds
[j
] -= minHrtd
;
1737 hData
-> mMaxHrtd
= maxHrtd
;
1740 // Store the OpenAL Soft HRTF data set.
1741 static int StoreMhr (const HrirDataT
* hData
, const char * filename
) {
1744 size_t step
, end
, n
, j
, i
;
1747 if ((fp
= fopen (filename
, "wb")) == NULL
) {
1748 fprintf (stderr
, "Error: Could not open MHR file '%s'.\n", filename
);
1751 if (! WriteAscii (MHR_FORMAT
, fp
, filename
))
1753 if (! WriteBin4 (BO_LITTLE
, 4, (uint4
) hData
-> mIrRate
, fp
, filename
))
1755 if (! WriteBin4 (BO_LITTLE
, 1, (uint4
) hData
-> mIrPoints
, fp
, filename
))
1757 if (! WriteBin4 (BO_LITTLE
, 1, (uint4
) hData
-> mEvCount
, fp
, filename
))
1759 for (e
= 0; e
< hData
-> mEvCount
; e
++) {
1760 if (! WriteBin4 (BO_LITTLE
, 1, (uint4
) hData
-> mAzCount
[e
], fp
, filename
))
1763 step
= hData
-> mIrSize
;
1764 end
= hData
-> mIrCount
* step
;
1765 n
= hData
-> mIrPoints
;
1767 for (j
= 0; j
< end
; j
+= step
) {
1769 for (i
= 0; i
< n
; i
++) {
1770 v
= HpTpdfDither (32767.0 * hData
-> mHrirs
[j
+ i
], & hpHist
);
1771 if (! WriteBin4 (BO_LITTLE
, 2, (uint4
) v
, fp
, filename
))
1775 for (j
= 0; j
< hData
-> mIrCount
; j
++) {
1776 // v = (int) fmin (round (44100.0 * hData -> mHrtds [j]), MAX_HRTD);
1777 v
= fmin (hData
-> mHrtds
[j
], MAX_HRTD
);
1778 if (! WriteBin4 (BO_LITTLE
, 1, (uint4
) v
, fp
, filename
))
1785 // Store the OpenAL Soft built-in table.
1786 static int StoreTable (const HrirDataT
* hData
, const char * filename
) {
1788 size_t step
, end
, n
, j
, i
;
1790 char text
[128 + 1];
1792 if ((fp
= fopen (filename
, "wb")) == NULL
) {
1793 fprintf (stderr
, "Error: Could not open table file '%s'.\n", filename
);
1796 snprintf (text
, 128, "/* Elevation metrics */\n"
1797 "static const ALubyte defaultAzCount[%u] = { ", hData
-> mEvCount
);
1798 if (! WriteAscii (text
, fp
, filename
))
1800 for (i
= 0; i
< hData
-> mEvCount
; i
++) {
1801 snprintf (text
, 128, "%u, ", hData
-> mAzCount
[i
]);
1802 if (! WriteAscii (text
, fp
, filename
))
1805 snprintf (text
, 128, "};\n"
1806 "static const ALushort defaultEvOffset[%u] = { ", hData
-> mEvCount
);
1807 if (! WriteAscii (text
, fp
, filename
))
1809 for (i
= 0; i
< hData
-> mEvCount
; i
++) {
1810 snprintf (text
, 128, "%u, ", hData
-> mEvOffset
[i
]);
1811 if (! WriteAscii (text
, fp
, filename
))
1814 step
= hData
-> mIrSize
;
1815 end
= hData
-> mIrCount
* step
;
1816 n
= hData
-> mIrPoints
;
1817 snprintf (text
, 128, "};\n\n"
1818 "/* HRIR Coefficients */\n"
1819 "static const ALshort defaultCoeffs[%u] =\n{\n", end
);
1820 if (! WriteAscii (text
, fp
, filename
))
1823 for (j
= 0; j
< end
; j
+= step
) {
1824 if (! WriteAscii (" ", fp
, filename
))
1827 for (i
= 0; i
< n
; i
++) {
1828 v
= HpTpdfDither (32767.0 * hData
-> mHrirs
[j
+ i
], & hpHist
);
1829 snprintf (text
, 128, "%+d, ", v
);
1830 if (! WriteAscii (text
, fp
, filename
))
1833 if (! WriteAscii ("\n", fp
, filename
))
1836 snprintf (text
, 128, "};\n\n"
1837 "/* HRIR Delays */\n"
1838 "static const ALubyte defaultDelays[%d] =\n{\n"
1839 " ", hData
-> mIrCount
);
1840 if (! WriteAscii (text
, fp
, filename
))
1842 for (j
= 0; j
< hData
-> mIrCount
; j
++) {
1843 v
= (int) fmin (round (44100.0 * hData
-> mHrtds
[j
]), MAX_HRTD
);
1844 snprintf (text
, 128, "%d, ", v
);
1845 if (! WriteAscii (text
, fp
, filename
))
1848 if (! WriteAscii ("\n};\n\n"
1849 "/* Default HRTF Definition */\n", fp
, filename
))
1851 snprintf (text
, 128, "static const struct Hrtf DefaultHrtf = {\n"
1852 " %u, %u, %u, defaultAzCount, defaultEvOffset,\n",
1853 hData
-> mIrRate
, hData
-> mIrPoints
, hData
-> mEvCount
);
1854 if (! WriteAscii (text
, fp
, filename
))
1856 if (! WriteAscii (" defaultCoeffs, defaultDelays, NULL\n"
1857 "};\n", fp
, filename
))
1863 // Process the data set definition to read and validate the data set metrics.
1864 static int ProcessMetrics (TokenReaderT
* tr
, const size_t fftSize
, const size_t truncSize
, HrirDataT
* hData
) {
1865 char ident
[MAX_IDENT_LEN
+ 1];
1870 int hasRate
= 0, hasPoints
= 0, hasAzimuths
= 0;
1871 int hasRadius
= 0, hasDistance
= 0;
1873 while (! (hasRate
&& hasPoints
&& hasAzimuths
&& hasRadius
&& hasDistance
)) {
1874 TrIndication (tr
, & line
, & col
);
1875 if (! TrReadIdent (tr
, MAX_IDENT_LEN
, ident
))
1877 if (strcasecmp (ident
, "rate") == 0) {
1879 TrErrorAt (tr
, line
, col
, "Redefinition of 'rate'.\n");
1882 if (! TrReadOperator (tr
, "="))
1884 if (! TrReadInt (tr
, MIN_RATE
, MAX_RATE
, & intVal
))
1886 hData
-> mIrRate
= intVal
;
1888 } else if (strcasecmp (ident
, "points") == 0) {
1890 TrErrorAt (tr
, line
, col
, "Redefinition of 'points'.\n");
1893 if (! TrReadOperator (tr
, "="))
1895 TrIndication (tr
, & line
, & col
);
1896 if (! TrReadInt (tr
, MIN_POINTS
, MAX_POINTS
, & intVal
))
1898 points
= (size_t) intVal
;
1899 if ((fftSize
> 0) && (points
> fftSize
)) {
1900 TrErrorAt (tr
, line
, col
, "Value exceeds the overriden FFT size.\n");
1903 if (points
< truncSize
) {
1904 TrErrorAt (tr
, line
, col
, "Value is below the truncation size.\n");
1907 hData
-> mIrPoints
= points
;
1908 hData
-> mFftSize
= fftSize
;
1911 while (points
< (4 * hData
-> mIrPoints
))
1913 hData
-> mFftSize
= points
;
1914 hData
-> mIrSize
= 1 + (points
/ 2);
1916 hData
-> mFftSize
= fftSize
;
1917 hData
-> mIrSize
= 1 + (fftSize
/ 2);
1918 if (points
> hData
-> mIrSize
)
1919 hData
-> mIrSize
= points
;
1922 } else if (strcasecmp (ident
, "azimuths") == 0) {
1924 TrErrorAt (tr
, line
, col
, "Redefinition of 'azimuths'.\n");
1927 if (! TrReadOperator (tr
, "="))
1929 hData
-> mIrCount
= 0;
1930 hData
-> mEvCount
= 0;
1931 hData
-> mEvOffset
[0] = 0;
1933 if (! TrReadInt (tr
, MIN_AZ_COUNT
, MAX_AZ_COUNT
, & intVal
))
1935 hData
-> mAzCount
[hData
-> mEvCount
] = intVal
;
1936 hData
-> mIrCount
+= intVal
;
1937 hData
-> mEvCount
++;
1938 if (! TrIsOperator (tr
, ","))
1940 if (hData
-> mEvCount
>= MAX_EV_COUNT
) {
1941 TrError (tr
, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT
);
1944 hData
-> mEvOffset
[hData
-> mEvCount
] = hData
-> mEvOffset
[hData
-> mEvCount
- 1] + intVal
;
1945 TrReadOperator (tr
, ",");
1947 if (hData
-> mEvCount
< MIN_EV_COUNT
) {
1948 TrErrorAt (tr
, line
, col
, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT
);
1952 } else if (strcasecmp (ident
, "radius") == 0) {
1954 TrErrorAt (tr
, line
, col
, "Redefinition of 'radius'.\n");
1957 if (! TrReadOperator (tr
, "="))
1959 if (! TrReadFloat (tr
, MIN_RADIUS
, MAX_RADIUS
, & fpVal
))
1961 hData
-> mRadius
= fpVal
;
1963 } else if (strcasecmp (ident
, "distance") == 0) {
1965 TrErrorAt (tr
, line
, col
, "Redefinition of 'distance'.\n");
1968 if (! TrReadOperator (tr
, "="))
1970 if (! TrReadFloat (tr
, MIN_DISTANCE
, MAX_DISTANCE
, & fpVal
))
1972 hData
-> mDistance
= fpVal
;
1975 TrErrorAt (tr
, line
, col
, "Expected a metric name.\n");
1978 TrSkipWhitespace (tr
);
1983 // Parse an index pair from the data set definition.
1984 static int ReadIndexPair (TokenReaderT
* tr
, const HrirDataT
* hData
, uint
* ei
, uint
* ai
) {
1987 if (! TrReadInt (tr
, 0, hData
-> mEvCount
, & intVal
))
1989 (* ei
) = (uint
) intVal
;
1990 if (! TrReadOperator (tr
, ","))
1992 if (! TrReadInt (tr
, 0, hData
-> mAzCount
[(* ei
)], & intVal
))
1994 (* ai
) = (uint
) intVal
;
1998 // Match the source format from a given identifier.
1999 static SourceFormatT
MatchSourceFormat (const char * ident
) {
2000 if (strcasecmp (ident
, "wave") == 0)
2002 else if (strcasecmp (ident
, "bin_le") == 0)
2004 else if (strcasecmp (ident
, "bin_be") == 0)
2006 else if (strcasecmp (ident
, "ascii") == 0)
2011 // Match the source element type from a given identifier.
2012 static ElementTypeT
MatchElementType (const char * ident
) {
2013 if (strcasecmp (ident
, "int") == 0)
2015 else if (strcasecmp (ident
, "fp") == 0)
2020 // Parse and validate a source reference from the data set definition.
2021 static int ReadSourceRef (TokenReaderT
* tr
, SourceRefT
* src
) {
2023 char ident
[MAX_IDENT_LEN
+ 1];
2026 TrIndication (tr
, & line
, & col
);
2027 if (! TrReadIdent (tr
, MAX_IDENT_LEN
, ident
))
2029 src
-> mFormat
= MatchSourceFormat (ident
);
2030 if (src
-> mFormat
== SF_NONE
) {
2031 TrErrorAt (tr
, line
, col
, "Expected a source format.\n");
2034 if (! TrReadOperator (tr
, "("))
2036 if (src
-> mFormat
== SF_WAVE
) {
2037 if (! TrReadInt (tr
, 0, MAX_WAVE_CHANNELS
, & intVal
))
2039 src
-> mType
= ET_NONE
;
2042 src
-> mChannel
= (uint
) intVal
;
2045 TrIndication (tr
, & line
, & col
);
2046 if (! TrReadIdent (tr
, MAX_IDENT_LEN
, ident
))
2048 src
-> mType
= MatchElementType (ident
);
2049 if (src
-> mType
== ET_NONE
) {
2050 TrErrorAt (tr
, line
, col
, "Expected a source element type.\n");
2053 if ((src
-> mFormat
== SF_BIN_LE
) || (src
-> mFormat
== SF_BIN_BE
)) {
2054 if (! TrReadOperator (tr
, ","))
2056 if (src
-> mType
== ET_INT
) {
2057 if (! TrReadInt (tr
, MIN_BIN_SIZE
, MAX_BIN_SIZE
, & intVal
))
2059 src
-> mSize
= intVal
;
2060 if (TrIsOperator (tr
, ",")) {
2061 TrReadOperator (tr
, ",");
2062 TrIndication (tr
, & line
, & col
);
2063 if (! TrReadInt (tr
, 0x80000000, 0x7FFFFFFF, & intVal
))
2065 if ((abs (intVal
) < MIN_BIN_BITS
) || (abs (intVal
) > (8 * src
-> mSize
))) {
2066 TrErrorAt (tr
, line
, col
, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS
, 8 * src
-> mSize
);
2069 src
-> mBits
= intVal
;
2071 src
-> mBits
= 8 * src
-> mSize
;
2074 TrIndication (tr
, & line
, & col
);
2075 if (! TrReadInt (tr
, 0x80000000, 0x7FFFFFFF, & intVal
))
2077 if ((intVal
!= 4) && (intVal
!= 8)) {
2078 TrErrorAt (tr
, line
, col
, "Expected a value of 4 or 8.\n");
2081 src
-> mSize
= (uint
) intVal
;
2084 } else if ((src
-> mFormat
== SF_ASCII
) && (src
-> mType
== ET_INT
)) {
2085 if (! TrReadOperator (tr
, ","))
2087 if (! TrReadInt (tr
, MIN_ASCII_BITS
, MAX_ASCII_BITS
, & intVal
))
2090 src
-> mBits
= intVal
;
2095 if (TrIsOperator (tr
, ";")) {
2096 TrReadOperator (tr
, ";");
2097 if (! TrReadInt (tr
, 0, 0x7FFFFFFF, & intVal
))
2099 src
-> mSkip
= (uint
) intVal
;
2104 if (! TrReadOperator (tr
, ")"))
2106 if (TrIsOperator (tr
, "@")) {
2107 TrReadOperator (tr
, "@");
2108 if (! TrReadInt (tr
, 0, 0x7FFFFFFF, & intVal
))
2110 src
-> mOffset
= (uint
) intVal
;
2114 if (! TrReadOperator (tr
, ":"))
2116 if (! TrReadString (tr
, MAX_PATH_LEN
, src
-> mPath
))
2121 // Process the list of sources in the data set definition.
2122 static int ProcessSources (TokenReaderT
* tr
, HrirDataT
* hData
) {
2123 uint
* setCount
= NULL
, * setFlag
= NULL
;
2124 double * hrir
= NULL
;
2125 uint line
, col
, ei
, ai
;
2129 setCount
= (uint
*) calloc (hData
-> mEvCount
, sizeof (uint
));
2130 setFlag
= (uint
*) calloc (hData
-> mIrCount
, sizeof (uint
));
2131 hrir
= CreateArray (hData
-> mIrPoints
);
2132 while (TrIsOperator (tr
, "[")) {
2133 TrIndication (tr
, & line
, & col
);
2134 TrReadOperator (tr
, "[");
2135 if (ReadIndexPair (tr
, hData
, & ei
, & ai
)) {
2136 if (TrReadOperator (tr
, "]")) {
2137 if (! setFlag
[hData
-> mEvOffset
[ei
] + ai
]) {
2138 if (TrReadOperator (tr
, "=")) {
2141 if (ReadSourceRef (tr
, & src
)) {
2142 if (LoadSource (& src
, hData
-> mIrRate
, hData
-> mIrPoints
, hrir
)) {
2143 AverageHrirMagnitude (hrir
, 1.0 / factor
, ei
, ai
, hData
);
2145 if (! TrIsOperator (tr
, "+"))
2147 TrReadOperator (tr
, "+");
2151 DestroyArray (hrir
);
2156 setFlag
[hData
-> mEvOffset
[ei
] + ai
] = 1;
2161 TrErrorAt (tr
, line
, col
, "Redefinition of source.\n");
2165 DestroyArray (hrir
);
2171 while ((ei
< hData
-> mEvCount
) && (setCount
[ei
] < 1))
2173 if (ei
< hData
-> mEvCount
) {
2174 hData
-> mEvStart
= ei
;
2175 while ((ei
< hData
-> mEvCount
) && (setCount
[ei
] == hData
-> mAzCount
[ei
]))
2177 if (ei
>= hData
-> mEvCount
) {
2178 if (! TrLoad (tr
)) {
2179 DestroyArray (hrir
);
2184 TrError (tr
, "Errant data at end of source list.\n");
2187 TrError (tr
, "Missing sources for elevation index %d.\n", ei
);
2190 TrError (tr
, "Missing source references.\n");
2192 DestroyArray (hrir
);
2198 /* Parse the data set definition and process the source data, storing the
2199 * resulting data set as desired. If the input name is NULL it will read
2200 * from standard input.
2202 static int ProcessDefinition (const char * inName
, const size_t fftSize
, const int equalize
, const int surface
, const double limit
, const size_t truncSize
, const OutputFormatT outFormat
, const char * outName
) {
2206 double * dfa
= NULL
;
2207 char rateStr
[8 + 1], expName
[MAX_PATH_LEN
];
2210 hData
.mIrPoints
= 0;
2216 hData
.mDistance
= 0;
2218 fprintf (stdout
, "Reading HRIR definition...\n");
2219 if (inName
!= NULL
) {
2220 fp
= fopen (inName
, "r");
2222 fprintf (stderr
, "Error: Could not open definition file '%s'\n", inName
);
2225 TrSetup (fp
, inName
, & tr
);
2228 TrSetup (fp
, "<stdin>", & tr
);
2230 if (! ProcessMetrics (& tr
, fftSize
, truncSize
, & hData
)) {
2235 hData
. mHrirs
= CreateArray (hData
. mIrCount
* hData
. mIrSize
);
2236 hData
. mHrtds
= CreateArray (hData
. mIrCount
);
2237 if (! ProcessSources (& tr
, & hData
)) {
2238 DestroyArray (hData
. mHrtds
);
2239 DestroyArray (hData
. mHrirs
);
2247 dfa
= CreateArray (1 + (hData
. mFftSize
/ 2));
2248 fprintf (stdout
, "Calculating diffuse-field average...\n");
2249 CalculateDiffuseFieldAverage (& hData
, surface
, limit
, dfa
);
2250 fprintf (stdout
, "Performing diffuse-field equalization...\n");
2251 DiffuseFieldEqualize (dfa
, & hData
);
2254 fprintf (stdout
, "Performing minimum phase reconstruction...\n");
2255 ReconstructHrirs (& hData
);
2256 fprintf (stdout
, "Truncating minimum-phase HRIRs...\n");
2257 hData
. mIrPoints
= truncSize
;
2258 fprintf (stdout
, "Synthesizing missing elevations...\n");
2259 SynthesizeHrirs (& hData
);
2260 fprintf (stdout
, "Normalizing final HRIRs...\n");
2261 NormalizeHrirs (& hData
);
2262 fprintf (stdout
, "Calculating impulse delays...\n");
2263 CalculateHrtds (& hData
);
2264 snprintf (rateStr
, 8, "%u", hData
. mIrRate
);
2265 StrSubst (outName
, "%r", rateStr
, MAX_PATH_LEN
, expName
);
2266 switch (outFormat
) {
2268 fprintf (stdout
, "Creating MHR data set file...\n");
2269 if (! StoreMhr (& hData
, expName
))
2273 fprintf (stderr
, "Creating OpenAL Soft table file...\n");
2274 if (! StoreTable (& hData
, expName
))
2280 DestroyArray (hData
. mHrtds
);
2281 DestroyArray (hData
. mHrirs
);
2285 // Standard command line dispatch.
2286 int main (const int argc
, const char * argv
[]) {
2287 const char * inName
= NULL
, * outName
= NULL
;
2288 OutputFormatT outFormat
;
2291 int equalize
, surface
;
2297 fprintf (stderr
, "Error: No command specified. See '%s -h' for help.\n", argv
[0]);
2300 if ((strcmp (argv
[1], "--help") == 0) || (strcmp (argv
[1], "-h") == 0)) {
2301 fprintf (stdout
, "HRTF Processing and Composition Utility\n\n");
2302 fprintf (stdout
, "Usage: %s <command> [<option>...]\n\n", argv
[0]);
2303 fprintf (stdout
, "Commands:\n");
2304 fprintf (stdout
, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n");
2305 fprintf (stdout
, " Defaults output to: ./oalsoft_hrtf_%%r.mhr\n");
2306 fprintf (stdout
, " -t, --make-tab Makes the built-in table used when compiling OpenAL Soft.\n");
2307 fprintf (stdout
, " Defaults output to: ./hrtf_tables.inc\n");
2308 fprintf (stdout
, " -h, --help Displays this help information.\n\n");
2309 fprintf (stdout
, "Options:\n");
2310 fprintf (stdout
, " -f=<points> Override the FFT window size (defaults to the first power-\n");
2311 fprintf (stdout
, " of-two that fits four times the number of HRIR points).\n");
2312 fprintf (stdout
, " -e={on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE
? "on" : "off"));
2313 fprintf (stdout
, " -s={on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE
? "on" : "off"));
2314 fprintf (stdout
, " -l={<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
2315 fprintf (stdout
, " average (default: %.2f).\n", DEFAULT_LIMIT
);
2316 fprintf (stdout
, " -w=<points> Specify the size of the truncation window that's applied\n");
2317 fprintf (stdout
, " after minimum-phase reconstruction (default: %d).\n", DEFAULT_TRUNCSIZE
);
2318 fprintf (stdout
, " -i=<filename> Specify an HRIR definition file to use (defaults to stdin).\n");
2319 fprintf (stdout
, " -o=<filename> Specify an output file. Overrides command-selected default.\n");
2320 fprintf (stdout
, " Use of '%%r' will be substituted with the data set sample rate.\n");
2323 if ((strcmp (argv
[1], "--make-mhr") == 0) || (strcmp (argv
[1], "-m") == 0)) {
2327 outName
= "./oalsoft_hrtf_%r.mhr";
2329 } else if ((strcmp (argv
[1], "--make-tab") == 0) || (strcmp (argv
[1], "-t") == 0)) {
2333 outName
= "./hrtf_tables.inc";
2334 outFormat
= OF_TABLE
;
2336 fprintf (stderr
, "Error: Invalid command '%s'.\n", argv
[1]);
2341 equalize
= DEFAULT_EQUALIZE
;
2342 surface
= DEFAULT_SURFACE
;
2343 limit
= DEFAULT_LIMIT
;
2344 truncSize
= DEFAULT_TRUNCSIZE
;
2345 while (argi
< argc
) {
2346 if (strncmp (argv
[argi
], "-f=", 3) == 0) {
2347 fftSize
= strtoul (& argv
[argi
] [3], & end
, 10);
2348 if ((end
[0] != '\0') || (fftSize
& (fftSize
- 1)) || (fftSize
< MIN_FFTSIZE
) || (fftSize
> MAX_FFTSIZE
)) {
2349 fprintf (stderr
, "Error: Expected a power-of-two value from %d to %d for '-f'.\n", MIN_FFTSIZE
, MAX_FFTSIZE
);
2352 } else if (strncmp (argv
[argi
], "-e=", 3) == 0) {
2353 if (strcmp (& argv
[argi
] [3], "on") == 0) {
2355 } else if (strcmp (& argv
[argi
] [3], "off") == 0) {
2358 fprintf (stderr
, "Error: Expected 'on' or 'off' for '-e'.\n");
2361 } else if (strncmp (argv
[argi
], "-s=", 3) == 0) {
2362 if (strcmp (& argv
[argi
] [3], "on") == 0) {
2364 } else if (strcmp (& argv
[argi
] [3], "off") == 0) {
2367 fprintf (stderr
, "Error: Expected 'on' or 'off' for '-s'.\n");
2370 } else if (strncmp (argv
[argi
], "-l=", 3) == 0) {
2371 if (strcmp (& argv
[argi
] [3], "none") == 0) {
2374 limit
= strtod (& argv
[argi
] [3], & end
);
2375 if ((end
[0] != '\0') || (limit
< MIN_LIMIT
) || (limit
> MAX_LIMIT
)) {
2376 fprintf (stderr
, "Error: Expected 'none' or a value from %.2f to %.2f for '-l'.\n", MIN_LIMIT
, MAX_LIMIT
);
2380 } else if (strncmp (argv
[argi
], "-w=", 3) == 0) {
2381 truncSize
= strtoul (& argv
[argi
] [3], & end
, 10);
2382 if ((end
[0] != '\0') || (truncSize
< MIN_TRUNCSIZE
) || (truncSize
> MAX_TRUNCSIZE
) || (truncSize
% MOD_TRUNCSIZE
)) {
2383 fprintf (stderr
, "Error: Expected a value from %d to %d in multiples of %d for '-w'.\n", MIN_TRUNCSIZE
, MAX_TRUNCSIZE
, MOD_TRUNCSIZE
);
2386 } else if (strncmp (argv
[argi
], "-i=", 3) == 0) {
2387 inName
= & argv
[argi
] [3];
2388 } else if (strncmp (argv
[argi
], "-o=", 3) == 0) {
2389 outName
= & argv
[argi
] [3];
2391 fprintf (stderr
, "Error: Invalid option '%s'.\n", argv
[argi
]);
2396 if (! ProcessDefinition (inName
, fftSize
, equalize
, surface
, limit
, truncSize
, outFormat
, outName
))
2398 fprintf (stdout
, "Operation completed.\n");