2 * HRTF utility for producing and demonstrating the process of creating an
3 * OpenAL Soft compatible HRIR data set.
5 * Copyright (C) 2011-2014 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
52 * The formulae for calculating the Kaiser window metrics are from the
55 * Discrete-Time Signal Processing
56 * Alan V. Oppenheim and Ronald W. Schafer
57 * Prentice-Hall Signal Processing Series
61 // Needed for 64-bit unsigned integer.
71 // Rely (if naively) on OpenAL's header for the types used for serialization.
76 #define M_PI (3.14159265358979323846)
80 #define HUGE_VAL (1.0 / 0.0)
83 // The epsilon used to maintain signal stability.
84 #define EPSILON (1e-15)
86 // Constants for accessing the token reader's ring buffer.
87 #define TR_RING_BITS (16)
88 #define TR_RING_SIZE (1 << TR_RING_BITS)
89 #define TR_RING_MASK (TR_RING_SIZE - 1)
91 // The token reader's load interval in bytes.
92 #define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
94 // The maximum identifier length used when processing the data set
96 #define MAX_IDENT_LEN (16)
98 // The maximum path length used when processing filenames.
99 #define MAX_PATH_LEN (256)
101 // The limits for the sample 'rate' metric in the data set definition and for
103 #define MIN_RATE (32000)
104 #define MAX_RATE (96000)
106 // The limits for the HRIR 'points' metric in the data set definition.
107 #define MIN_POINTS (16)
108 #define MAX_POINTS (8192)
110 // The limits to the number of 'azimuths' listed in the data set definition.
111 #define MIN_EV_COUNT (5)
112 #define MAX_EV_COUNT (128)
114 // The limits for each of the 'azimuths' listed in the data set definition.
115 #define MIN_AZ_COUNT (1)
116 #define MAX_AZ_COUNT (128)
118 // The limits for the listener's head 'radius' in the data set definition.
119 #define MIN_RADIUS (0.05)
120 #define MAX_RADIUS (0.15)
122 // The limits for the 'distance' from source to listener in the definition
124 #define MIN_DISTANCE (0.5)
125 #define MAX_DISTANCE (2.5)
127 // The maximum number of channels that can be addressed for a WAVE file
128 // source listed in the data set definition.
129 #define MAX_WAVE_CHANNELS (65535)
131 // The limits to the byte size for a binary source listed in the definition
133 #define MIN_BIN_SIZE (2)
134 #define MAX_BIN_SIZE (4)
136 // The minimum number of significant bits for binary sources listed in the
137 // data set definition. The maximum is calculated from the byte size.
138 #define MIN_BIN_BITS (16)
140 // The limits to the number of significant bits for an ASCII source listed in
141 // the data set definition.
142 #define MIN_ASCII_BITS (16)
143 #define MAX_ASCII_BITS (32)
145 // The limits to the FFT window size override on the command line.
146 #define MIN_FFTSIZE (512)
147 #define MAX_FFTSIZE (16384)
149 // The limits to the equalization range limit on the command line.
150 #define MIN_LIMIT (2.0)
151 #define MAX_LIMIT (120.0)
153 // The limits to the truncation window size on the command line.
154 #define MIN_TRUNCSIZE (8)
155 #define MAX_TRUNCSIZE (128)
157 // The limits to the custom head radius on the command line.
158 #define MIN_CUSTOM_RADIUS (0.05)
159 #define MAX_CUSTOM_RADIUS (0.15)
161 // The truncation window size must be a multiple of the below value to allow
162 // for vectorized convolution.
163 #define MOD_TRUNCSIZE (8)
165 // The defaults for the command line options.
166 #define DEFAULT_EQUALIZE (1)
167 #define DEFAULT_SURFACE (1)
168 #define DEFAULT_LIMIT (24.0)
169 #define DEFAULT_TRUNCSIZE (32)
170 #define DEFAULT_HEAD_MODEL (HM_DATASET)
171 #define DEFAULT_CUSTOM_RADIUS (0.0)
173 // The four-character-codes for RIFF/RIFX WAVE file chunks.
174 #define FOURCC_RIFF (0x46464952) // 'RIFF'
175 #define FOURCC_RIFX (0x58464952) // 'RIFX'
176 #define FOURCC_WAVE (0x45564157) // 'WAVE'
177 #define FOURCC_FMT (0x20746D66) // 'fmt '
178 #define FOURCC_DATA (0x61746164) // 'data'
179 #define FOURCC_LIST (0x5453494C) // 'LIST'
180 #define FOURCC_WAVL (0x6C766177) // 'wavl'
181 #define FOURCC_SLNT (0x746E6C73) // 'slnt'
183 // The supported wave formats.
184 #define WAVE_FORMAT_PCM (0x0001)
185 #define WAVE_FORMAT_IEEE_FLOAT (0x0003)
186 #define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
188 // The maximum propagation delay value supported by OpenAL Soft.
189 #define MAX_HRTD (63.0)
191 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
192 // response protocol 01.
193 #define MHR_FORMAT ("MinPHR01")
195 // Byte order for the serialization routines.
202 // Source format for the references listed in the data set definition.
205 SF_WAVE
, // RIFF/RIFX WAVE file.
206 SF_BIN_LE
, // Little-endian binary file.
207 SF_BIN_BE
, // Big-endian binary file.
208 SF_ASCII
// ASCII text file.
211 // Element types for the references listed in the data set definition.
214 ET_INT
, // Integer elements.
215 ET_FP
// Floating-point elements.
218 // Head model used for calculating the impulse delays.
221 HM_DATASET
, // Measure the onset from the dataset.
222 HM_SPHERE
// Calculate the onset using a spherical head model.
225 // Desired output format from the command line.
228 OF_MHR
, // OpenAL Soft MHR data set file.
229 OF_TABLE
// OpenAL Soft built-in table file (used when compiling).
232 // Unsigned integer type.
233 typedef unsigned int uint
;
235 // Serialization types. The trailing digit indicates the number of bits.
236 typedef ALubyte uint8
;
239 typedef ALuint uint32
;
240 typedef ALuint64SOFT uint64
;
242 typedef enum ByteOrderT ByteOrderT
;
243 typedef enum SourceFormatT SourceFormatT
;
244 typedef enum ElementTypeT ElementTypeT
;
245 typedef enum HeadModelT HeadModelT
;
246 typedef enum OutputFormatT OutputFormatT
;
248 typedef struct TokenReaderT TokenReaderT
;
249 typedef struct SourceRefT SourceRefT
;
250 typedef struct HrirDataT HrirDataT
;
251 typedef struct ResamplerT ResamplerT
;
253 // Token reader state for parsing the data set definition.
254 struct TokenReaderT
{
259 char mRing
[TR_RING_SIZE
];
264 // Source reference state used when loading sources.
266 SourceFormatT mFormat
;
273 char mPath
[MAX_PATH_LEN
+ 1];
276 // The HRIR metrics and data set used when loading, processing, and storing
277 // the resulting HRTF.
286 mAzCount
[MAX_EV_COUNT
],
287 mEvOffset
[MAX_EV_COUNT
];
295 // The resampler metrics and FIR filter.
304 /* Token reader routines for parsing text files. Whitespace is not
305 * significant. It can process tokens as identifiers, numbers (integer and
306 * floating-point), strings, and operators. Strings must be encapsulated by
307 * double-quotes and cannot span multiple lines.
310 // Setup the reader on the given file. The filename can be NULL if no error
311 // output is desired.
312 static void TrSetup (FILE * fp
, const char * filename
, TokenReaderT
* tr
) {
313 const char * name
= NULL
;
318 // If a filename was given, store a pointer to the base name.
319 if (filename
!= NULL
) {
320 while ((ch
= (* filename
)) != '\0') {
321 if ((ch
== '/') || (ch
== '\\'))
333 // Prime the reader's ring buffer, and return a result indicating that there
334 // is text to process.
335 static int TrLoad (TokenReaderT
* tr
) {
336 size_t toLoad
, in
, count
;
338 toLoad
= TR_RING_SIZE
- (tr
-> mIn
- tr
-> mOut
);
339 if ((toLoad
>= TR_LOAD_SIZE
) && (! feof (tr
-> mFile
))) {
340 // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
341 toLoad
= TR_LOAD_SIZE
;
342 in
= tr
-> mIn
& TR_RING_MASK
;
343 count
= TR_RING_SIZE
- in
;
344 if (count
< toLoad
) {
345 tr
-> mIn
+= fread (& tr
-> mRing
[in
], 1, count
, tr
-> mFile
);
346 tr
-> mIn
+= fread (& tr
-> mRing
[0], 1, toLoad
- count
, tr
-> mFile
);
348 tr
-> mIn
+= fread (& tr
-> mRing
[in
], 1, toLoad
, tr
-> mFile
);
350 if (tr
-> mOut
>= TR_RING_SIZE
) {
351 tr
-> mOut
-= TR_RING_SIZE
;
352 tr
-> mIn
-= TR_RING_SIZE
;
355 if (tr
-> mIn
> tr
-> mOut
)
360 // Error display routine. Only displays when the base name is not NULL.
361 static void TrErrorVA (const TokenReaderT
* tr
, uint line
, uint column
, const char * format
, va_list argPtr
) {
362 if (tr
-> mName
!= NULL
) {
363 fprintf (stderr
, "Error (%s:%u:%u): ", tr
-> mName
, line
, column
);
364 vfprintf (stderr
, format
, argPtr
);
368 // Used to display an error at a saved line/column.
369 static void TrErrorAt (const TokenReaderT
* tr
, uint line
, uint column
, const char * format
, ...) {
372 va_start (argPtr
, format
);
373 TrErrorVA (tr
, line
, column
, format
, argPtr
);
377 // Used to display an error at the current line/column.
378 static void TrError (const TokenReaderT
* tr
, const char * format
, ...) {
381 va_start (argPtr
, format
);
382 TrErrorVA (tr
, tr
-> mLine
, tr
-> mColumn
, format
, argPtr
);
386 // Skips to the next line.
387 static void TrSkipLine (TokenReaderT
* tr
) {
390 while (TrLoad (tr
)) {
391 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
402 // Skips to the next token.
403 static int TrSkipWhitespace (TokenReaderT
* tr
) {
406 while (TrLoad (tr
)) {
407 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
416 } else if (ch
== '#') {
425 // Get the line and/or column of the next token (or the end of input).
426 static void TrIndication (TokenReaderT
* tr
, uint
* line
, uint
* column
) {
427 TrSkipWhitespace (tr
);
429 (* line
) = tr
-> mLine
;
431 (* column
) = tr
-> mColumn
;
434 // Checks to see if a token is the given operator. It does not display any
435 // errors and will not proceed to the next token.
436 static int TrIsOperator (TokenReaderT
* tr
, const char * op
) {
440 if (! TrSkipWhitespace (tr
))
444 while ((op
[len
] != '\0') && (out
< tr
-> mIn
)) {
445 ch
= tr
-> mRing
[out
& TR_RING_MASK
];
451 if (op
[len
] == '\0')
456 /* The TrRead*() routines obtain the value of a matching token type. They
457 * display type, form, and boundary errors and will proceed to the next
461 // Reads and validates an identifier token.
462 static int TrReadIdent (TokenReaderT
* tr
, const uint maxLen
, char * ident
) {
467 if (TrSkipWhitespace (tr
)) {
469 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
470 if ((ch
== '_') || isalpha (ch
)) {
479 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
480 } while ((ch
== '_') || isdigit (ch
) || isalpha (ch
));
481 tr
-> mColumn
+= len
;
483 TrErrorAt (tr
, tr
-> mLine
, col
, "Identifier is too long.\n");
490 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected an identifier.\n");
494 // Reads and validates (including bounds) an integer token.
495 static int TrReadInt (TokenReaderT
* tr
, const int loBound
, const int hiBound
, int * value
) {
496 uint col
, digis
, len
;
497 char ch
, temp
[64 + 1];
500 if (TrSkipWhitespace (tr
)) {
503 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
504 if ((ch
== '+') || (ch
== '-')) {
510 while (TrLoad (tr
)) {
511 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
520 tr
-> mColumn
+= len
;
521 if ((digis
> 0) && (ch
!= '.') && (! isalpha (ch
))) {
523 TrErrorAt (tr
, tr
-> mLine
, col
, "Integer is too long.");
527 (* value
) = strtol (temp
, NULL
, 10);
528 if (((* value
) < loBound
) || ((* value
) > hiBound
)) {
529 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a value from %d to %d.\n", loBound
, hiBound
);
535 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected an integer.\n");
539 // Reads and validates (including bounds) a float token.
540 static int TrReadFloat (TokenReaderT
* tr
, const double loBound
, const double hiBound
, double * value
) {
541 uint col
, digis
, len
;
542 char ch
, temp
[64 + 1];
545 if (TrSkipWhitespace (tr
)) {
548 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
549 if ((ch
== '+') || (ch
== '-')) {
555 while (TrLoad (tr
)) {
556 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
571 while (TrLoad (tr
)) {
572 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
582 if ((ch
== 'E') || (ch
== 'e')) {
588 if ((ch
== '+') || (ch
== '-')) {
594 while (TrLoad (tr
)) {
595 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
605 tr
-> mColumn
+= len
;
606 if ((digis
> 0) && (ch
!= '.') && (! isalpha (ch
))) {
608 TrErrorAt (tr
, tr
-> mLine
, col
, "Float is too long.");
612 (* value
) = strtod (temp
, NULL
);
613 if (((* value
) < loBound
) || ((* value
) > hiBound
)) {
614 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a value from %f to %f.\n", loBound
, hiBound
);
620 tr
-> mColumn
+= len
;
623 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a float.\n");
627 // Reads and validates a string token.
628 static int TrReadString (TokenReaderT
* tr
, const uint maxLen
, char * text
) {
633 if (TrSkipWhitespace (tr
)) {
635 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
639 while (TrLoad (tr
)) {
640 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
645 TrErrorAt (tr
, tr
-> mLine
, col
, "Unterminated string at end of line.\n");
653 tr
-> mColumn
+= 1 + len
;
654 TrErrorAt (tr
, tr
-> mLine
, col
, "Unterminated string at end of input.\n");
657 tr
-> mColumn
+= 2 + len
;
659 TrErrorAt (tr
, tr
-> mLine
, col
, "String is too long.\n");
666 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected a string.\n");
670 // Reads and validates the given operator.
671 static int TrReadOperator (TokenReaderT
* tr
, const char * op
) {
676 if (TrSkipWhitespace (tr
)) {
679 while ((op
[len
] != '\0') && TrLoad (tr
)) {
680 ch
= tr
-> mRing
[tr
-> mOut
& TR_RING_MASK
];
686 tr
-> mColumn
+= len
;
687 if (op
[len
] == '\0')
690 TrErrorAt (tr
, tr
-> mLine
, col
, "Expected '%s' operator.\n", op
);
694 /* Performs a string substitution. Any case-insensitive occurrences of the
695 * pattern string are replaced with the replacement string. The result is
696 * truncated if necessary.
698 static int StrSubst (const char * in
, const char * pat
, const char * rep
, const size_t maxLen
, char * out
) {
699 size_t inLen
, patLen
, repLen
;
704 patLen
= strlen (pat
);
705 repLen
= strlen (rep
);
709 while ((si
< inLen
) && (di
< maxLen
)) {
710 if (patLen
<= (inLen
- si
)) {
711 if (strncasecmp (& in
[si
], pat
, patLen
) == 0) {
712 if (repLen
> (maxLen
- di
)) {
713 repLen
= maxLen
- di
;
716 strncpy (& out
[di
], rep
, repLen
);
728 return (! truncated
);
731 // Provide missing math routines for MSVC versions < 1800 (Visual Studio 2013).
732 #if defined(_MSC_VER) && _MSC_VER < 1800
733 static double round (double val
) {
735 return (ceil (val
- 0.5));
736 return (floor (val
+ 0.5));
739 static double fmin (double a
, double b
) {
740 return ((a
< b
) ? a
: b
);
743 static double fmax (double a
, double b
) {
744 return ((a
> b
) ? a
: b
);
748 // Simple clamp routine.
749 static double Clamp (const double val
, const double lower
, const double upper
) {
750 return (fmin (fmax (val
, lower
), upper
));
753 // Performs linear interpolation.
754 static double Lerp (const double a
, const double b
, const double f
) {
755 return (a
+ (f
* (b
- a
)));
758 // Performs a high-passed triangular probability density function dither from
759 // a double to an integer. It assumes the input sample is already scaled.
760 static int HpTpdfDither (const double in
, int * hpHist
) {
761 const double PRNG_SCALE
= 1.0 / (RAND_MAX
+ 1.0);
766 out
= round (in
+ (PRNG_SCALE
* (prn
- (* hpHist
))));
771 // Allocates an array of doubles.
772 static double *CreateArray(const size_t n
)
776 a
= calloc(n
, sizeof(double));
779 fprintf(stderr
, "Error: Out of memory.\n");
785 // Frees an array of doubles.
786 static void DestroyArray(double *a
)
789 // Complex number routines. All outputs must be non-NULL.
791 // Magnitude/absolute value.
792 static double ComplexAbs (const double r
, const double i
) {
793 return (sqrt ((r
* r
) + (i
* i
)));
797 static void ComplexMul (const double aR
, const double aI
, const double bR
, const double bI
, double * outR
, double * outI
) {
798 (* outR
) = (aR
* bR
) - (aI
* bI
);
799 (* outI
) = (aI
* bR
) + (aR
* bI
);
803 static void ComplexExp (const double inR
, const double inI
, double * outR
, double * outI
) {
807 (* outR
) = e
* cos (inI
);
808 (* outI
) = e
* sin (inI
);
811 /* Fast Fourier transform routines. The number of points must be a power of
812 * two. In-place operation is possible only if both the real and imaginary
813 * parts are in-place together.
816 // Performs bit-reversal ordering.
817 static void FftArrange (const uint n
, const double * inR
, const double * inI
, double * outR
, double * outI
) {
821 if ((inR
== outR
) && (inI
== outI
)) {
822 // Handle in-place arrangement.
824 for (k
= 0; k
< n
; k
++) {
834 while (rk
& (m
>>= 1))
839 // Handle copy arrangement.
841 for (k
= 0; k
< n
; k
++) {
845 while (rk
& (m
>>= 1))
852 // Performs the summation.
853 static void FftSummation (const uint n
, const double s
, double * re
, double * im
) {
856 double vR
, vI
, wR
, wI
;
861 for (m
= 1, m2
= 2; m
< n
; m
<<= 1, m2
<<= 1) {
862 // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
863 vR
= sin (0.5 * pi
/ m
);
866 // w = Complex (1.0, 0.0)
869 for (i
= 0; i
< m
; i
++) {
870 for (k
= i
; k
< n
; k
+= m2
) {
872 // t = ComplexMul (w, out [km2])
873 tR
= (wR
* re
[mk
]) - (wI
* im
[mk
]);
874 tI
= (wR
* im
[mk
]) + (wI
* re
[mk
]);
875 // out [mk] = ComplexSub (out [k], t)
876 re
[mk
] = re
[k
] - tR
;
877 im
[mk
] = im
[k
] - tI
;
878 // out [k] = ComplexAdd (out [k], t)
882 // t = ComplexMul (v, w)
883 tR
= (vR
* wR
) - (vI
* wI
);
884 tI
= (vR
* wI
) + (vI
* wR
);
885 // w = ComplexAdd (w, t)
892 // Performs a forward FFT.
893 static void FftForward (const uint n
, const double * inR
, const double * inI
, double * outR
, double * outI
) {
894 FftArrange (n
, inR
, inI
, outR
, outI
);
895 FftSummation (n
, 1.0, outR
, outI
);
898 // Performs an inverse FFT.
899 static void FftInverse (const uint n
, const double * inR
, const double * inI
, double * outR
, double * outI
) {
903 FftArrange (n
, inR
, inI
, outR
, outI
);
904 FftSummation (n
, -1.0, outR
, outI
);
906 for (i
= 0; i
< n
; i
++) {
912 /* Calculate the complex helical sequence (or discrete-time analytical
913 * signal) of the given input using the Hilbert transform. Given the
914 * negative natural logarithm of a signal's magnitude response, the imaginary
915 * components can be used as the angles for minimum-phase reconstruction.
917 static void Hilbert (const uint n
, const double * in
, double * outR
, double * outI
) {
921 // Handle in-place operation.
922 for (i
= 0; i
< n
; i
++)
925 // Handle copy operation.
926 for (i
= 0; i
< n
; i
++) {
931 FftForward (n
, outR
, outI
, outR
, outI
);
932 /* Currently the Fourier routines operate only on point counts that are
933 * powers of two. If that changes and n is odd, the following conditional
934 * should be: i < (n + 1) / 2.
936 for (i
= 1; i
< (n
/ 2); i
++) {
940 // If n is odd, the following increment should be skipped.
942 for (; i
< n
; i
++) {
946 FftInverse (n
, outR
, outI
, outR
, outI
);
949 /* Calculate the magnitude response of the given input. This is used in
950 * place of phase decomposition, since the phase residuals are discarded for
951 * minimum phase reconstruction. The mirrored half of the response is also
954 static void MagnitudeResponse (const uint n
, const double * inR
, const double * inI
, double * out
) {
955 const uint m
= 1 + (n
/ 2);
958 for (i
= 0; i
< m
; i
++)
959 out
[i
] = fmax (ComplexAbs (inR
[i
], inI
[i
]), EPSILON
);
962 /* Apply a range limit (in dB) to the given magnitude response. This is used
963 * to adjust the effects of the diffuse-field average on the equalization
966 static void LimitMagnitudeResponse (const uint n
, const double limit
, const double * in
, double * out
) {
967 const uint m
= 1 + (n
/ 2);
969 uint i
, lower
, upper
;
972 halfLim
= limit
/ 2.0;
973 // Convert the response to dB.
974 for (i
= 0; i
< m
; i
++)
975 out
[i
] = 20.0 * log10 (in
[i
]);
976 // Use six octaves to calculate the average magnitude of the signal.
977 lower
= ((uint
) ceil (n
/ pow (2.0, 8.0))) - 1;
978 upper
= ((uint
) floor (n
/ pow (2.0, 2.0))) - 1;
980 for (i
= lower
; i
<= upper
; i
++)
982 ave
/= upper
- lower
+ 1;
983 // Keep the response within range of the average magnitude.
984 for (i
= 0; i
< m
; i
++)
985 out
[i
] = Clamp (out
[i
], ave
- halfLim
, ave
+ halfLim
);
986 // Convert the response back to linear magnitude.
987 for (i
= 0; i
< m
; i
++)
988 out
[i
] = pow (10.0, out
[i
] / 20.0);
991 /* Reconstructs the minimum-phase component for the given magnitude response
992 * of a signal. This is equivalent to phase recomposition, sans the missing
993 * residuals (which were discarded). The mirrored half of the response is
996 static void MinimumPhase (const uint n
, const double * in
, double * outR
, double * outI
) {
997 const uint m
= 1 + (n
/ 2);
998 double * mags
= NULL
;
1002 mags
= CreateArray (n
);
1003 for (i
= 0; i
< m
; i
++) {
1004 mags
[i
] = fmax (in
[i
], EPSILON
);
1005 outR
[i
] = -log (mags
[i
]);
1007 for (; i
< n
; i
++) {
1008 mags
[i
] = mags
[n
- i
];
1009 outR
[i
] = outR
[n
- i
];
1011 Hilbert (n
, outR
, outR
, outI
);
1012 // Remove any DC offset the filter has.
1015 for (i
= 1; i
< n
; i
++) {
1016 ComplexExp (0.0, outI
[i
], & aR
, & aI
);
1017 ComplexMul (mags
[i
], 0.0, aR
, aI
, & outR
[i
], & outI
[i
]);
1019 DestroyArray (mags
);
1022 /* This is the normalized cardinal sine (sinc) function.
1024 * sinc(x) = { 1, x = 0
1025 * { sin(pi x) / (pi x), otherwise.
1027 static double Sinc (const double x
) {
1028 if (fabs (x
) < EPSILON
)
1030 return (sin (M_PI
* x
) / (M_PI
* x
));
1033 /* The zero-order modified Bessel function of the first kind, used for the
1036 * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
1037 * = sum_{k=0}^inf ((x / 2)^k / k!)^2
1039 static double BesselI_0 (const double x
) {
1040 double term
, sum
, x2
, y
, last_sum
;
1043 // Start at k=1 since k=0 is trivial.
1048 // Let the integration converge until the term of the sum is no longer
1056 } while (sum
!= last_sum
);
1060 /* Calculate a Kaiser window from the given beta value and a normalized k
1063 * w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
1066 * Where k can be calculated as:
1068 * k = i / l, where -l <= i <= l.
1072 * k = 2 i / M - 1, where 0 <= i <= M.
1074 static double Kaiser (const double b
, const double k
) {
1077 k2
= Clamp (k
, -1.0, 1.0);
1078 if ((k
< -1.0) || (k
> 1.0))
1081 return (BesselI_0 (b
* sqrt (1.0 - k2
)) / BesselI_0 (b
));
1084 // Calculates the greatest common divisor of a and b.
1085 static uint
Gcd (const uint a
, const uint b
) {
1098 /* Calculates the size (order) of the Kaiser window. Rejection is in dB and
1099 * the transition width is normalized frequency (0.5 is nyquist).
1101 * M = { ceil((r - 7.95) / (2.285 2 pi f_t)), r > 21
1102 * { ceil(5.79 / 2 pi f_t), r <= 21.
1105 static uint
CalcKaiserOrder (const double rejection
, const double transition
) {
1108 w_t
= 2.0 * M_PI
* transition
;
1109 if (rejection
> 21.0)
1110 return ((uint
) ceil ((rejection
- 7.95) / (2.285 * w_t
)));
1111 return ((uint
) ceil (5.79 / w_t
));
1114 // Calculates the beta value of the Kaiser window. Rejection is in dB.
1115 static double CalcKaiserBeta (const double rejection
) {
1116 if (rejection
> 50.0)
1117 return (0.1102 * (rejection
- 8.7));
1118 else if (rejection
>= 21.0)
1119 return ((0.5842 * pow (rejection
- 21.0, 0.4)) +
1120 (0.07886 * (rejection
- 21.0)));
1125 /* Calculates a point on the Kaiser-windowed sinc filter for the given half-
1126 * width, beta, gain, and cutoff. The point is specified in non-normalized
1127 * samples, from 0 to M, where M = (2 l + 1).
1129 * w(k) 2 p f_t sinc(2 f_t x)
1131 * x -- centered sample index (i - l)
1132 * k -- normalized and centered window index (x / l)
1133 * w(k) -- window function (Kaiser)
1134 * p -- gain compensation factor when sampling
1135 * f_t -- normalized center frequency (or cutoff; 0.5 is nyquist)
1137 static double SincFilter (const int l
, const double b
, const double gain
, const double cutoff
, const int i
) {
1138 return (Kaiser (b
, ((double) (i
- l
)) / l
) * 2.0 * gain
* cutoff
* Sinc (2.0 * cutoff
* (i
- l
)));
1141 /* This is a polyphase sinc-filtered resampler.
1143 * Upsample Downsample
1145 * p/q = 3/2 p/q = 3/5
1147 * M-+-+-+-> M-+-+-+->
1148 * -------------------+ ---------------------+
1149 * p s * f f f f|f| | p s * f f f f f |
1150 * | 0 * 0 0 0|0|0 | | 0 * 0 0 0 0|0| |
1151 * v 0 * 0 0|0|0 0 | v 0 * 0 0 0|0|0 |
1152 * s * f|f|f f f | s * f f|f|f f |
1153 * 0 * |0|0 0 0 0 | 0 * 0|0|0 0 0 |
1154 * --------+=+--------+ 0 * |0|0 0 0 0 |
1155 * d . d .|d|. d . d ----------+=+--------+
1156 * d . . . .|d|. . . .
1160 * P_f(i,j) = q i mod p + pj
1161 * P_s(i,j) = floor(q i / p) - j
1162 * d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
1163 * { f[P_f(i,j)] s[P_s(i,j)], P_f(i,j) < M
1164 * { 0, P_f(i,j) >= M. }
1167 // Calculate the resampling metrics and build the Kaiser-windowed sinc filter
1168 // that's used to cut frequencies above the destination nyquist.
1169 static void ResamplerSetup (ResamplerT
* rs
, const uint srcRate
, const uint dstRate
) {
1171 double cutoff
, width
, beta
;
1174 gcd
= Gcd (srcRate
, dstRate
);
1175 rs
-> mP
= dstRate
/ gcd
;
1176 rs
-> mQ
= srcRate
/ gcd
;
1177 /* The cutoff is adjusted by half the transition width, so the transition
1178 * ends before the nyquist (0.5). Both are scaled by the downsampling
1181 if (rs
-> mP
> rs
-> mQ
) {
1182 cutoff
= 0.45 / rs
-> mP
;
1183 width
= 0.1 / rs
-> mP
;
1185 cutoff
= 0.45 / rs
-> mQ
;
1186 width
= 0.1 / rs
-> mQ
;
1188 // A rejection of -180 dB is used for the stop band.
1189 l
= CalcKaiserOrder (180.0, width
) / 2;
1190 beta
= CalcKaiserBeta (180.0);
1191 rs
-> mM
= (2 * l
) + 1;
1193 rs
-> mF
= CreateArray (rs
-> mM
);
1194 for (i
= 0; i
< ((int) rs
-> mM
); i
++)
1195 rs
-> mF
[i
] = SincFilter ((int) l
, beta
, rs
-> mP
, cutoff
, i
);
1198 // Clean up after the resampler.
1199 static void ResamplerClear (ResamplerT
* rs
) {
1200 DestroyArray (rs
-> mF
);
1204 // Perform the upsample-filter-downsample resampling operation using a
1205 // polyphase filter implementation.
1206 static void ResamplerRun (ResamplerT
* rs
, const uint inN
, const double * in
, const uint outN
, double * out
) {
1207 const uint p
= rs
-> mP
, q
= rs
-> mQ
, m
= rs
-> mM
, l
= rs
-> mL
;
1208 const double * f
= rs
-> mF
;
1209 double * work
= NULL
;
1217 // Handle in-place operation.
1219 work
= CreateArray (outN
);
1222 // Resample the input.
1223 for (i
= 0; i
< outN
; i
++) {
1225 // Input starts at l to compensate for the filter delay. This will
1226 // drop any build-up from the first half of the filter.
1227 j_f
= (l
+ (q
* i
)) % p
;
1228 j_s
= (l
+ (q
* i
)) / p
;
1230 // Only take input when 0 <= j_s < inN. This single unsigned
1231 // comparison catches both cases.
1233 r
+= f
[j_f
] * in
[j_s
];
1239 // Clean up after in-place operation.
1241 for (i
= 0; i
< outN
; i
++)
1243 DestroyArray (work
);
1247 // Read a binary value of the specified byte order and byte size from a file,
1248 // storing it as a 32-bit unsigned integer.
1249 static int ReadBin4 (FILE * fp
, const char * filename
, const ByteOrderT order
, const uint bytes
, uint32
* out
) {
1254 if (fread (in
, 1, bytes
, fp
) != bytes
) {
1255 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1261 for (i
= 0; i
< bytes
; i
++)
1262 accum
= (accum
<< 8) | in
[bytes
- i
- 1];
1265 for (i
= 0; i
< bytes
; i
++)
1266 accum
= (accum
<< 8) | in
[i
];
1275 // Read a binary value of the specified byte order from a file, storing it as
1276 // a 64-bit unsigned integer.
1277 static int ReadBin8 (FILE * fp
, const char * filename
, const ByteOrderT order
, uint64
* out
) {
1282 if (fread (in
, 1, 8, fp
) != 8) {
1283 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1289 for (i
= 0; i
< 8; i
++)
1290 accum
= (accum
<< 8) | in
[8 - i
- 1];
1293 for (i
= 0; i
< 8; i
++)
1294 accum
= (accum
<< 8) | in
[i
];
1303 // Write an ASCII string to a file.
1304 static int WriteAscii (const char * out
, FILE * fp
, const char * filename
) {
1308 if (fwrite (out
, 1, len
, fp
) != len
) {
1310 fprintf (stderr
, "Error: Bad write to file '%s'.\n", filename
);
1316 // Write a binary value of the given byte order and byte size to a file,
1317 // loading it from a 32-bit unsigned integer.
1318 static int WriteBin4 (const ByteOrderT order
, const uint bytes
, const uint32 in
, FILE * fp
, const char * filename
) {
1324 for (i
= 0; i
< bytes
; i
++)
1325 out
[i
] = (in
>> (i
* 8)) & 0x000000FF;
1328 for (i
= 0; i
< bytes
; i
++)
1329 out
[bytes
- i
- 1] = (in
>> (i
* 8)) & 0x000000FF;
1334 if (fwrite (out
, 1, bytes
, fp
) != bytes
) {
1335 fprintf (stderr
, "Error: Bad write to file '%s'.\n", filename
);
1341 /* Read a binary value of the specified type, byte order, and byte size from
1342 * a file, converting it to a double. For integer types, the significant
1343 * bits are used to normalize the result. The sign of bits determines
1344 * whether they are padded toward the MSB (negative) or LSB (positive).
1345 * Floating-point types are not normalized.
1347 static int ReadBinAsDouble (FILE * fp
, const char * filename
, const ByteOrderT order
, const ElementTypeT type
, const uint bytes
, const int bits
, double * out
) {
1360 if (! ReadBin8 (fp
, filename
, order
, & v8
. ui
))
1365 if (! ReadBin4 (fp
, filename
, order
, bytes
, & v4
. ui
))
1367 if (type
== ET_FP
) {
1368 (* out
) = (double) v4
. f
;
1371 v4
. ui
>>= (8 * bytes
) - ((uint
) bits
);
1373 v4
. ui
&= (0xFFFFFFFF >> (32 + bits
));
1374 if (v4
. ui
& ((uint
) (1 << (abs (bits
) - 1))))
1375 v4
. ui
|= (0xFFFFFFFF << abs (bits
));
1376 (* out
) = v4
. i
/ ((double) (1 << (abs (bits
) - 1)));
1382 /* Read an ascii value of the specified type from a file, converting it to a
1383 * double. For integer types, the significant bits are used to normalize the
1384 * result. The sign of the bits should always be positive. This also skips
1385 * up to one separator character before the element itself.
1387 static int ReadAsciiAsDouble (TokenReaderT
* tr
, const char * filename
, const ElementTypeT type
, const uint bits
, double * out
) {
1390 if (TrIsOperator (tr
, ","))
1391 TrReadOperator (tr
, ",");
1392 else if (TrIsOperator (tr
, ":"))
1393 TrReadOperator (tr
, ":");
1394 else if (TrIsOperator (tr
, ";"))
1395 TrReadOperator (tr
, ";");
1396 else if (TrIsOperator (tr
, "|"))
1397 TrReadOperator (tr
, "|");
1398 if (type
== ET_FP
) {
1399 if (! TrReadFloat (tr
, -HUGE_VAL
, HUGE_VAL
, out
)) {
1400 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1404 if (! TrReadInt (tr
, -(1 << (bits
- 1)), (1 << (bits
- 1)) - 1, & v
)) {
1405 fprintf (stderr
, "Error: Bad read from file '%s'.\n", filename
);
1408 (* out
) = v
/ ((double) ((1 << (bits
- 1)) - 1));
1413 // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
1414 // the source parameters and data set metrics.
1415 static int ReadWaveFormat (FILE * fp
, const ByteOrderT order
, const uint hrirRate
, SourceRefT
* src
) {
1416 uint32 fourCC
, chunkSize
;
1417 uint32 format
, channels
, rate
, dummy
, block
, size
, bits
;
1422 fseek (fp
, (long) chunkSize
, SEEK_CUR
);
1423 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1424 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & chunkSize
)))
1426 } while (fourCC
!= FOURCC_FMT
);
1427 if ((! ReadBin4 (fp
, src
-> mPath
, order
, 2, & format
)) ||
1428 (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & channels
)) ||
1429 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & rate
)) ||
1430 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & dummy
)) ||
1431 (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & block
)))
1434 if (chunkSize
> 14) {
1435 if (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & size
))
1443 if (format
== WAVE_FORMAT_EXTENSIBLE
) {
1444 fseek (fp
, 2, SEEK_CUR
);
1445 if (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & bits
))
1449 fseek (fp
, 4, SEEK_CUR
);
1450 if (! ReadBin4 (fp
, src
-> mPath
, order
, 2, & format
))
1452 fseek (fp
, (long) (chunkSize
- 26), SEEK_CUR
);
1456 fseek (fp
, (long) (chunkSize
- 16), SEEK_CUR
);
1458 fseek (fp
, (long) (chunkSize
- 14), SEEK_CUR
);
1460 if ((format
!= WAVE_FORMAT_PCM
) && (format
!= WAVE_FORMAT_IEEE_FLOAT
)) {
1461 fprintf (stderr
, "Error: Unsupported WAVE format in file '%s'.\n", src
-> mPath
);
1464 if (src
-> mChannel
>= channels
) {
1465 fprintf (stderr
, "Error: Missing source channel in WAVE file '%s'.\n", src
-> mPath
);
1468 if (rate
!= hrirRate
) {
1469 fprintf (stderr
, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src
-> mPath
);
1472 if (format
== WAVE_FORMAT_PCM
) {
1473 if ((size
< 2) || (size
> 4)) {
1474 fprintf (stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
-> mPath
);
1477 if ((bits
< 16) || (bits
> (8 * size
))) {
1478 fprintf (stderr
, "Error: Bad significant bits in WAVE file '%s'.\n", src
-> mPath
);
1481 src
-> mType
= ET_INT
;
1483 if ((size
!= 4) && (size
!= 8)) {
1484 fprintf (stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
-> mPath
);
1487 src
-> mType
= ET_FP
;
1489 src
-> mSize
= size
;
1490 src
-> mBits
= (int) bits
;
1491 src
-> mSkip
= channels
;
1495 // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
1496 static int ReadWaveData (FILE * fp
, const SourceRefT
* src
, const ByteOrderT order
, const uint n
, double * hrir
) {
1497 int pre
, post
, skip
;
1500 pre
= (int) (src
-> mSize
* src
-> mChannel
);
1501 post
= (int) (src
-> mSize
* (src
-> mSkip
- src
-> mChannel
- 1));
1503 for (i
= 0; i
< n
; i
++) {
1506 fseek (fp
, skip
, SEEK_CUR
);
1507 if (! ReadBinAsDouble (fp
, src
-> mPath
, order
, src
-> mType
, src
-> mSize
, src
-> mBits
, & hrir
[i
]))
1512 fseek (fp
, skip
, SEEK_CUR
);
1516 // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
1518 static int ReadWaveList (FILE * fp
, const SourceRefT
* src
, const ByteOrderT order
, const uint n
, double * hrir
) {
1519 uint32 fourCC
, chunkSize
, listSize
, count
;
1520 uint block
, skip
, offset
, i
;
1524 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1525 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & chunkSize
)))
1527 if (fourCC
== FOURCC_DATA
) {
1528 block
= src
-> mSize
* src
-> mSkip
;
1529 count
= chunkSize
/ block
;
1530 if (count
< (src
-> mOffset
+ n
)) {
1531 fprintf (stderr
, "Error: Bad read from file '%s'.\n", src
-> mPath
);
1534 fseek (fp
, (long) (src
-> mOffset
* block
), SEEK_CUR
);
1535 if (! ReadWaveData (fp
, src
, order
, n
, & hrir
[0]))
1538 } else if (fourCC
== FOURCC_LIST
) {
1539 if (! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
))
1542 if (fourCC
== FOURCC_WAVL
)
1546 fseek (fp
, (long) chunkSize
, SEEK_CUR
);
1548 listSize
= chunkSize
;
1549 block
= src
-> mSize
* src
-> mSkip
;
1550 skip
= src
-> mOffset
;
1553 while ((offset
< n
) && (listSize
> 8)) {
1554 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1555 (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & chunkSize
)))
1557 listSize
-= 8 + chunkSize
;
1558 if (fourCC
== FOURCC_DATA
) {
1559 count
= chunkSize
/ block
;
1561 fseek (fp
, (long) (skip
* block
), SEEK_CUR
);
1562 chunkSize
-= skip
* block
;
1565 if (count
> (n
- offset
))
1567 if (! ReadWaveData (fp
, src
, order
, count
, & hrir
[offset
]))
1569 chunkSize
-= count
* block
;
1571 lastSample
= hrir
[offset
- 1];
1576 } else if (fourCC
== FOURCC_SLNT
) {
1577 if (! ReadBin4 (fp
, src
-> mPath
, order
, 4, & count
))
1583 if (count
> (n
- offset
))
1585 for (i
= 0; i
< count
; i
++)
1586 hrir
[offset
+ i
] = lastSample
;
1594 fseek (fp
, (long) chunkSize
, SEEK_CUR
);
1597 fprintf (stderr
, "Error: Bad read from file '%s'.\n", src
-> mPath
);
1603 // Load a source HRIR from a RIFF/RIFX WAVE file.
1604 static int LoadWaveSource (FILE * fp
, SourceRefT
* src
, const uint hrirRate
, const uint n
, double * hrir
) {
1605 uint32 fourCC
, dummy
;
1608 if ((! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
)) ||
1609 (! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & dummy
)))
1611 if (fourCC
== FOURCC_RIFF
) {
1613 } else if (fourCC
== FOURCC_RIFX
) {
1616 fprintf (stderr
, "Error: No RIFF/RIFX chunk in file '%s'.\n", src
-> mPath
);
1619 if (! ReadBin4 (fp
, src
-> mPath
, BO_LITTLE
, 4, & fourCC
))
1621 if (fourCC
!= FOURCC_WAVE
) {
1622 fprintf (stderr
, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src
-> mPath
);
1625 if (! ReadWaveFormat (fp
, order
, hrirRate
, src
))
1627 if (! ReadWaveList (fp
, src
, order
, n
, hrir
))
1632 // Load a source HRIR from a binary file.
1633 static int LoadBinarySource (FILE * fp
, const SourceRefT
* src
, const ByteOrderT order
, const uint n
, double * hrir
) {
1636 fseek (fp
, (long) src
-> mOffset
, SEEK_SET
);
1637 for (i
= 0; i
< n
; i
++) {
1638 if (! ReadBinAsDouble (fp
, src
-> mPath
, order
, src
-> mType
, src
-> mSize
, src
-> mBits
, & hrir
[i
]))
1640 if (src
-> mSkip
> 0)
1641 fseek (fp
, (long) src
-> mSkip
, SEEK_CUR
);
1646 // Load a source HRIR from an ASCII text file containing a list of elements
1647 // separated by whitespace or common list operators (',', ';', ':', '|').
1648 static int LoadAsciiSource (FILE * fp
, const SourceRefT
* src
, const uint n
, double * hrir
) {
1653 TrSetup (fp
, NULL
, & tr
);
1654 for (i
= 0; i
< src
-> mOffset
; i
++) {
1655 if (! ReadAsciiAsDouble (& tr
, src
-> mPath
, src
-> mType
, (uint
) src
-> mBits
, & dummy
))
1658 for (i
= 0; i
< n
; i
++) {
1659 if (! ReadAsciiAsDouble (& tr
, src
-> mPath
, src
-> mType
, (uint
) src
-> mBits
, & hrir
[i
]))
1661 for (j
= 0; j
< src
-> mSkip
; j
++) {
1662 if (! ReadAsciiAsDouble (& tr
, src
-> mPath
, src
-> mType
, (uint
) src
-> mBits
, & dummy
))
1669 // Load a source HRIR from a supported file type.
1670 static int LoadSource (SourceRefT
* src
, const uint hrirRate
, const uint n
, double * hrir
) {
1674 if (src
-> mFormat
== SF_ASCII
)
1675 fp
= fopen (src
-> mPath
, "r");
1677 fp
= fopen (src
-> mPath
, "rb");
1679 fprintf (stderr
, "Error: Could not open source file '%s'.\n", src
-> mPath
);
1682 if (src
-> mFormat
== SF_WAVE
)
1683 result
= LoadWaveSource (fp
, src
, hrirRate
, n
, hrir
);
1684 else if (src
-> mFormat
== SF_BIN_LE
)
1685 result
= LoadBinarySource (fp
, src
, BO_LITTLE
, n
, hrir
);
1686 else if (src
-> mFormat
== SF_BIN_BE
)
1687 result
= LoadBinarySource (fp
, src
, BO_BIG
, n
, hrir
);
1689 result
= LoadAsciiSource (fp
, src
, n
, hrir
);
1694 // Calculate the onset time of an HRIR and average it with any existing
1695 // timing for its elevation and azimuth.
1696 static void AverageHrirOnset (const double * hrir
, const double f
, const uint ei
, const uint ai
, const HrirDataT
* hData
) {
1701 n
= hData
-> mIrPoints
;
1702 for (i
= 0; i
< n
; i
++)
1703 mag
= fmax (fabs (hrir
[i
]), mag
);
1705 for (i
= 0; i
< n
; i
++) {
1706 if (fabs (hrir
[i
]) >= mag
)
1709 j
= hData
-> mEvOffset
[ei
] + ai
;
1710 hData
-> mHrtds
[j
] = Lerp (hData
-> mHrtds
[j
], ((double) i
) / hData
-> mIrRate
, f
);
1713 // Calculate the magnitude response of an HRIR and average it with any
1714 // existing responses for its elevation and azimuth.
1715 static void AverageHrirMagnitude (const double * hrir
, const double f
, const uint ei
, const uint ai
, const HrirDataT
* hData
) {
1716 double * re
= NULL
, * im
= NULL
;
1719 n
= hData
-> mFftSize
;
1720 re
= CreateArray (n
);
1721 im
= CreateArray (n
);
1722 for (i
= 0; i
< hData
-> mIrPoints
; i
++) {
1726 for (; i
< n
; i
++) {
1730 FftForward (n
, re
, im
, re
, im
);
1731 MagnitudeResponse (n
, re
, im
, re
);
1733 j
= (hData
-> mEvOffset
[ei
] + ai
) * hData
-> mIrSize
;
1734 for (i
= 0; i
< m
; i
++)
1735 hData
-> mHrirs
[j
+ i
] = Lerp (hData
-> mHrirs
[j
+ i
], re
[i
], f
);
1740 /* Calculate the contribution of each HRIR to the diffuse-field average based
1741 * on the area of its surface patch. All patches are centered at the HRIR
1742 * coordinates on the unit sphere and are measured by solid angle.
1744 static void CalculateDfWeights (const HrirDataT
* hData
, double * weights
) {
1746 double evs
, sum
, ev
, up_ev
, down_ev
, solidAngle
;
1748 evs
= 90.0 / (hData
-> mEvCount
- 1);
1750 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++) {
1751 // For each elevation, calculate the upper and lower limits of the
1753 ev
= -90.0 + (ei
* 2.0 * evs
);
1754 if (ei
< (hData
-> mEvCount
- 1))
1755 up_ev
= (ev
+ evs
) * M_PI
/ 180.0;
1759 down_ev
= (ev
- evs
) * M_PI
/ 180.0;
1761 down_ev
= -M_PI
/ 2.0;
1762 // Calculate the area of the patch band.
1763 solidAngle
= 2.0 * M_PI
* (sin (up_ev
) - sin (down_ev
));
1764 // Each weight is the area of one patch.
1765 weights
[ei
] = solidAngle
/ hData
-> mAzCount
[ei
];
1766 // Sum the total surface area covered by the HRIRs.
1769 // Normalize the weights given the total surface coverage.
1770 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++)
1771 weights
[ei
] /= sum
;
1774 /* Calculate the diffuse-field average from the given magnitude responses of
1775 * the HRIR set. Weighting can be applied to compensate for the varying
1776 * surface area covered by each HRIR. The final average can then be limited
1777 * by the specified magnitude range (in positive dB; 0.0 to skip).
1779 static void CalculateDiffuseFieldAverage (const HrirDataT
* hData
, const int weighted
, const double limit
, double * dfa
) {
1780 double * weights
= NULL
;
1781 uint ei
, ai
, count
, step
, start
, end
, m
, j
, i
;
1784 weights
= CreateArray (hData
-> mEvCount
);
1786 // Use coverage weighting to calculate the average.
1787 CalculateDfWeights (hData
, weights
);
1789 // If coverage weighting is not used, the weights still need to be
1790 // averaged by the number of HRIRs.
1792 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++)
1793 count
+= hData
-> mAzCount
[ei
];
1794 for (ei
= hData
-> mEvStart
; ei
< hData
-> mEvCount
; ei
++)
1795 weights
[ei
] = 1.0 / count
;
1797 ei
= hData
-> mEvStart
;
1799 step
= hData
-> mIrSize
;
1800 start
= hData
-> mEvOffset
[ei
] * step
;
1801 end
= hData
-> mIrCount
* step
;
1802 m
= 1 + (hData
-> mFftSize
/ 2);
1803 for (i
= 0; i
< m
; i
++)
1805 for (j
= start
; j
< end
; j
+= step
) {
1806 // Get the weight for this HRIR's contribution.
1807 weight
= weights
[ei
];
1808 // Add this HRIR's weighted power average to the total.
1809 for (i
= 0; i
< m
; i
++)
1810 dfa
[i
] += weight
* hData
-> mHrirs
[j
+ i
] * hData
-> mHrirs
[j
+ i
];
1811 // Determine the next weight to use.
1813 if (ai
>= hData
-> mAzCount
[ei
]) {
1818 // Finish the average calculation and keep it from being too small.
1819 for (i
= 0; i
< m
; i
++)
1820 dfa
[i
] = fmax (sqrt (dfa
[i
]), EPSILON
);
1821 // Apply a limit to the magnitude range of the diffuse-field average if
1824 LimitMagnitudeResponse (hData
-> mFftSize
, limit
, dfa
, dfa
);
1825 DestroyArray (weights
);
1828 // Perform diffuse-field equalization on the magnitude responses of the HRIR
1829 // set using the given average response.
1830 static void DiffuseFieldEqualize (const double * dfa
, const HrirDataT
* hData
) {
1831 uint step
, start
, end
, m
, j
, i
;
1833 step
= hData
-> mIrSize
;
1834 start
= hData
-> mEvOffset
[hData
-> mEvStart
] * step
;
1835 end
= hData
-> mIrCount
* step
;
1836 m
= 1 + (hData
-> mFftSize
/ 2);
1837 for (j
= start
; j
< end
; j
+= step
) {
1838 for (i
= 0; i
< m
; i
++)
1839 hData
-> mHrirs
[j
+ i
] /= dfa
[i
];
1843 // Perform minimum-phase reconstruction using the magnitude responses of the
1845 static void ReconstructHrirs (const HrirDataT
* hData
) {
1846 double * re
= NULL
, * im
= NULL
;
1847 uint step
, start
, end
, n
, j
, i
;
1849 step
= hData
-> mIrSize
;
1850 start
= hData
-> mEvOffset
[hData
-> mEvStart
] * step
;
1851 end
= hData
-> mIrCount
* step
;
1852 n
= hData
-> mFftSize
;
1853 re
= CreateArray (n
);
1854 im
= CreateArray (n
);
1855 for (j
= start
; j
< end
; j
+= step
) {
1856 MinimumPhase (n
, & hData
-> mHrirs
[j
], re
, im
);
1857 FftInverse (n
, re
, im
, re
, im
);
1858 for (i
= 0; i
< hData
-> mIrPoints
; i
++)
1859 hData
-> mHrirs
[j
+ i
] = re
[i
];
1865 // Resamples the HRIRs for use at the given sampling rate.
1866 static void ResampleHrirs (const uint rate
, HrirDataT
* hData
) {
1868 uint n
, step
, start
, end
, j
;
1870 ResamplerSetup (& rs
, hData
-> mIrRate
, rate
);
1871 n
= hData
-> mIrPoints
;
1872 step
= hData
-> mIrSize
;
1873 start
= hData
-> mEvOffset
[hData
-> mEvStart
] * step
;
1874 end
= hData
-> mIrCount
* step
;
1875 for (j
= start
; j
< end
; j
+= step
)
1876 ResamplerRun (& rs
, n
, & hData
-> mHrirs
[j
], n
, & hData
-> mHrirs
[j
]);
1877 ResamplerClear (& rs
);
1878 hData
-> mIrRate
= rate
;
1881 /* Given an elevation index and an azimuth, calculate the indices of the two
1882 * HRIRs that bound the coordinate along with a factor for calculating the
1883 * continous HRIR using interpolation.
1885 static void CalcAzIndices (const HrirDataT
* hData
, const uint ei
, const double az
, uint
* j0
, uint
* j1
, double * jf
) {
1889 af
= ((2.0 * M_PI
) + az
) * hData
-> mAzCount
[ei
] / (2.0 * M_PI
);
1890 ai
= ((uint
) af
) % hData
-> mAzCount
[ei
];
1892 (* j0
) = hData
-> mEvOffset
[ei
] + ai
;
1893 (* j1
) = hData
-> mEvOffset
[ei
] + ((ai
+ 1) % hData
-> mAzCount
[ei
]);
1897 // Synthesize any missing onset timings at the bottom elevations. This just
1898 // blends between slightly exaggerated known onsets. Not an accurate model.
1899 static void SynthesizeOnsets (HrirDataT
* hData
) {
1900 uint oi
, e
, a
, j0
, j1
;
1903 oi
= hData
-> mEvStart
;
1905 for (a
= 0; a
< hData
-> mAzCount
[oi
]; a
++)
1906 t
+= hData
-> mHrtds
[hData
-> mEvOffset
[oi
] + a
];
1907 hData
-> mHrtds
[0] = 1.32e-4 + (t
/ hData
-> mAzCount
[oi
]);
1908 for (e
= 1; e
< hData
-> mEvStart
; e
++) {
1909 of
= ((double) e
) / hData
-> mEvStart
;
1910 for (a
= 0; a
< hData
-> mAzCount
[e
]; a
++) {
1911 CalcAzIndices (hData
, oi
, a
* 2.0 * M_PI
/ hData
-> mAzCount
[e
], & j0
, & j1
, & jf
);
1912 hData
-> mHrtds
[hData
-> mEvOffset
[e
] + a
] = Lerp (hData
-> mHrtds
[0], Lerp (hData
-> mHrtds
[j0
], hData
-> mHrtds
[j1
], jf
), of
);
1917 /* Attempt to synthesize any missing HRIRs at the bottom elevations. Right
1918 * now this just blends the lowest elevation HRIRs together and applies some
1919 * attenuation and high frequency damping. It is a simple, if inaccurate
1922 static void SynthesizeHrirs (HrirDataT
* hData
) {
1923 uint oi
, a
, e
, step
, n
, i
, j
;
1927 double lp
[4], s0
, s1
;
1929 if (hData
-> mEvStart
<= 0)
1931 step
= hData
-> mIrSize
;
1932 oi
= hData
-> mEvStart
;
1933 n
= hData
-> mIrPoints
;
1934 for (i
= 0; i
< n
; i
++)
1935 hData
-> mHrirs
[i
] = 0.0;
1936 for (a
= 0; a
< hData
-> mAzCount
[oi
]; a
++) {
1937 j
= (hData
-> mEvOffset
[oi
] + a
) * step
;
1938 for (i
= 0; i
< n
; i
++)
1939 hData
-> mHrirs
[i
] += hData
-> mHrirs
[j
+ i
] / hData
-> mAzCount
[oi
];
1941 for (e
= 1; e
< hData
-> mEvStart
; e
++) {
1942 of
= ((double) e
) / hData
-> mEvStart
;
1943 b
= (1.0 - of
) * (3.5e-6 * hData
-> mIrRate
);
1944 for (a
= 0; a
< hData
-> mAzCount
[e
]; a
++) {
1945 j
= (hData
-> mEvOffset
[e
] + a
) * step
;
1946 CalcAzIndices (hData
, oi
, a
* 2.0 * M_PI
/ hData
-> mAzCount
[e
], & j0
, & j1
, & jf
);
1953 for (i
= 0; i
< n
; i
++) {
1954 s0
= hData
-> mHrirs
[i
];
1955 s1
= Lerp (hData
-> mHrirs
[j0
+ i
], hData
-> mHrirs
[j1
+ i
], jf
);
1956 s0
= Lerp (s0
, s1
, of
);
1957 lp
[0] = Lerp (s0
, lp
[0], b
);
1958 lp
[1] = Lerp (lp
[0], lp
[1], b
);
1959 lp
[2] = Lerp (lp
[1], lp
[2], b
);
1960 lp
[3] = Lerp (lp
[2], lp
[3], b
);
1961 hData
-> mHrirs
[j
+ i
] = lp
[3];
1965 b
= 3.5e-6 * hData
-> mIrRate
;
1970 for (i
= 0; i
< n
; i
++) {
1971 s0
= hData
-> mHrirs
[i
];
1972 lp
[0] = Lerp (s0
, lp
[0], b
);
1973 lp
[1] = Lerp (lp
[0], lp
[1], b
);
1974 lp
[2] = Lerp (lp
[1], lp
[2], b
);
1975 lp
[3] = Lerp (lp
[2], lp
[3], b
);
1976 hData
-> mHrirs
[i
] = lp
[3];
1978 hData
-> mEvStart
= 0;
1981 // The following routines assume a full set of HRIRs for all elevations.
1983 // Normalize the HRIR set and slightly attenuate the result.
1984 static void NormalizeHrirs (const HrirDataT
* hData
) {
1985 uint step
, end
, n
, j
, i
;
1988 step
= hData
-> mIrSize
;
1989 end
= hData
-> mIrCount
* step
;
1990 n
= hData
-> mIrPoints
;
1992 for (j
= 0; j
< end
; j
+= step
) {
1993 for (i
= 0; i
< n
; i
++)
1994 maxLevel
= fmax (fabs (hData
-> mHrirs
[j
+ i
]), maxLevel
);
1996 maxLevel
= 1.01 * maxLevel
;
1997 for (j
= 0; j
< end
; j
+= step
) {
1998 for (i
= 0; i
< n
; i
++)
1999 hData
-> mHrirs
[j
+ i
] /= maxLevel
;
2003 // Calculate the left-ear time delay using a spherical head model.
2004 static double CalcLTD (const double ev
, const double az
, const double rad
, const double dist
) {
2005 double azp
, dlp
, l
, al
;
2007 azp
= asin (cos (ev
) * sin (az
));
2008 dlp
= sqrt ((dist
* dist
) + (rad
* rad
) + (2.0 * dist
* rad
* sin (azp
)));
2009 l
= sqrt ((dist
* dist
) - (rad
* rad
));
2010 al
= (0.5 * M_PI
) + azp
;
2012 dlp
= l
+ (rad
* (al
- acos (rad
/ dist
)));
2013 return (dlp
/ 343.3);
2016 // Calculate the effective head-related time delays for each minimum-phase
2018 static void CalculateHrtds (const HeadModelT model
, const double radius
, HrirDataT
* hData
) {
2019 double minHrtd
, maxHrtd
;
2025 for (e
= 0; e
< hData
-> mEvCount
; e
++) {
2026 for (a
= 0; a
< hData
-> mAzCount
[e
]; a
++) {
2027 j
= hData
-> mEvOffset
[e
] + a
;
2028 if (model
== HM_DATASET
) {
2029 t
= hData
-> mHrtds
[j
] * radius
/ hData
-> mRadius
;
2031 t
= CalcLTD ((-90.0 + (e
* 180.0 / (hData
-> mEvCount
- 1))) * M_PI
/ 180.0,
2032 (a
* 360.0 / hData
-> mAzCount
[e
]) * M_PI
/ 180.0,
2033 radius
, hData
-> mDistance
);
2035 hData
-> mHrtds
[j
] = t
;
2036 maxHrtd
= fmax (t
, maxHrtd
);
2037 minHrtd
= fmin (t
, minHrtd
);
2041 for (j
= 0; j
< hData
-> mIrCount
; j
++)
2042 hData
-> mHrtds
[j
] -= minHrtd
;
2043 hData
-> mMaxHrtd
= maxHrtd
;
2046 // Store the OpenAL Soft HRTF data set.
2047 static int StoreMhr (const HrirDataT
* hData
, const char * filename
) {
2049 uint e
, step
, end
, n
, j
, i
;
2052 if ((fp
= fopen (filename
, "wb")) == NULL
) {
2053 fprintf (stderr
, "Error: Could not open MHR file '%s'.\n", filename
);
2056 if (! WriteAscii (MHR_FORMAT
, fp
, filename
))
2058 if (! WriteBin4 (BO_LITTLE
, 4, (uint32
) hData
-> mIrRate
, fp
, filename
))
2060 if (! WriteBin4 (BO_LITTLE
, 1, (uint32
) hData
-> mIrPoints
, fp
, filename
))
2062 if (! WriteBin4 (BO_LITTLE
, 1, (uint32
) hData
-> mEvCount
, fp
, filename
))
2064 for (e
= 0; e
< hData
-> mEvCount
; e
++) {
2065 if (! WriteBin4 (BO_LITTLE
, 1, (uint32
) hData
-> mAzCount
[e
], fp
, filename
))
2068 step
= hData
-> mIrSize
;
2069 end
= hData
-> mIrCount
* step
;
2070 n
= hData
-> mIrPoints
;
2072 for (j
= 0; j
< end
; j
+= step
) {
2074 for (i
= 0; i
< n
; i
++) {
2075 v
= HpTpdfDither (32767.0 * hData
-> mHrirs
[j
+ i
], & hpHist
);
2076 if (! WriteBin4 (BO_LITTLE
, 2, (uint32
) v
, fp
, filename
))
2080 for (j
= 0; j
< hData
-> mIrCount
; j
++) {
2081 v
= (int) fmin (round (hData
-> mIrRate
* hData
-> mHrtds
[j
]), MAX_HRTD
);
2082 if (! WriteBin4 (BO_LITTLE
, 1, (uint32
) v
, fp
, filename
))
2089 // Store the OpenAL Soft built-in table.
2090 static int StoreTable (const HrirDataT
* hData
, const char * filename
) {
2092 uint step
, end
, n
, j
, i
;
2094 char text
[128 + 1];
2096 if ((fp
= fopen (filename
, "wb")) == NULL
) {
2097 fprintf (stderr
, "Error: Could not open table file '%s'.\n", filename
);
2100 snprintf (text
, 128, "/* Elevation metrics */\n"
2101 "static const ALubyte defaultAzCount[%u] = { ", hData
-> mEvCount
);
2102 if (! WriteAscii (text
, fp
, filename
))
2104 for (i
= 0; i
< hData
-> mEvCount
; i
++) {
2105 snprintf (text
, 128, "%u, ", hData
-> mAzCount
[i
]);
2106 if (! WriteAscii (text
, fp
, filename
))
2109 snprintf (text
, 128, "};\n"
2110 "static const ALushort defaultEvOffset[%u] = { ", hData
-> mEvCount
);
2111 if (! WriteAscii (text
, fp
, filename
))
2113 for (i
= 0; i
< hData
-> mEvCount
; i
++) {
2114 snprintf (text
, 128, "%u, ", hData
-> mEvOffset
[i
]);
2115 if (! WriteAscii (text
, fp
, filename
))
2118 step
= hData
-> mIrSize
;
2119 end
= hData
-> mIrCount
* step
;
2120 n
= hData
-> mIrPoints
;
2121 snprintf (text
, 128, "};\n\n"
2122 "/* HRIR Coefficients */\n"
2123 "static const ALshort defaultCoeffs[%u] =\n{\n", hData
-> mIrCount
* n
);
2124 if (! WriteAscii (text
, fp
, filename
))
2127 for (j
= 0; j
< end
; j
+= step
) {
2128 if (! WriteAscii (" ", fp
, filename
))
2131 for (i
= 0; i
< n
; i
++) {
2132 v
= HpTpdfDither (32767.0 * hData
-> mHrirs
[j
+ i
], & hpHist
);
2133 snprintf (text
, 128, " %+d,", v
);
2134 if (! WriteAscii (text
, fp
, filename
))
2137 if (! WriteAscii ("\n", fp
, filename
))
2140 snprintf (text
, 128, "};\n\n"
2141 "/* HRIR Delays */\n"
2142 "static const ALubyte defaultDelays[%u] =\n{\n"
2143 " ", hData
-> mIrCount
);
2144 if (! WriteAscii (text
, fp
, filename
))
2146 for (j
= 0; j
< hData
-> mIrCount
; j
++) {
2147 v
= (int) fmin (round (hData
-> mIrRate
* hData
-> mHrtds
[j
]), MAX_HRTD
);
2148 snprintf (text
, 128, " %d,", v
);
2149 if (! WriteAscii (text
, fp
, filename
))
2152 if (! WriteAscii ("\n};\n\n"
2153 "/* Default HRTF Definition */\n", fp
, filename
))
2155 snprintf (text
, 128, "static const struct Hrtf DefaultHrtf = {\n"
2156 " %u, %u, %u, defaultAzCount, defaultEvOffset,\n",
2157 hData
-> mIrRate
, hData
-> mIrPoints
, hData
-> mEvCount
);
2158 if (! WriteAscii (text
, fp
, filename
))
2160 if (! WriteAscii (" defaultCoeffs, defaultDelays, NULL\n"
2161 "};\n", fp
, filename
))
2167 // Process the data set definition to read and validate the data set metrics.
2168 static int ProcessMetrics (TokenReaderT
* tr
, const uint fftSize
, const uint truncSize
, HrirDataT
* hData
) {
2169 char ident
[MAX_IDENT_LEN
+ 1];
2174 int hasRate
= 0, hasPoints
= 0, hasAzimuths
= 0;
2175 int hasRadius
= 0, hasDistance
= 0;
2177 while (! (hasRate
&& hasPoints
&& hasAzimuths
&& hasRadius
&& hasDistance
)) {
2178 TrIndication (tr
, & line
, & col
);
2179 if (! TrReadIdent (tr
, MAX_IDENT_LEN
, ident
))
2181 if (strcasecmp (ident
, "rate") == 0) {
2183 TrErrorAt (tr
, line
, col
, "Redefinition of 'rate'.\n");
2186 if (! TrReadOperator (tr
, "="))
2188 if (! TrReadInt (tr
, MIN_RATE
, MAX_RATE
, & intVal
))
2190 hData
-> mIrRate
= (uint
) intVal
;
2192 } else if (strcasecmp (ident
, "points") == 0) {
2194 TrErrorAt (tr
, line
, col
, "Redefinition of 'points'.\n");
2197 if (! TrReadOperator (tr
, "="))
2199 TrIndication (tr
, & line
, & col
);
2200 if (! TrReadInt (tr
, MIN_POINTS
, MAX_POINTS
, & intVal
))
2202 points
= (uint
) intVal
;
2203 if ((fftSize
> 0) && (points
> fftSize
)) {
2204 TrErrorAt (tr
, line
, col
, "Value exceeds the overriden FFT size.\n");
2207 if (points
< truncSize
) {
2208 TrErrorAt (tr
, line
, col
, "Value is below the truncation size.\n");
2211 hData
-> mIrPoints
= points
;
2212 hData
-> mFftSize
= fftSize
;
2215 while (points
< (4 * hData
-> mIrPoints
))
2217 hData
-> mFftSize
= points
;
2218 hData
-> mIrSize
= 1 + (points
/ 2);
2220 hData
-> mFftSize
= fftSize
;
2221 hData
-> mIrSize
= 1 + (fftSize
/ 2);
2222 if (points
> hData
-> mIrSize
)
2223 hData
-> mIrSize
= points
;
2226 } else if (strcasecmp (ident
, "azimuths") == 0) {
2228 TrErrorAt (tr
, line
, col
, "Redefinition of 'azimuths'.\n");
2231 if (! TrReadOperator (tr
, "="))
2233 hData
-> mIrCount
= 0;
2234 hData
-> mEvCount
= 0;
2235 hData
-> mEvOffset
[0] = 0;
2237 if (! TrReadInt (tr
, MIN_AZ_COUNT
, MAX_AZ_COUNT
, & intVal
))
2239 hData
-> mAzCount
[hData
-> mEvCount
] = (uint
) intVal
;
2240 hData
-> mIrCount
+= (uint
) intVal
;
2241 hData
-> mEvCount
++;
2242 if (! TrIsOperator (tr
, ","))
2244 if (hData
-> mEvCount
>= MAX_EV_COUNT
) {
2245 TrError (tr
, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT
);
2248 hData
-> mEvOffset
[hData
-> mEvCount
] = hData
-> mEvOffset
[hData
-> mEvCount
- 1] + ((uint
) intVal
);
2249 TrReadOperator (tr
, ",");
2251 if (hData
-> mEvCount
< MIN_EV_COUNT
) {
2252 TrErrorAt (tr
, line
, col
, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT
);
2256 } else if (strcasecmp (ident
, "radius") == 0) {
2258 TrErrorAt (tr
, line
, col
, "Redefinition of 'radius'.\n");
2261 if (! TrReadOperator (tr
, "="))
2263 if (! TrReadFloat (tr
, MIN_RADIUS
, MAX_RADIUS
, & fpVal
))
2265 hData
-> mRadius
= fpVal
;
2267 } else if (strcasecmp (ident
, "distance") == 0) {
2269 TrErrorAt (tr
, line
, col
, "Redefinition of 'distance'.\n");
2272 if (! TrReadOperator (tr
, "="))
2274 if (! TrReadFloat (tr
, MIN_DISTANCE
, MAX_DISTANCE
, & fpVal
))
2276 hData
-> mDistance
= fpVal
;
2279 TrErrorAt (tr
, line
, col
, "Expected a metric name.\n");
2282 TrSkipWhitespace (tr
);
2287 // Parse an index pair from the data set definition.
2288 static int ReadIndexPair (TokenReaderT
* tr
, const HrirDataT
* hData
, uint
* ei
, uint
* ai
) {
2291 if (! TrReadInt (tr
, 0, (int) hData
-> mEvCount
, & intVal
))
2293 (* ei
) = (uint
) intVal
;
2294 if (! TrReadOperator (tr
, ","))
2296 if (! TrReadInt (tr
, 0, (int) hData
-> mAzCount
[(* ei
)], & intVal
))
2298 (* ai
) = (uint
) intVal
;
2302 // Match the source format from a given identifier.
2303 static SourceFormatT
MatchSourceFormat (const char * ident
) {
2304 if (strcasecmp (ident
, "wave") == 0)
2306 else if (strcasecmp (ident
, "bin_le") == 0)
2308 else if (strcasecmp (ident
, "bin_be") == 0)
2310 else if (strcasecmp (ident
, "ascii") == 0)
2315 // Match the source element type from a given identifier.
2316 static ElementTypeT
MatchElementType (const char * ident
) {
2317 if (strcasecmp (ident
, "int") == 0)
2319 else if (strcasecmp (ident
, "fp") == 0)
2324 // Parse and validate a source reference from the data set definition.
2325 static int ReadSourceRef (TokenReaderT
* tr
, SourceRefT
* src
) {
2327 char ident
[MAX_IDENT_LEN
+ 1];
2330 TrIndication (tr
, & line
, & col
);
2331 if (! TrReadIdent (tr
, MAX_IDENT_LEN
, ident
))
2333 src
-> mFormat
= MatchSourceFormat (ident
);
2334 if (src
-> mFormat
== SF_NONE
) {
2335 TrErrorAt (tr
, line
, col
, "Expected a source format.\n");
2338 if (! TrReadOperator (tr
, "("))
2340 if (src
-> mFormat
== SF_WAVE
) {
2341 if (! TrReadInt (tr
, 0, MAX_WAVE_CHANNELS
, & intVal
))
2343 src
-> mType
= ET_NONE
;
2346 src
-> mChannel
= (uint
) intVal
;
2349 TrIndication (tr
, & line
, & col
);
2350 if (! TrReadIdent (tr
, MAX_IDENT_LEN
, ident
))
2352 src
-> mType
= MatchElementType (ident
);
2353 if (src
-> mType
== ET_NONE
) {
2354 TrErrorAt (tr
, line
, col
, "Expected a source element type.\n");
2357 if ((src
-> mFormat
== SF_BIN_LE
) || (src
-> mFormat
== SF_BIN_BE
)) {
2358 if (! TrReadOperator (tr
, ","))
2360 if (src
-> mType
== ET_INT
) {
2361 if (! TrReadInt (tr
, MIN_BIN_SIZE
, MAX_BIN_SIZE
, & intVal
))
2363 src
-> mSize
= (uint
) intVal
;
2364 if (TrIsOperator (tr
, ",")) {
2365 TrReadOperator (tr
, ",");
2366 TrIndication (tr
, & line
, & col
);
2367 if (! TrReadInt (tr
, -2147483647 - 1, 2147483647, & intVal
))
2369 if ((abs (intVal
) < MIN_BIN_BITS
) || (((uint
) abs (intVal
)) > (8 * src
-> mSize
))) {
2370 TrErrorAt (tr
, line
, col
, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS
, 8 * src
-> mSize
);
2373 src
-> mBits
= intVal
;
2375 src
-> mBits
= (int) (8 * src
-> mSize
);
2378 TrIndication (tr
, & line
, & col
);
2379 if (! TrReadInt (tr
, -2147483647 - 1, 2147483647, & intVal
))
2381 if ((intVal
!= 4) && (intVal
!= 8)) {
2382 TrErrorAt (tr
, line
, col
, "Expected a value of 4 or 8.\n");
2385 src
-> mSize
= (uint
) intVal
;
2388 } else if ((src
-> mFormat
== SF_ASCII
) && (src
-> mType
== ET_INT
)) {
2389 if (! TrReadOperator (tr
, ","))
2391 if (! TrReadInt (tr
, MIN_ASCII_BITS
, MAX_ASCII_BITS
, & intVal
))
2394 src
-> mBits
= intVal
;
2399 if (TrIsOperator (tr
, ";")) {
2400 TrReadOperator (tr
, ";");
2401 if (! TrReadInt (tr
, 0, 0x7FFFFFFF, & intVal
))
2403 src
-> mSkip
= (uint
) intVal
;
2408 if (! TrReadOperator (tr
, ")"))
2410 if (TrIsOperator (tr
, "@")) {
2411 TrReadOperator (tr
, "@");
2412 if (! TrReadInt (tr
, 0, 0x7FFFFFFF, & intVal
))
2414 src
-> mOffset
= (uint
) intVal
;
2418 if (! TrReadOperator (tr
, ":"))
2420 if (! TrReadString (tr
, MAX_PATH_LEN
, src
-> mPath
))
2425 // Process the list of sources in the data set definition.
2426 static int ProcessSources (const HeadModelT model
, TokenReaderT
* tr
, HrirDataT
* hData
) {
2427 uint
* setCount
= NULL
, * setFlag
= NULL
;
2428 double * hrir
= NULL
;
2429 uint line
, col
, ei
, ai
;
2433 setCount
= (uint
*) calloc (hData
-> mEvCount
, sizeof (uint
));
2434 setFlag
= (uint
*) calloc (hData
-> mIrCount
, sizeof (uint
));
2435 hrir
= CreateArray (hData
-> mIrPoints
);
2436 while (TrIsOperator (tr
, "[")) {
2437 TrIndication (tr
, & line
, & col
);
2438 TrReadOperator (tr
, "[");
2439 if (ReadIndexPair (tr
, hData
, & ei
, & ai
)) {
2440 if (TrReadOperator (tr
, "]")) {
2441 if (! setFlag
[hData
-> mEvOffset
[ei
] + ai
]) {
2442 if (TrReadOperator (tr
, "=")) {
2445 if (ReadSourceRef (tr
, & src
)) {
2446 if (LoadSource (& src
, hData
-> mIrRate
, hData
-> mIrPoints
, hrir
)) {
2447 if (model
== HM_DATASET
)
2448 AverageHrirOnset (hrir
, 1.0 / factor
, ei
, ai
, hData
);
2449 AverageHrirMagnitude (hrir
, 1.0 / factor
, ei
, ai
, hData
);
2451 if (! TrIsOperator (tr
, "+"))
2453 TrReadOperator (tr
, "+");
2457 DestroyArray (hrir
);
2462 setFlag
[hData
-> mEvOffset
[ei
] + ai
] = 1;
2467 TrErrorAt (tr
, line
, col
, "Redefinition of source.\n");
2471 DestroyArray (hrir
);
2477 while ((ei
< hData
-> mEvCount
) && (setCount
[ei
] < 1))
2479 if (ei
< hData
-> mEvCount
) {
2480 hData
-> mEvStart
= ei
;
2481 while ((ei
< hData
-> mEvCount
) && (setCount
[ei
] == hData
-> mAzCount
[ei
]))
2483 if (ei
>= hData
-> mEvCount
) {
2484 if (! TrLoad (tr
)) {
2485 DestroyArray (hrir
);
2490 TrError (tr
, "Errant data at end of source list.\n");
2493 TrError (tr
, "Missing sources for elevation index %d.\n", ei
);
2496 TrError (tr
, "Missing source references.\n");
2498 DestroyArray (hrir
);
2504 /* Parse the data set definition and process the source data, storing the
2505 * resulting data set as desired. If the input name is NULL it will read
2506 * from standard input.
2508 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 HeadModelT model
, const double radius
, const OutputFormatT outFormat
, const char * outName
) {
2512 double * dfa
= NULL
;
2513 char rateStr
[8 + 1], expName
[MAX_PATH_LEN
];
2515 hData
. mIrRate
= 0;
2516 hData
. mIrPoints
= 0;
2517 hData
. mFftSize
= 0;
2518 hData
. mIrSize
= 0;
2519 hData
. mIrCount
= 0;
2520 hData
. mEvCount
= 0;
2521 hData
. mRadius
= 0;
2522 hData
. mDistance
= 0;
2523 fprintf (stdout
, "Reading HRIR definition...\n");
2524 if (inName
!= NULL
) {
2525 fp
= fopen (inName
, "r");
2527 fprintf (stderr
, "Error: Could not open definition file '%s'\n", inName
);
2530 TrSetup (fp
, inName
, & tr
);
2533 TrSetup (fp
, "<stdin>", & tr
);
2535 if (! ProcessMetrics (& tr
, fftSize
, truncSize
, & hData
)) {
2540 hData
. mHrirs
= CreateArray (hData
. mIrCount
* hData
. mIrSize
);
2541 hData
. mHrtds
= CreateArray (hData
. mIrCount
);
2542 if (! ProcessSources (model
, & tr
, & hData
)) {
2543 DestroyArray (hData
. mHrtds
);
2544 DestroyArray (hData
. mHrirs
);
2552 dfa
= CreateArray (1 + (hData
. mFftSize
/ 2));
2553 fprintf (stdout
, "Calculating diffuse-field average...\n");
2554 CalculateDiffuseFieldAverage (& hData
, surface
, limit
, dfa
);
2555 fprintf (stdout
, "Performing diffuse-field equalization...\n");
2556 DiffuseFieldEqualize (dfa
, & hData
);
2559 fprintf (stdout
, "Performing minimum phase reconstruction...\n");
2560 ReconstructHrirs (& hData
);
2561 if ((outRate
!= 0) && (outRate
!= hData
. mIrRate
)) {
2562 fprintf (stdout
, "Resampling HRIRs...\n");
2563 ResampleHrirs (outRate
, & hData
);
2565 fprintf (stdout
, "Truncating minimum-phase HRIRs...\n");
2566 hData
. mIrPoints
= truncSize
;
2567 fprintf (stdout
, "Synthesizing missing elevations...\n");
2568 if (model
== HM_DATASET
)
2569 SynthesizeOnsets (& hData
);
2570 SynthesizeHrirs (& hData
);
2571 fprintf (stdout
, "Normalizing final HRIRs...\n");
2572 NormalizeHrirs (& hData
);
2573 fprintf (stdout
, "Calculating impulse delays...\n");
2574 CalculateHrtds (model
, (radius
> DEFAULT_CUSTOM_RADIUS
) ? radius
: hData
. mRadius
, & hData
);
2575 snprintf (rateStr
, 8, "%u", hData
. mIrRate
);
2576 StrSubst (outName
, "%r", rateStr
, MAX_PATH_LEN
, expName
);
2577 switch (outFormat
) {
2579 fprintf (stdout
, "Creating MHR data set file...\n");
2580 if (! StoreMhr (& hData
, expName
))
2584 fprintf (stderr
, "Creating OpenAL Soft table file...\n");
2585 if (! StoreTable (& hData
, expName
))
2591 DestroyArray (hData
. mHrtds
);
2592 DestroyArray (hData
. mHrirs
);
2596 // Standard command line dispatch.
2597 int main (const int argc
, const char * argv
[]) {
2598 const char * inName
= NULL
, * outName
= NULL
;
2599 OutputFormatT outFormat
;
2601 uint outRate
, fftSize
;
2602 int equalize
, surface
;
2610 fprintf (stderr
, "Error: No command specified. See '%s -h' for help.\n", argv
[0]);
2613 if ((strcmp (argv
[1], "--help") == 0) || (strcmp (argv
[1], "-h") == 0)) {
2614 fprintf (stdout
, "HRTF Processing and Composition Utility\n\n");
2615 fprintf (stdout
, "Usage: %s <command> [<option>...]\n\n", argv
[0]);
2616 fprintf (stdout
, "Commands:\n");
2617 fprintf (stdout
, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n");
2618 fprintf (stdout
, " Defaults output to: ./oalsoft_hrtf_%%r.mhr\n");
2619 fprintf (stdout
, " -t, --make-tab Makes the built-in table used when compiling OpenAL Soft.\n");
2620 fprintf (stdout
, " Defaults output to: ./hrtf_tables.inc\n");
2621 fprintf (stdout
, " -h, --help Displays this help information.\n\n");
2622 fprintf (stdout
, "Options:\n");
2623 fprintf (stdout
, " -r=<rate> Change the data set sample rate to the specified value and\n");
2624 fprintf (stdout
, " resample the HRIRs accordingly.\n");
2625 fprintf (stdout
, " -f=<points> Override the FFT window size (defaults to the first power-\n");
2626 fprintf (stdout
, " of-two that fits four times the number of HRIR points).\n");
2627 fprintf (stdout
, " -e={on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE
? "on" : "off"));
2628 fprintf (stdout
, " -s={on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE
? "on" : "off"));
2629 fprintf (stdout
, " -l={<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
2630 fprintf (stdout
, " average (default: %.2f).\n", DEFAULT_LIMIT
);
2631 fprintf (stdout
, " -w=<points> Specify the size of the truncation window that's applied\n");
2632 fprintf (stdout
, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE
);
2633 fprintf (stdout
, " -d={dataset| Specify the model used for calculating the head-delay timing\n");
2634 fprintf (stdout
, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL
== HM_DATASET
) ? "dataset" : "sphere"));
2635 fprintf (stdout
, " -c=<size> Use a customized head radius measured ear-to-ear in meters.\n");
2636 fprintf (stdout
, " -i=<filename> Specify an HRIR definition file to use (defaults to stdin).\n");
2637 fprintf (stdout
, " -o=<filename> Specify an output file. Overrides command-selected default.\n");
2638 fprintf (stdout
, " Use of '%%r' will be substituted with the data set sample rate.\n");
2641 if ((strcmp (argv
[1], "--make-mhr") == 0) || (strcmp (argv
[1], "-m") == 0)) {
2645 outName
= "./oalsoft_hrtf_%r.mhr";
2647 } else if ((strcmp (argv
[1], "--make-tab") == 0) || (strcmp (argv
[1], "-t") == 0)) {
2651 outName
= "./hrtf_tables.inc";
2652 outFormat
= OF_TABLE
;
2654 fprintf (stderr
, "Error: Invalid command '%s'.\n", argv
[1]);
2660 equalize
= DEFAULT_EQUALIZE
;
2661 surface
= DEFAULT_SURFACE
;
2662 limit
= DEFAULT_LIMIT
;
2663 truncSize
= DEFAULT_TRUNCSIZE
;
2664 model
= DEFAULT_HEAD_MODEL
;
2665 radius
= DEFAULT_CUSTOM_RADIUS
;
2666 while (argi
< argc
) {
2667 if (strncmp (argv
[argi
], "-r=", 3) == 0) {
2668 outRate
= strtoul (& argv
[argi
] [3], & end
, 10);
2669 if ((end
[0] != '\0') || (outRate
< MIN_RATE
) || (outRate
> MAX_RATE
)) {
2670 fprintf (stderr
, "Error: Expected a value from %u to %u for '-r'.\n", MIN_RATE
, MAX_RATE
);
2673 } else if (strncmp (argv
[argi
], "-f=", 3) == 0) {
2674 fftSize
= strtoul (& argv
[argi
] [3], & end
, 10);
2675 if ((end
[0] != '\0') || (fftSize
& (fftSize
- 1)) || (fftSize
< MIN_FFTSIZE
) || (fftSize
> MAX_FFTSIZE
)) {
2676 fprintf (stderr
, "Error: Expected a power-of-two value from %u to %u for '-f'.\n", MIN_FFTSIZE
, MAX_FFTSIZE
);
2679 } else if (strncmp (argv
[argi
], "-e=", 3) == 0) {
2680 if (strcmp (& argv
[argi
] [3], "on") == 0) {
2682 } else if (strcmp (& argv
[argi
] [3], "off") == 0) {
2685 fprintf (stderr
, "Error: Expected 'on' or 'off' for '-e'.\n");
2688 } else if (strncmp (argv
[argi
], "-s=", 3) == 0) {
2689 if (strcmp (& argv
[argi
] [3], "on") == 0) {
2691 } else if (strcmp (& argv
[argi
] [3], "off") == 0) {
2694 fprintf (stderr
, "Error: Expected 'on' or 'off' for '-s'.\n");
2697 } else if (strncmp (argv
[argi
], "-l=", 3) == 0) {
2698 if (strcmp (& argv
[argi
] [3], "none") == 0) {
2701 limit
= strtod (& argv
[argi
] [3], & end
);
2702 if ((end
[0] != '\0') || (limit
< MIN_LIMIT
) || (limit
> MAX_LIMIT
)) {
2703 fprintf (stderr
, "Error: Expected 'none' or a value from %.2f to %.2f for '-l'.\n", MIN_LIMIT
, MAX_LIMIT
);
2707 } else if (strncmp (argv
[argi
], "-w=", 3) == 0) {
2708 truncSize
= strtoul (& argv
[argi
] [3], & end
, 10);
2709 if ((end
[0] != '\0') || (truncSize
< MIN_TRUNCSIZE
) || (truncSize
> MAX_TRUNCSIZE
) || (truncSize
% MOD_TRUNCSIZE
)) {
2710 fprintf (stderr
, "Error: Expected a value from %u to %u in multiples of %u for '-w'.\n", MIN_TRUNCSIZE
, MAX_TRUNCSIZE
, MOD_TRUNCSIZE
);
2713 } else if (strncmp (argv
[argi
], "-d=", 3) == 0) {
2714 if (strcmp (& argv
[argi
] [3], "dataset") == 0) {
2716 } else if (strcmp (& argv
[argi
] [3], "sphere") == 0) {
2719 fprintf (stderr
, "Error: Expected 'dataset' or 'sphere' for '-d'.\n");
2722 } else if (strncmp (argv
[argi
], "-c=", 3) == 0) {
2723 radius
= strtod (& argv
[argi
] [3], & end
);
2724 if ((end
[0] != '\0') || (radius
< MIN_CUSTOM_RADIUS
) || (radius
> MAX_CUSTOM_RADIUS
)) {
2725 fprintf (stderr
, "Error: Expected a value from %.2f to %.2f for '-c'.\n", MIN_CUSTOM_RADIUS
, MAX_CUSTOM_RADIUS
);
2728 } else if (strncmp (argv
[argi
], "-i=", 3) == 0) {
2729 inName
= & argv
[argi
] [3];
2730 } else if (strncmp (argv
[argi
], "-o=", 3) == 0) {
2731 outName
= & argv
[argi
] [3];
2733 fprintf (stderr
, "Error: Invalid option '%s'.\n", argv
[argi
]);
2738 if (! ProcessDefinition (inName
, outRate
, fftSize
, equalize
, surface
, limit
, truncSize
, model
, radius
, outFormat
, outName
))
2740 fprintf (stdout
, "Operation completed.\n");