2 * HRTF utility for producing and demonstrating the process of creating an
3 * OpenAL Soft compatible HRIR data set.
5 * Copyright (C) 2011-2017 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
80 // Rely (if naively) on OpenAL's header for the types used for serialization.
85 #define M_PI (3.14159265358979323846)
89 #define HUGE_VAL (1.0 / 0.0)
94 #define WIN32_LEAN_AND_MEAN
97 static char *ToUTF8(const wchar_t *from
)
101 if((len
=WideCharToMultiByte(CP_UTF8
, 0, from
, -1, NULL
, 0, NULL
, NULL
)) > 0)
103 out
= calloc(sizeof(*out
), len
);
104 WideCharToMultiByte(CP_UTF8
, 0, from
, -1, out
, len
, NULL
, NULL
);
110 static WCHAR
*FromUTF8(const char *str
)
115 if((len
=MultiByteToWideChar(CP_UTF8
, 0, str
, -1, NULL
, 0)) > 0)
117 out
= calloc(sizeof(WCHAR
), len
);
118 MultiByteToWideChar(CP_UTF8
, 0, str
, -1, out
, len
);
125 static FILE *my_fopen(const char *fname
, const char *mode
)
127 WCHAR
*wname
=NULL
, *wmode
=NULL
;
130 wname
= FromUTF8(fname
);
131 wmode
= FromUTF8(mode
);
133 fprintf(stderr
, "Failed to convert UTF-8 filename: \"%s\"\n", fname
);
135 fprintf(stderr
, "Failed to convert UTF-8 mode: \"%s\"\n", mode
);
137 file
= _wfopen(wname
, wmode
);
144 #define fopen my_fopen
147 // The epsilon used to maintain signal stability.
148 #define EPSILON (1e-9)
150 // Constants for accessing the token reader's ring buffer.
151 #define TR_RING_BITS (16)
152 #define TR_RING_SIZE (1 << TR_RING_BITS)
153 #define TR_RING_MASK (TR_RING_SIZE - 1)
155 // The token reader's load interval in bytes.
156 #define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
158 // The maximum identifier length used when processing the data set
160 #define MAX_IDENT_LEN (16)
162 // The maximum path length used when processing filenames.
163 #define MAX_PATH_LEN (256)
165 // The limits for the sample 'rate' metric in the data set definition and for
167 #define MIN_RATE (32000)
168 #define MAX_RATE (96000)
170 // The limits for the HRIR 'points' metric in the data set definition.
171 #define MIN_POINTS (16)
172 #define MAX_POINTS (8192)
174 // The limits to the number of 'azimuths' listed in the data set definition.
175 #define MIN_EV_COUNT (5)
176 #define MAX_EV_COUNT (128)
178 // The limits for each of the 'azimuths' listed in the data set definition.
179 #define MIN_AZ_COUNT (1)
180 #define MAX_AZ_COUNT (128)
182 // The limits for the listener's head 'radius' in the data set definition.
183 #define MIN_RADIUS (0.05)
184 #define MAX_RADIUS (0.15)
186 // The limits for the 'distance' from source to listener in the definition
188 #define MIN_DISTANCE (0.5)
189 #define MAX_DISTANCE (2.5)
191 // The maximum number of channels that can be addressed for a WAVE file
192 // source listed in the data set definition.
193 #define MAX_WAVE_CHANNELS (65535)
195 // The limits to the byte size for a binary source listed in the definition
197 #define MIN_BIN_SIZE (2)
198 #define MAX_BIN_SIZE (4)
200 // The minimum number of significant bits for binary sources listed in the
201 // data set definition. The maximum is calculated from the byte size.
202 #define MIN_BIN_BITS (16)
204 // The limits to the number of significant bits for an ASCII source listed in
205 // the data set definition.
206 #define MIN_ASCII_BITS (16)
207 #define MAX_ASCII_BITS (32)
209 // The limits to the FFT window size override on the command line.
210 #define MIN_FFTSIZE (65536)
211 #define MAX_FFTSIZE (131072)
213 // The limits to the equalization range limit on the command line.
214 #define MIN_LIMIT (2.0)
215 #define MAX_LIMIT (120.0)
217 // The limits to the truncation window size on the command line.
218 #define MIN_TRUNCSIZE (16)
219 #define MAX_TRUNCSIZE (512)
221 // The limits to the custom head radius on the command line.
222 #define MIN_CUSTOM_RADIUS (0.05)
223 #define MAX_CUSTOM_RADIUS (0.15)
225 // The truncation window size must be a multiple of the below value to allow
226 // for vectorized convolution.
227 #define MOD_TRUNCSIZE (8)
229 // The defaults for the command line options.
230 #define DEFAULT_FFTSIZE (65536)
231 #define DEFAULT_EQUALIZE (1)
232 #define DEFAULT_SURFACE (1)
233 #define DEFAULT_LIMIT (24.0)
234 #define DEFAULT_TRUNCSIZE (32)
235 #define DEFAULT_HEAD_MODEL (HM_DATASET)
236 #define DEFAULT_CUSTOM_RADIUS (0.0)
238 // The four-character-codes for RIFF/RIFX WAVE file chunks.
239 #define FOURCC_RIFF (0x46464952) // 'RIFF'
240 #define FOURCC_RIFX (0x58464952) // 'RIFX'
241 #define FOURCC_WAVE (0x45564157) // 'WAVE'
242 #define FOURCC_FMT (0x20746D66) // 'fmt '
243 #define FOURCC_DATA (0x61746164) // 'data'
244 #define FOURCC_LIST (0x5453494C) // 'LIST'
245 #define FOURCC_WAVL (0x6C766177) // 'wavl'
246 #define FOURCC_SLNT (0x746E6C73) // 'slnt'
248 // The supported wave formats.
249 #define WAVE_FORMAT_PCM (0x0001)
250 #define WAVE_FORMAT_IEEE_FLOAT (0x0003)
251 #define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
253 // The maximum propagation delay value supported by OpenAL Soft.
254 #define MAX_HRTD (63.0)
256 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
257 // response protocol 01.
258 #define MHR_FORMAT ("MinPHR01")
260 #define MHR_FORMAT_EXPERIMENTAL ("MinPHRTEMPDONOTUSE")
262 // Sample and channel type enum values
263 typedef enum SampleTypeT
{
268 typedef enum ChannelTypeT
{
273 // Byte order for the serialization routines.
274 typedef enum ByteOrderT
{
280 // Source format for the references listed in the data set definition.
281 typedef enum SourceFormatT
{
283 SF_WAVE
, // RIFF/RIFX WAVE file.
284 SF_BIN_LE
, // Little-endian binary file.
285 SF_BIN_BE
, // Big-endian binary file.
286 SF_ASCII
// ASCII text file.
289 // Element types for the references listed in the data set definition.
290 typedef enum ElementTypeT
{
292 ET_INT
, // Integer elements.
293 ET_FP
// Floating-point elements.
296 // Head model used for calculating the impulse delays.
297 typedef enum HeadModelT
{
299 HM_DATASET
, // Measure the onset from the dataset.
300 HM_SPHERE
// Calculate the onset using a spherical head model.
303 // Desired output format from the command line.
304 typedef enum OutputFormatT
{
306 OF_MHR
// OpenAL Soft MHR data set file.
309 // Unsigned integer type.
310 typedef unsigned int uint
;
312 // Serialization types. The trailing digit indicates the number of bits.
313 typedef ALubyte uint8
;
315 typedef ALuint uint32
;
316 typedef ALuint64SOFT uint64
;
318 // Token reader state for parsing the data set definition.
319 typedef struct TokenReaderT
{
324 char mRing
[TR_RING_SIZE
];
329 // Source reference state used when loading sources.
330 typedef struct SourceRefT
{
331 SourceFormatT mFormat
;
338 char mPath
[MAX_PATH_LEN
+1];
341 // The HRIR metrics and data set used when loading, processing, and storing
342 // the resulting HRTF.
343 typedef struct HrirDataT
{
345 SampleTypeT mSampleType
;
346 ChannelTypeT mChannelType
;
353 uint mAzCount
[MAX_EV_COUNT
];
354 uint mEvOffset
[MAX_EV_COUNT
];
362 // The resampler metrics and FIR filter.
363 typedef struct ResamplerT
{
369 /****************************************
370 *** Complex number type and routines ***
371 ****************************************/
377 static Complex
MakeComplex(double r
, double i
)
379 Complex c
= { r
, i
};
383 static Complex
c_add(Complex a
, Complex b
)
386 r
.Real
= a
.Real
+ b
.Real
;
387 r
.Imag
= a
.Imag
+ b
.Imag
;
391 static Complex
c_sub(Complex a
, Complex b
)
394 r
.Real
= a
.Real
- b
.Real
;
395 r
.Imag
= a
.Imag
- b
.Imag
;
399 static Complex
c_mul(Complex a
, Complex b
)
402 r
.Real
= a
.Real
*b
.Real
- a
.Imag
*b
.Imag
;
403 r
.Imag
= a
.Imag
*b
.Real
+ a
.Real
*b
.Imag
;
407 static Complex
c_muls(Complex a
, double s
)
415 static double c_abs(Complex a
)
417 return sqrt(a
.Real
*a
.Real
+ a
.Imag
*a
.Imag
);
420 static Complex
c_exp(Complex a
)
423 double e
= exp(a
.Real
);
424 r
.Real
= e
* cos(a
.Imag
);
425 r
.Imag
= e
* sin(a
.Imag
);
429 /*****************************
430 *** Token reader routines ***
431 *****************************/
433 /* Whitespace is not significant. It can process tokens as identifiers, numbers
434 * (integer and floating-point), strings, and operators. Strings must be
435 * encapsulated by double-quotes and cannot span multiple lines.
438 // Setup the reader on the given file. The filename can be NULL if no error
439 // output is desired.
440 static void TrSetup(FILE *fp
, const char *filename
, TokenReaderT
*tr
)
442 const char *name
= NULL
;
446 const char *slash
= strrchr(filename
, '/');
449 const char *bslash
= strrchr(slash
+1, '\\');
450 if(bslash
) name
= bslash
+1;
455 const char *bslash
= strrchr(filename
, '\\');
456 if(bslash
) name
= bslash
+1;
457 else name
= filename
;
469 // Prime the reader's ring buffer, and return a result indicating that there
470 // is text to process.
471 static int TrLoad(TokenReaderT
*tr
)
473 size_t toLoad
, in
, count
;
475 toLoad
= TR_RING_SIZE
- (tr
->mIn
- tr
->mOut
);
476 if(toLoad
>= TR_LOAD_SIZE
&& !feof(tr
->mFile
))
478 // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
479 toLoad
= TR_LOAD_SIZE
;
480 in
= tr
->mIn
&TR_RING_MASK
;
481 count
= TR_RING_SIZE
- in
;
484 tr
->mIn
+= fread(&tr
->mRing
[in
], 1, count
, tr
->mFile
);
485 tr
->mIn
+= fread(&tr
->mRing
[0], 1, toLoad
-count
, tr
->mFile
);
488 tr
->mIn
+= fread(&tr
->mRing
[in
], 1, toLoad
, tr
->mFile
);
490 if(tr
->mOut
>= TR_RING_SIZE
)
492 tr
->mOut
-= TR_RING_SIZE
;
493 tr
->mIn
-= TR_RING_SIZE
;
496 if(tr
->mIn
> tr
->mOut
)
501 // Error display routine. Only displays when the base name is not NULL.
502 static void TrErrorVA(const TokenReaderT
*tr
, uint line
, uint column
, const char *format
, va_list argPtr
)
506 fprintf(stderr
, "Error (%s:%u:%u): ", tr
->mName
, line
, column
);
507 vfprintf(stderr
, format
, argPtr
);
510 // Used to display an error at a saved line/column.
511 static void TrErrorAt(const TokenReaderT
*tr
, uint line
, uint column
, const char *format
, ...)
515 va_start(argPtr
, format
);
516 TrErrorVA(tr
, line
, column
, format
, argPtr
);
520 // Used to display an error at the current line/column.
521 static void TrError(const TokenReaderT
*tr
, const char *format
, ...)
525 va_start(argPtr
, format
);
526 TrErrorVA(tr
, tr
->mLine
, tr
->mColumn
, format
, argPtr
);
530 // Skips to the next line.
531 static void TrSkipLine(TokenReaderT
*tr
)
537 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
549 // Skips to the next token.
550 static int TrSkipWhitespace(TokenReaderT
*tr
)
556 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
576 // Get the line and/or column of the next token (or the end of input).
577 static void TrIndication(TokenReaderT
*tr
, uint
*line
, uint
*column
)
579 TrSkipWhitespace(tr
);
580 if(line
) *line
= tr
->mLine
;
581 if(column
) *column
= tr
->mColumn
;
584 // Checks to see if a token is the given operator. It does not display any
585 // errors and will not proceed to the next token.
586 static int TrIsOperator(TokenReaderT
*tr
, const char *op
)
591 if(!TrSkipWhitespace(tr
))
595 while(op
[len
] != '\0' && out
< tr
->mIn
)
597 ch
= tr
->mRing
[out
&TR_RING_MASK
];
598 if(ch
!= op
[len
]) break;
607 /* The TrRead*() routines obtain the value of a matching token type. They
608 * display type, form, and boundary errors and will proceed to the next
612 // Reads and validates an identifier token.
613 static int TrReadIdent(TokenReaderT
*tr
, const uint maxLen
, char *ident
)
619 if(TrSkipWhitespace(tr
))
622 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
623 if(ch
== '_' || isalpha(ch
))
633 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
634 } while(ch
== '_' || isdigit(ch
) || isalpha(ch
));
642 TrErrorAt(tr
, tr
->mLine
, col
, "Identifier is too long.\n");
646 TrErrorAt(tr
, tr
->mLine
, col
, "Expected an identifier.\n");
650 // Reads and validates (including bounds) an integer token.
651 static int TrReadInt(TokenReaderT
*tr
, const int loBound
, const int hiBound
, int *value
)
653 uint col
, digis
, len
;
657 if(TrSkipWhitespace(tr
))
661 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
662 if(ch
== '+' || ch
== '-')
671 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
672 if(!isdigit(ch
)) break;
680 if(digis
> 0 && ch
!= '.' && !isalpha(ch
))
684 TrErrorAt(tr
, tr
->mLine
, col
, "Integer is too long.");
688 *value
= strtol(temp
, NULL
, 10);
689 if(*value
< loBound
|| *value
> hiBound
)
691 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a value from %d to %d.\n", loBound
, hiBound
);
697 TrErrorAt(tr
, tr
->mLine
, col
, "Expected an integer.\n");
701 // Reads and validates (including bounds) a float token.
702 static int TrReadFloat(TokenReaderT
*tr
, const double loBound
, const double hiBound
, double *value
)
704 uint col
, digis
, len
;
708 if(TrSkipWhitespace(tr
))
712 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
713 if(ch
== '+' || ch
== '-')
723 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
724 if(!isdigit(ch
)) break;
740 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
741 if(!isdigit(ch
)) break;
750 if(ch
== 'E' || ch
== 'e')
757 if(ch
== '+' || ch
== '-')
766 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
767 if(!isdigit(ch
)) break;
776 if(digis
> 0 && ch
!= '.' && !isalpha(ch
))
780 TrErrorAt(tr
, tr
->mLine
, col
, "Float is too long.");
784 *value
= strtod(temp
, NULL
);
785 if(*value
< loBound
|| *value
> hiBound
)
787 TrErrorAt (tr
, tr
->mLine
, col
, "Expected a value from %f to %f.\n", loBound
, hiBound
);
796 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a float.\n");
800 // Reads and validates a string token.
801 static int TrReadString(TokenReaderT
*tr
, const uint maxLen
, char *text
)
807 if(TrSkipWhitespace(tr
))
810 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
817 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
823 TrErrorAt (tr
, tr
->mLine
, col
, "Unterminated string at end of line.\n");
832 tr
->mColumn
+= 1 + len
;
833 TrErrorAt(tr
, tr
->mLine
, col
, "Unterminated string at end of input.\n");
836 tr
->mColumn
+= 2 + len
;
839 TrErrorAt (tr
, tr
->mLine
, col
, "String is too long.\n");
846 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a string.\n");
850 // Reads and validates the given operator.
851 static int TrReadOperator(TokenReaderT
*tr
, const char *op
)
857 if(TrSkipWhitespace(tr
))
861 while(op
[len
] != '\0' && TrLoad(tr
))
863 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
864 if(ch
!= op
[len
]) break;
872 TrErrorAt(tr
, tr
->mLine
, col
, "Expected '%s' operator.\n", op
);
876 /* Performs a string substitution. Any case-insensitive occurrences of the
877 * pattern string are replaced with the replacement string. The result is
878 * truncated if necessary.
880 static int StrSubst(const char *in
, const char *pat
, const char *rep
, const size_t maxLen
, char *out
)
882 size_t inLen
, patLen
, repLen
;
887 patLen
= strlen(pat
);
888 repLen
= strlen(rep
);
892 while(si
< inLen
&& di
< maxLen
)
894 if(patLen
<= inLen
-si
)
896 if(strncasecmp(&in
[si
], pat
, patLen
) == 0)
898 if(repLen
> maxLen
-di
)
900 repLen
= maxLen
- di
;
903 strncpy(&out
[di
], rep
, repLen
);
919 /*********************
920 *** Math routines ***
921 *********************/
923 // Provide missing math routines for MSVC versions < 1800 (Visual Studio 2013).
924 #if defined(_MSC_VER) && _MSC_VER < 1800
925 static double round(double val
)
928 return ceil(val
-0.5);
929 return floor(val
+0.5);
932 static double fmin(double a
, double b
)
934 return (a
<b
) ? a
: b
;
937 static double fmax(double a
, double b
)
939 return (a
>b
) ? a
: b
;
943 // Simple clamp routine.
944 static double Clamp(const double val
, const double lower
, const double upper
)
946 return fmin(fmax(val
, lower
), upper
);
949 // Performs linear interpolation.
950 static double Lerp(const double a
, const double b
, const double f
)
952 return a
+ (f
* (b
- a
));
955 static inline uint
dither_rng(uint
*seed
)
957 *seed
= (*seed
* 96314165) + 907633515;
961 // Performs a triangular probability density function dither. It assumes the
962 // input sample is already scaled.
963 static inline double TpdfDither(const double in
, uint
*seed
)
965 static const double PRNG_SCALE
= 1.0 / UINT_MAX
;
968 prn0
= dither_rng(seed
);
969 prn1
= dither_rng(seed
);
970 return round(in
+ (prn0
*PRNG_SCALE
- prn1
*PRNG_SCALE
));
973 // Allocates an array of doubles.
974 static double *CreateArray(size_t n
)
979 a
= calloc(n
, sizeof(double));
982 fprintf(stderr
, "Error: Out of memory.\n");
988 // Frees an array of doubles.
989 static void DestroyArray(double *a
)
992 /* Fast Fourier transform routines. The number of points must be a power of
993 * two. In-place operation is possible only if both the real and imaginary
994 * parts are in-place together.
997 // Performs bit-reversal ordering.
998 static void FftArrange(const uint n
, const Complex
*in
, Complex
*out
)
1004 // Handle in-place arrangement.
1006 for(k
= 0;k
< n
;k
++)
1010 Complex temp
= in
[rk
];
1022 // Handle copy arrangement.
1024 for(k
= 0;k
< n
;k
++)
1035 // Performs the summation.
1036 static void FftSummation(const int n
, const double s
, Complex
*cplx
)
1043 for(m
= 1, m2
= 2;m
< n
; m
<<= 1, m2
<<= 1)
1045 // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
1046 double sm
= sin(0.5 * pi
/ m
);
1047 Complex v
= MakeComplex(-2.0*sm
*sm
, -sin(pi
/ m
));
1048 Complex w
= MakeComplex(1.0, 0.0);
1049 for(i
= 0;i
< m
;i
++)
1051 for(k
= i
;k
< n
;k
+= m2
)
1055 t
= c_mul(w
, cplx
[mk
]);
1056 cplx
[mk
] = c_sub(cplx
[k
], t
);
1057 cplx
[k
] = c_add(cplx
[k
], t
);
1059 w
= c_add(w
, c_mul(v
, w
));
1064 // Performs a forward FFT.
1065 static void FftForward(const uint n
, const Complex
*in
, Complex
*out
)
1067 FftArrange(n
, in
, out
);
1068 FftSummation(n
, 1.0, out
);
1071 // Performs an inverse FFT.
1072 static void FftInverse(const uint n
, const Complex
*in
, Complex
*out
)
1077 FftArrange(n
, in
, out
);
1078 FftSummation(n
, -1.0, out
);
1080 for(i
= 0;i
< n
;i
++)
1081 out
[i
] = c_muls(out
[i
], f
);
1084 /* Calculate the complex helical sequence (or discrete-time analytical signal)
1085 * of the given input using the Hilbert transform. Given the natural logarithm
1086 * of a signal's magnitude response, the imaginary components can be used as
1087 * the angles for minimum-phase reconstruction.
1089 static void Hilbert(const uint n
, const Complex
*in
, Complex
*out
)
1095 // Handle in-place operation.
1096 for(i
= 0;i
< n
;i
++)
1101 // Handle copy operation.
1102 for(i
= 0;i
< n
;i
++)
1103 out
[i
] = MakeComplex(in
[i
].Real
, 0.0);
1105 FftInverse(n
, out
, out
);
1106 for(i
= 1;i
< (n
+1)/2;i
++)
1107 out
[i
] = c_muls(out
[i
], 2.0);
1108 /* Increment i if n is even. */
1111 out
[i
] = MakeComplex(0.0, 0.0);
1112 FftForward(n
, out
, out
);
1115 /* Calculate the magnitude response of the given input. This is used in
1116 * place of phase decomposition, since the phase residuals are discarded for
1117 * minimum phase reconstruction. The mirrored half of the response is also
1120 static void MagnitudeResponse(const uint n
, const Complex
*in
, double *out
)
1122 const uint m
= 1 + (n
/ 2);
1124 for(i
= 0;i
< m
;i
++)
1125 out
[i
] = fmax(c_abs(in
[i
]), EPSILON
);
1128 /* Apply a range limit (in dB) to the given magnitude response. This is used
1129 * to adjust the effects of the diffuse-field average on the equalization
1132 static void LimitMagnitudeResponse(const uint n
, const double limit
, const double *in
, double *out
)
1134 const uint m
= 1 + (n
/ 2);
1136 uint i
, lower
, upper
;
1139 halfLim
= limit
/ 2.0;
1140 // Convert the response to dB.
1141 for(i
= 0;i
< m
;i
++)
1142 out
[i
] = 20.0 * log10(in
[i
]);
1143 // Use six octaves to calculate the average magnitude of the signal.
1144 lower
= ((uint
)ceil(n
/ pow(2.0, 8.0))) - 1;
1145 upper
= ((uint
)floor(n
/ pow(2.0, 2.0))) - 1;
1147 for(i
= lower
;i
<= upper
;i
++)
1149 ave
/= upper
- lower
+ 1;
1150 // Keep the response within range of the average magnitude.
1151 for(i
= 0;i
< m
;i
++)
1152 out
[i
] = Clamp(out
[i
], ave
- halfLim
, ave
+ halfLim
);
1153 // Convert the response back to linear magnitude.
1154 for(i
= 0;i
< m
;i
++)
1155 out
[i
] = pow(10.0, out
[i
] / 20.0);
1158 /* Reconstructs the minimum-phase component for the given magnitude response
1159 * of a signal. This is equivalent to phase recomposition, sans the missing
1160 * residuals (which were discarded). The mirrored half of the response is
1163 static void MinimumPhase(const uint n
, const double *in
, Complex
*out
)
1165 const uint m
= 1 + (n
/ 2);
1169 mags
= CreateArray(n
);
1170 for(i
= 0;i
< m
;i
++)
1172 mags
[i
] = fmax(EPSILON
, in
[i
]);
1173 out
[i
] = MakeComplex(log(mags
[i
]), 0.0);
1177 mags
[i
] = mags
[n
- i
];
1178 out
[i
] = out
[n
- i
];
1180 Hilbert(n
, out
, out
);
1181 // Remove any DC offset the filter has.
1183 for(i
= 0;i
< n
;i
++)
1185 Complex a
= c_exp(MakeComplex(0.0, out
[i
].Imag
));
1186 out
[i
] = c_mul(MakeComplex(mags
[i
], 0.0), a
);
1192 /***************************
1193 *** Resampler functions ***
1194 ***************************/
1196 /* This is the normalized cardinal sine (sinc) function.
1198 * sinc(x) = { 1, x = 0
1199 * { sin(pi x) / (pi x), otherwise.
1201 static double Sinc(const double x
)
1203 if(fabs(x
) < EPSILON
)
1205 return sin(M_PI
* x
) / (M_PI
* x
);
1208 /* The zero-order modified Bessel function of the first kind, used for the
1211 * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
1212 * = sum_{k=0}^inf ((x / 2)^k / k!)^2
1214 static double BesselI_0(const double x
)
1216 double term
, sum
, x2
, y
, last_sum
;
1219 // Start at k=1 since k=0 is trivial.
1225 // Let the integration converge until the term of the sum is no longer
1233 } while(sum
!= last_sum
);
1237 /* Calculate a Kaiser window from the given beta value and a normalized k
1240 * w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
1243 * Where k can be calculated as:
1245 * k = i / l, where -l <= i <= l.
1249 * k = 2 i / M - 1, where 0 <= i <= M.
1251 static double Kaiser(const double b
, const double k
)
1253 if(!(k
>= -1.0 && k
<= 1.0))
1255 return BesselI_0(b
* sqrt(1.0 - k
*k
)) / BesselI_0(b
);
1258 // Calculates the greatest common divisor of a and b.
1259 static uint
Gcd(uint x
, uint y
)
1270 /* Calculates the size (order) of the Kaiser window. Rejection is in dB and
1271 * the transition width is normalized frequency (0.5 is nyquist).
1273 * M = { ceil((r - 7.95) / (2.285 2 pi f_t)), r > 21
1274 * { ceil(5.79 / 2 pi f_t), r <= 21.
1277 static uint
CalcKaiserOrder(const double rejection
, const double transition
)
1279 double w_t
= 2.0 * M_PI
* transition
;
1280 if(rejection
> 21.0)
1281 return (uint
)ceil((rejection
- 7.95) / (2.285 * w_t
));
1282 return (uint
)ceil(5.79 / w_t
);
1285 // Calculates the beta value of the Kaiser window. Rejection is in dB.
1286 static double CalcKaiserBeta(const double rejection
)
1288 if(rejection
> 50.0)
1289 return 0.1102 * (rejection
- 8.7);
1290 if(rejection
>= 21.0)
1291 return (0.5842 * pow(rejection
- 21.0, 0.4)) +
1292 (0.07886 * (rejection
- 21.0));
1296 /* Calculates a point on the Kaiser-windowed sinc filter for the given half-
1297 * width, beta, gain, and cutoff. The point is specified in non-normalized
1298 * samples, from 0 to M, where M = (2 l + 1).
1300 * w(k) 2 p f_t sinc(2 f_t x)
1302 * x -- centered sample index (i - l)
1303 * k -- normalized and centered window index (x / l)
1304 * w(k) -- window function (Kaiser)
1305 * p -- gain compensation factor when sampling
1306 * f_t -- normalized center frequency (or cutoff; 0.5 is nyquist)
1308 static double SincFilter(const int l
, const double b
, const double gain
, const double cutoff
, const int i
)
1310 return Kaiser(b
, (double)(i
- l
) / l
) * 2.0 * gain
* cutoff
* Sinc(2.0 * cutoff
* (i
- l
));
1313 /* This is a polyphase sinc-filtered resampler.
1315 * Upsample Downsample
1317 * p/q = 3/2 p/q = 3/5
1319 * M-+-+-+-> M-+-+-+->
1320 * -------------------+ ---------------------+
1321 * p s * f f f f|f| | p s * f f f f f |
1322 * | 0 * 0 0 0|0|0 | | 0 * 0 0 0 0|0| |
1323 * v 0 * 0 0|0|0 0 | v 0 * 0 0 0|0|0 |
1324 * s * f|f|f f f | s * f f|f|f f |
1325 * 0 * |0|0 0 0 0 | 0 * 0|0|0 0 0 |
1326 * --------+=+--------+ 0 * |0|0 0 0 0 |
1327 * d . d .|d|. d . d ----------+=+--------+
1328 * d . . . .|d|. . . .
1332 * P_f(i,j) = q i mod p + pj
1333 * P_s(i,j) = floor(q i / p) - j
1334 * d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
1335 * { f[P_f(i,j)] s[P_s(i,j)], P_f(i,j) < M
1336 * { 0, P_f(i,j) >= M. }
1339 // Calculate the resampling metrics and build the Kaiser-windowed sinc filter
1340 // that's used to cut frequencies above the destination nyquist.
1341 static void ResamplerSetup(ResamplerT
*rs
, const uint srcRate
, const uint dstRate
)
1343 double cutoff
, width
, beta
;
1347 gcd
= Gcd(srcRate
, dstRate
);
1348 rs
->mP
= dstRate
/ gcd
;
1349 rs
->mQ
= srcRate
/ gcd
;
1350 /* The cutoff is adjusted by half the transition width, so the transition
1351 * ends before the nyquist (0.5). Both are scaled by the downsampling
1356 cutoff
= 0.475 / rs
->mP
;
1357 width
= 0.05 / rs
->mP
;
1361 cutoff
= 0.475 / rs
->mQ
;
1362 width
= 0.05 / rs
->mQ
;
1364 // A rejection of -180 dB is used for the stop band. Round up when
1365 // calculating the left offset to avoid increasing the transition width.
1366 l
= (CalcKaiserOrder(180.0, width
)+1) / 2;
1367 beta
= CalcKaiserBeta(180.0);
1370 rs
->mF
= CreateArray(rs
->mM
);
1371 for(i
= 0;i
< ((int)rs
->mM
);i
++)
1372 rs
->mF
[i
] = SincFilter((int)l
, beta
, rs
->mP
, cutoff
, i
);
1375 // Clean up after the resampler.
1376 static void ResamplerClear(ResamplerT
*rs
)
1378 DestroyArray(rs
->mF
);
1382 // Perform the upsample-filter-downsample resampling operation using a
1383 // polyphase filter implementation.
1384 static void ResamplerRun(ResamplerT
*rs
, const uint inN
, const double *in
, const uint outN
, double *out
)
1386 const uint p
= rs
->mP
, q
= rs
->mQ
, m
= rs
->mM
, l
= rs
->mL
;
1387 const double *f
= rs
->mF
;
1395 // Handle in-place operation.
1397 work
= CreateArray(outN
);
1400 // Resample the input.
1401 for(i
= 0;i
< outN
;i
++)
1404 // Input starts at l to compensate for the filter delay. This will
1405 // drop any build-up from the first half of the filter.
1406 j_f
= (l
+ (q
* i
)) % p
;
1407 j_s
= (l
+ (q
* i
)) / p
;
1410 // Only take input when 0 <= j_s < inN. This single unsigned
1411 // comparison catches both cases.
1413 r
+= f
[j_f
] * in
[j_s
];
1419 // Clean up after in-place operation.
1422 for(i
= 0;i
< outN
;i
++)
1428 /*************************
1429 *** File source input ***
1430 *************************/
1432 // Read a binary value of the specified byte order and byte size from a file,
1433 // storing it as a 32-bit unsigned integer.
1434 static int ReadBin4(FILE *fp
, const char *filename
, const ByteOrderT order
, const uint bytes
, uint32
*out
)
1440 if(fread(in
, 1, bytes
, fp
) != bytes
)
1442 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1449 for(i
= 0;i
< bytes
;i
++)
1450 accum
= (accum
<<8) | in
[bytes
- i
- 1];
1453 for(i
= 0;i
< bytes
;i
++)
1454 accum
= (accum
<<8) | in
[i
];
1463 // Read a binary value of the specified byte order from a file, storing it as
1464 // a 64-bit unsigned integer.
1465 static int ReadBin8(FILE *fp
, const char *filename
, const ByteOrderT order
, uint64
*out
)
1471 if(fread(in
, 1, 8, fp
) != 8)
1473 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1480 for(i
= 0;i
< 8;i
++)
1481 accum
= (accum
<<8) | in
[8 - i
- 1];
1484 for(i
= 0;i
< 8;i
++)
1485 accum
= (accum
<<8) | in
[i
];
1494 /* Read a binary value of the specified type, byte order, and byte size from
1495 * a file, converting it to a double. For integer types, the significant
1496 * bits are used to normalize the result. The sign of bits determines
1497 * whether they are padded toward the MSB (negative) or LSB (positive).
1498 * Floating-point types are not normalized.
1500 static int ReadBinAsDouble(FILE *fp
, const char *filename
, const ByteOrderT order
, const ElementTypeT type
, const uint bytes
, const int bits
, double *out
)
1515 if(!ReadBin8(fp
, filename
, order
, &v8
.ui
))
1522 if(!ReadBin4(fp
, filename
, order
, bytes
, &v4
.ui
))
1529 v4
.ui
>>= (8*bytes
) - ((uint
)bits
);
1531 v4
.ui
&= (0xFFFFFFFF >> (32+bits
));
1533 if(v4
.ui
&(uint
)(1<<(abs(bits
)-1)))
1534 v4
.ui
|= (0xFFFFFFFF << abs (bits
));
1535 *out
= v4
.i
/ (double)(1<<(abs(bits
)-1));
1541 /* Read an ascii value of the specified type from a file, converting it to a
1542 * double. For integer types, the significant bits are used to normalize the
1543 * result. The sign of the bits should always be positive. This also skips
1544 * up to one separator character before the element itself.
1546 static int ReadAsciiAsDouble(TokenReaderT
*tr
, const char *filename
, const ElementTypeT type
, const uint bits
, double *out
)
1548 if(TrIsOperator(tr
, ","))
1549 TrReadOperator(tr
, ",");
1550 else if(TrIsOperator(tr
, ":"))
1551 TrReadOperator(tr
, ":");
1552 else if(TrIsOperator(tr
, ";"))
1553 TrReadOperator(tr
, ";");
1554 else if(TrIsOperator(tr
, "|"))
1555 TrReadOperator(tr
, "|");
1559 if(!TrReadFloat(tr
, -HUGE_VAL
, HUGE_VAL
, out
))
1561 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1568 if(!TrReadInt(tr
, -(1<<(bits
-1)), (1<<(bits
-1))-1, &v
))
1570 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1573 *out
= v
/ (double)((1<<(bits
-1))-1);
1578 // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
1579 // the source parameters and data set metrics.
1580 static int ReadWaveFormat(FILE *fp
, const ByteOrderT order
, const uint hrirRate
, SourceRefT
*src
)
1582 uint32 fourCC
, chunkSize
;
1583 uint32 format
, channels
, rate
, dummy
, block
, size
, bits
;
1588 fseek (fp
, (long) chunkSize
, SEEK_CUR
);
1589 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1590 !ReadBin4(fp
, src
->mPath
, order
, 4, &chunkSize
))
1592 } while(fourCC
!= FOURCC_FMT
);
1593 if(!ReadBin4(fp
, src
->mPath
, order
, 2, & format
) ||
1594 !ReadBin4(fp
, src
->mPath
, order
, 2, & channels
) ||
1595 !ReadBin4(fp
, src
->mPath
, order
, 4, & rate
) ||
1596 !ReadBin4(fp
, src
->mPath
, order
, 4, & dummy
) ||
1597 !ReadBin4(fp
, src
->mPath
, order
, 2, & block
))
1602 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &size
))
1610 if(format
== WAVE_FORMAT_EXTENSIBLE
)
1612 fseek(fp
, 2, SEEK_CUR
);
1613 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &bits
))
1617 fseek(fp
, 4, SEEK_CUR
);
1618 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &format
))
1620 fseek(fp
, (long)(chunkSize
- 26), SEEK_CUR
);
1626 fseek(fp
, (long)(chunkSize
- 16), SEEK_CUR
);
1628 fseek(fp
, (long)(chunkSize
- 14), SEEK_CUR
);
1630 if(format
!= WAVE_FORMAT_PCM
&& format
!= WAVE_FORMAT_IEEE_FLOAT
)
1632 fprintf(stderr
, "Error: Unsupported WAVE format in file '%s'.\n", src
->mPath
);
1635 if(src
->mChannel
>= channels
)
1637 fprintf(stderr
, "Error: Missing source channel in WAVE file '%s'.\n", src
->mPath
);
1640 if(rate
!= hrirRate
)
1642 fprintf(stderr
, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src
->mPath
);
1645 if(format
== WAVE_FORMAT_PCM
)
1647 if(size
< 2 || size
> 4)
1649 fprintf(stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
->mPath
);
1652 if(bits
< 16 || bits
> (8*size
))
1654 fprintf (stderr
, "Error: Bad significant bits in WAVE file '%s'.\n", src
->mPath
);
1657 src
->mType
= ET_INT
;
1661 if(size
!= 4 && size
!= 8)
1663 fprintf(stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
->mPath
);
1669 src
->mBits
= (int)bits
;
1670 src
->mSkip
= channels
;
1674 // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
1675 static int ReadWaveData(FILE *fp
, const SourceRefT
*src
, const ByteOrderT order
, const uint n
, double *hrir
)
1677 int pre
, post
, skip
;
1680 pre
= (int)(src
->mSize
* src
->mChannel
);
1681 post
= (int)(src
->mSize
* (src
->mSkip
- src
->mChannel
- 1));
1683 for(i
= 0;i
< n
;i
++)
1687 fseek(fp
, skip
, SEEK_CUR
);
1688 if(!ReadBinAsDouble(fp
, src
->mPath
, order
, src
->mType
, src
->mSize
, src
->mBits
, &hrir
[i
]))
1693 fseek(fp
, skip
, SEEK_CUR
);
1697 // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
1699 static int ReadWaveList(FILE *fp
, const SourceRefT
*src
, const ByteOrderT order
, const uint n
, double *hrir
)
1701 uint32 fourCC
, chunkSize
, listSize
, count
;
1702 uint block
, skip
, offset
, i
;
1706 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, & fourCC
) ||
1707 !ReadBin4(fp
, src
->mPath
, order
, 4, & chunkSize
))
1710 if(fourCC
== FOURCC_DATA
)
1712 block
= src
->mSize
* src
->mSkip
;
1713 count
= chunkSize
/ block
;
1714 if(count
< (src
->mOffset
+ n
))
1716 fprintf(stderr
, "Error: Bad read from file '%s'.\n", src
->mPath
);
1719 fseek(fp
, (long)(src
->mOffset
* block
), SEEK_CUR
);
1720 if(!ReadWaveData(fp
, src
, order
, n
, &hrir
[0]))
1724 else if(fourCC
== FOURCC_LIST
)
1726 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
))
1729 if(fourCC
== FOURCC_WAVL
)
1733 fseek(fp
, (long)chunkSize
, SEEK_CUR
);
1735 listSize
= chunkSize
;
1736 block
= src
->mSize
* src
->mSkip
;
1737 skip
= src
->mOffset
;
1740 while(offset
< n
&& listSize
> 8)
1742 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1743 !ReadBin4(fp
, src
->mPath
, order
, 4, &chunkSize
))
1745 listSize
-= 8 + chunkSize
;
1746 if(fourCC
== FOURCC_DATA
)
1748 count
= chunkSize
/ block
;
1751 fseek(fp
, (long)(skip
* block
), SEEK_CUR
);
1752 chunkSize
-= skip
* block
;
1755 if(count
> (n
- offset
))
1757 if(!ReadWaveData(fp
, src
, order
, count
, &hrir
[offset
]))
1759 chunkSize
-= count
* block
;
1761 lastSample
= hrir
[offset
- 1];
1769 else if(fourCC
== FOURCC_SLNT
)
1771 if(!ReadBin4(fp
, src
->mPath
, order
, 4, &count
))
1778 if(count
> (n
- offset
))
1780 for(i
= 0; i
< count
; i
++)
1781 hrir
[offset
+ i
] = lastSample
;
1791 fseek(fp
, (long)chunkSize
, SEEK_CUR
);
1795 fprintf(stderr
, "Error: Bad read from file '%s'.\n", src
->mPath
);
1801 // Load a source HRIR from a RIFF/RIFX WAVE file.
1802 static int LoadWaveSource(FILE *fp
, SourceRefT
*src
, const uint hrirRate
, const uint n
, double *hrir
)
1804 uint32 fourCC
, dummy
;
1807 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1808 !ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &dummy
))
1810 if(fourCC
== FOURCC_RIFF
)
1812 else if(fourCC
== FOURCC_RIFX
)
1816 fprintf(stderr
, "Error: No RIFF/RIFX chunk in file '%s'.\n", src
->mPath
);
1820 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
))
1822 if(fourCC
!= FOURCC_WAVE
)
1824 fprintf(stderr
, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src
->mPath
);
1827 if(!ReadWaveFormat(fp
, order
, hrirRate
, src
))
1829 if(!ReadWaveList(fp
, src
, order
, n
, hrir
))
1834 // Load a source HRIR from a binary file.
1835 static int LoadBinarySource(FILE *fp
, const SourceRefT
*src
, const ByteOrderT order
, const uint n
, double *hrir
)
1839 fseek(fp
, (long)src
->mOffset
, SEEK_SET
);
1840 for(i
= 0;i
< n
;i
++)
1842 if(!ReadBinAsDouble(fp
, src
->mPath
, order
, src
->mType
, src
->mSize
, src
->mBits
, &hrir
[i
]))
1845 fseek(fp
, (long)src
->mSkip
, SEEK_CUR
);
1850 // Load a source HRIR from an ASCII text file containing a list of elements
1851 // separated by whitespace or common list operators (',', ';', ':', '|').
1852 static int LoadAsciiSource(FILE *fp
, const SourceRefT
*src
, const uint n
, double *hrir
)
1858 TrSetup(fp
, NULL
, &tr
);
1859 for(i
= 0;i
< src
->mOffset
;i
++)
1861 if(!ReadAsciiAsDouble(&tr
, src
->mPath
, src
->mType
, (uint
)src
->mBits
, &dummy
))
1864 for(i
= 0;i
< n
;i
++)
1866 if(!ReadAsciiAsDouble(&tr
, src
->mPath
, src
->mType
, (uint
)src
->mBits
, &hrir
[i
]))
1868 for(j
= 0;j
< src
->mSkip
;j
++)
1870 if(!ReadAsciiAsDouble(&tr
, src
->mPath
, src
->mType
, (uint
)src
->mBits
, &dummy
))
1877 // Load a source HRIR from a supported file type.
1878 static int LoadSource(SourceRefT
*src
, const uint hrirRate
, const uint n
, double *hrir
)
1883 if (src
->mFormat
== SF_ASCII
)
1884 fp
= fopen(src
->mPath
, "r");
1886 fp
= fopen(src
->mPath
, "rb");
1889 fprintf(stderr
, "Error: Could not open source file '%s'.\n", src
->mPath
);
1892 if(src
->mFormat
== SF_WAVE
)
1893 result
= LoadWaveSource(fp
, src
, hrirRate
, n
, hrir
);
1894 else if(src
->mFormat
== SF_BIN_LE
)
1895 result
= LoadBinarySource(fp
, src
, BO_LITTLE
, n
, hrir
);
1896 else if(src
->mFormat
== SF_BIN_BE
)
1897 result
= LoadBinarySource(fp
, src
, BO_BIG
, n
, hrir
);
1899 result
= LoadAsciiSource(fp
, src
, n
, hrir
);
1905 /***************************
1906 *** File storage output ***
1907 ***************************/
1909 // Write an ASCII string to a file.
1910 static int WriteAscii(const char *out
, FILE *fp
, const char *filename
)
1915 if(fwrite(out
, 1, len
, fp
) != len
)
1918 fprintf(stderr
, "Error: Bad write to file '%s'.\n", filename
);
1924 // Write a binary value of the given byte order and byte size to a file,
1925 // loading it from a 32-bit unsigned integer.
1926 static int WriteBin4(const ByteOrderT order
, const uint bytes
, const uint32 in
, FILE *fp
, const char *filename
)
1934 for(i
= 0;i
< bytes
;i
++)
1935 out
[i
] = (in
>>(i
*8)) & 0x000000FF;
1938 for(i
= 0;i
< bytes
;i
++)
1939 out
[bytes
- i
- 1] = (in
>>(i
*8)) & 0x000000FF;
1944 if(fwrite(out
, 1, bytes
, fp
) != bytes
)
1946 fprintf(stderr
, "Error: Bad write to file '%s'.\n", filename
);
1952 // Store the OpenAL Soft HRTF data set.
1953 static int StoreMhr(const HrirDataT
*hData
, const int experimental
, const char *filename
)
1955 uint e
, step
, end
, n
, j
, i
;
1960 if((fp
=fopen(filename
, "wb")) == NULL
)
1962 fprintf(stderr
, "Error: Could not open MHR file '%s'.\n", filename
);
1965 if(!WriteAscii(experimental
? MHR_FORMAT_EXPERIMENTAL
: MHR_FORMAT
, fp
, filename
))
1967 if(!WriteBin4(BO_LITTLE
, 4, (uint32
)hData
->mIrRate
, fp
, filename
))
1971 if(!WriteBin4(BO_LITTLE
, 1, (uint32
)hData
->mSampleType
, fp
, filename
))
1973 if(!WriteBin4(BO_LITTLE
, 1, (uint32
)hData
->mChannelType
, fp
, filename
))
1976 if(!WriteBin4(BO_LITTLE
, 1, (uint32
)hData
->mIrPoints
, fp
, filename
))
1978 if(!WriteBin4(BO_LITTLE
, 1, (uint32
)hData
->mEvCount
, fp
, filename
))
1980 for(e
= 0;e
< hData
->mEvCount
;e
++)
1982 if(!WriteBin4(BO_LITTLE
, 1, (uint32
)hData
->mAzCount
[e
], fp
, filename
))
1985 step
= hData
->mIrSize
;
1986 end
= hData
->mIrCount
* step
;
1987 n
= hData
->mIrPoints
;
1988 dither_seed
= 22222;
1989 for(j
= 0;j
< end
;j
+= step
)
1991 const double scale
= (!experimental
|| hData
->mSampleType
== ST_S16
) ? 32767.0 :
1992 ((hData
->mSampleType
== ST_S24
) ? 8388607.0 : 0.0);
1993 const int bps
= (!experimental
|| hData
->mSampleType
== ST_S16
) ? 2 :
1994 ((hData
->mSampleType
== ST_S24
) ? 3 : 0);
1995 double out
[MAX_TRUNCSIZE
];
1996 for(i
= 0;i
< n
;i
++)
1997 out
[i
] = TpdfDither(scale
* hData
->mHrirs
[j
+i
], &dither_seed
);
1998 for(i
= 0;i
< n
;i
++)
2000 v
= (int)Clamp(out
[i
], -scale
-1.0, scale
);
2001 if(!WriteBin4(BO_LITTLE
, bps
, (uint32
)v
, fp
, filename
))
2005 for(j
= 0;j
< hData
->mIrCount
;j
++)
2007 v
= (int)fmin(round(hData
->mIrRate
* hData
->mHrtds
[j
]), MAX_HRTD
);
2008 if(!WriteBin4(BO_LITTLE
, 1, (uint32
)v
, fp
, filename
))
2016 /***********************
2017 *** HRTF processing ***
2018 ***********************/
2020 // Calculate the onset time of an HRIR and average it with any existing
2021 // timing for its elevation and azimuth.
2022 static void AverageHrirOnset(const double *hrir
, const double f
, const uint ei
, const uint ai
, const HrirDataT
*hData
)
2028 n
= hData
->mIrPoints
;
2029 for(i
= 0;i
< n
;i
++)
2030 mag
= fmax(fabs(hrir
[i
]), mag
);
2032 for(i
= 0;i
< n
;i
++)
2034 if(fabs(hrir
[i
]) >= mag
)
2037 j
= hData
->mEvOffset
[ei
] + ai
;
2038 hData
->mHrtds
[j
] = Lerp(hData
->mHrtds
[j
], ((double)i
) / hData
->mIrRate
, f
);
2041 // Calculate the magnitude response of an HRIR and average it with any
2042 // existing responses for its elevation and azimuth.
2043 static void AverageHrirMagnitude(const double *hrir
, const double f
, const uint ei
, const uint ai
, const HrirDataT
*hData
)
2049 n
= hData
->mFftSize
;
2050 cplx
= calloc(sizeof(*cplx
), n
);
2051 mags
= calloc(sizeof(*mags
), n
);
2052 for(i
= 0;i
< hData
->mIrPoints
;i
++)
2053 cplx
[i
] = MakeComplex(hrir
[i
], 0.0);
2055 cplx
[i
] = MakeComplex(0.0, 0.0);
2056 FftForward(n
, cplx
, cplx
);
2057 MagnitudeResponse(n
, cplx
, mags
);
2059 j
= (hData
->mEvOffset
[ei
] + ai
) * hData
->mIrSize
;
2060 for(i
= 0;i
< m
;i
++)
2061 hData
->mHrirs
[j
+i
] = Lerp(hData
->mHrirs
[j
+i
], mags
[i
], f
);
2066 /* Calculate the contribution of each HRIR to the diffuse-field average based
2067 * on the area of its surface patch. All patches are centered at the HRIR
2068 * coordinates on the unit sphere and are measured by solid angle.
2070 static void CalculateDfWeights(const HrirDataT
*hData
, double *weights
)
2072 double evs
, sum
, ev
, up_ev
, down_ev
, solidAngle
;
2075 evs
= 90.0 / (hData
->mEvCount
- 1);
2077 for(ei
= hData
->mEvStart
;ei
< hData
->mEvCount
;ei
++)
2079 // For each elevation, calculate the upper and lower limits of the
2081 ev
= -90.0 + (ei
* 2.0 * evs
);
2082 if(ei
< (hData
->mEvCount
- 1))
2083 up_ev
= (ev
+ evs
) * M_PI
/ 180.0;
2087 down_ev
= (ev
- evs
) * M_PI
/ 180.0;
2089 down_ev
= -M_PI
/ 2.0;
2090 // Calculate the area of the patch band.
2091 solidAngle
= 2.0 * M_PI
* (sin(up_ev
) - sin(down_ev
));
2092 // Each weight is the area of one patch.
2093 weights
[ei
] = solidAngle
/ hData
->mAzCount
[ei
];
2094 // Sum the total surface area covered by the HRIRs.
2097 // Normalize the weights given the total surface coverage.
2098 for(ei
= hData
->mEvStart
;ei
< hData
->mEvCount
;ei
++)
2102 /* Calculate the diffuse-field average from the given magnitude responses of
2103 * the HRIR set. Weighting can be applied to compensate for the varying
2104 * surface area covered by each HRIR. The final average can then be limited
2105 * by the specified magnitude range (in positive dB; 0.0 to skip).
2107 static void CalculateDiffuseFieldAverage(const HrirDataT
*hData
, const int weighted
, const double limit
, double *dfa
)
2109 uint ei
, ai
, count
, step
, start
, end
, m
, j
, i
;
2112 weights
= CreateArray(hData
->mEvCount
);
2115 // Use coverage weighting to calculate the average.
2116 CalculateDfWeights(hData
, weights
);
2120 // If coverage weighting is not used, the weights still need to be
2121 // averaged by the number of HRIRs.
2123 for(ei
= hData
->mEvStart
;ei
< hData
->mEvCount
;ei
++)
2124 count
+= hData
->mAzCount
[ei
];
2125 for(ei
= hData
->mEvStart
;ei
< hData
->mEvCount
;ei
++)
2126 weights
[ei
] = 1.0 / count
;
2128 ei
= hData
->mEvStart
;
2130 step
= hData
->mIrSize
;
2131 start
= hData
->mEvOffset
[ei
] * step
;
2132 end
= hData
->mIrCount
* step
;
2133 m
= 1 + (hData
->mFftSize
/ 2);
2134 for(i
= 0;i
< m
;i
++)
2136 for(j
= start
;j
< end
;j
+= step
)
2138 // Get the weight for this HRIR's contribution.
2139 double weight
= weights
[ei
];
2140 // Add this HRIR's weighted power average to the total.
2141 for(i
= 0;i
< m
;i
++)
2142 dfa
[i
] += weight
* hData
->mHrirs
[j
+i
] * hData
->mHrirs
[j
+i
];
2143 // Determine the next weight to use.
2145 if(ai
>= hData
->mAzCount
[ei
])
2151 // Finish the average calculation and keep it from being too small.
2152 for(i
= 0;i
< m
;i
++)
2153 dfa
[i
] = fmax(sqrt(dfa
[i
]), EPSILON
);
2154 // Apply a limit to the magnitude range of the diffuse-field average if
2157 LimitMagnitudeResponse(hData
->mFftSize
, limit
, dfa
, dfa
);
2158 DestroyArray(weights
);
2161 // Perform diffuse-field equalization on the magnitude responses of the HRIR
2162 // set using the given average response.
2163 static void DiffuseFieldEqualize(const double *dfa
, const HrirDataT
*hData
)
2165 uint step
, start
, end
, m
, j
, i
;
2167 step
= hData
->mIrSize
;
2168 start
= hData
->mEvOffset
[hData
->mEvStart
] * step
;
2169 end
= hData
->mIrCount
* step
;
2170 m
= 1 + (hData
->mFftSize
/ 2);
2171 for(j
= start
;j
< end
;j
+= step
)
2173 for(i
= 0;i
< m
;i
++)
2174 hData
->mHrirs
[j
+i
] /= dfa
[i
];
2178 // Perform minimum-phase reconstruction using the magnitude responses of the
2180 static void ReconstructHrirs(const HrirDataT
*hData
)
2182 uint step
, start
, end
, n
, j
, i
;
2183 uint pcdone
, lastpc
;
2186 pcdone
= lastpc
= 0;
2187 printf("%3d%% done.", pcdone
);
2190 step
= hData
->mIrSize
;
2191 start
= hData
->mEvOffset
[hData
->mEvStart
] * step
;
2192 end
= hData
->mIrCount
* step
;
2193 n
= hData
->mFftSize
;
2194 cplx
= calloc(sizeof(*cplx
), n
);
2195 for(j
= start
;j
< end
;j
+= step
)
2197 MinimumPhase(n
, &hData
->mHrirs
[j
], cplx
);
2198 FftInverse(n
, cplx
, cplx
);
2199 for(i
= 0;i
< hData
->mIrPoints
;i
++)
2200 hData
->mHrirs
[j
+i
] = cplx
[i
].Real
;
2201 pcdone
= (j
+step
-start
) * 100 / (end
-start
);
2202 if(pcdone
!= lastpc
)
2205 printf("\r%3d%% done.", pcdone
);
2213 // Resamples the HRIRs for use at the given sampling rate.
2214 static void ResampleHrirs(const uint rate
, HrirDataT
*hData
)
2216 uint n
, step
, start
, end
, j
;
2219 ResamplerSetup(&rs
, hData
->mIrRate
, rate
);
2220 n
= hData
->mIrPoints
;
2221 step
= hData
->mIrSize
;
2222 start
= hData
->mEvOffset
[hData
->mEvStart
] * step
;
2223 end
= hData
->mIrCount
* step
;
2224 for(j
= start
;j
< end
;j
+= step
)
2225 ResamplerRun(&rs
, n
, &hData
->mHrirs
[j
], n
, &hData
->mHrirs
[j
]);
2226 ResamplerClear(&rs
);
2227 hData
->mIrRate
= rate
;
2230 /* Given an elevation index and an azimuth, calculate the indices of the two
2231 * HRIRs that bound the coordinate along with a factor for calculating the
2232 * continous HRIR using interpolation.
2234 static void CalcAzIndices(const HrirDataT
*hData
, const uint ei
, const double az
, uint
*j0
, uint
*j1
, double *jf
)
2239 af
= ((2.0*M_PI
) + az
) * hData
->mAzCount
[ei
] / (2.0*M_PI
);
2240 ai
= ((uint
)af
) % hData
->mAzCount
[ei
];
2243 *j0
= hData
->mEvOffset
[ei
] + ai
;
2244 *j1
= hData
->mEvOffset
[ei
] + ((ai
+1) % hData
->mAzCount
[ei
]);
2248 // Synthesize any missing onset timings at the bottom elevations. This just
2249 // blends between slightly exaggerated known onsets. Not an accurate model.
2250 static void SynthesizeOnsets(HrirDataT
*hData
)
2252 uint oi
, e
, a
, j0
, j1
;
2255 oi
= hData
->mEvStart
;
2257 for(a
= 0;a
< hData
->mAzCount
[oi
];a
++)
2258 t
+= hData
->mHrtds
[hData
->mEvOffset
[oi
] + a
];
2259 hData
->mHrtds
[0] = 1.32e-4 + (t
/ hData
->mAzCount
[oi
]);
2260 for(e
= 1;e
< hData
->mEvStart
;e
++)
2262 of
= ((double)e
) / hData
->mEvStart
;
2263 for(a
= 0;a
< hData
->mAzCount
[e
];a
++)
2265 CalcAzIndices(hData
, oi
, a
* 2.0 * M_PI
/ hData
->mAzCount
[e
], &j0
, &j1
, &jf
);
2266 hData
->mHrtds
[hData
->mEvOffset
[e
] + a
] = Lerp(hData
->mHrtds
[0], Lerp(hData
->mHrtds
[j0
], hData
->mHrtds
[j1
], jf
), of
);
2271 /* Attempt to synthesize any missing HRIRs at the bottom elevations. Right
2272 * now this just blends the lowest elevation HRIRs together and applies some
2273 * attenuation and high frequency damping. It is a simple, if inaccurate
2276 static void SynthesizeHrirs (HrirDataT
*hData
)
2278 uint oi
, a
, e
, step
, n
, i
, j
;
2279 double lp
[4], s0
, s1
;
2284 if(hData
->mEvStart
<= 0)
2286 step
= hData
->mIrSize
;
2287 oi
= hData
->mEvStart
;
2288 n
= hData
->mIrPoints
;
2289 for(i
= 0;i
< n
;i
++)
2290 hData
->mHrirs
[i
] = 0.0;
2291 for(a
= 0;a
< hData
->mAzCount
[oi
];a
++)
2293 j
= (hData
->mEvOffset
[oi
] + a
) * step
;
2294 for(i
= 0;i
< n
;i
++)
2295 hData
->mHrirs
[i
] += hData
->mHrirs
[j
+i
] / hData
->mAzCount
[oi
];
2297 for(e
= 1;e
< hData
->mEvStart
;e
++)
2299 of
= ((double)e
) / hData
->mEvStart
;
2300 b
= (1.0 - of
) * (3.5e-6 * hData
->mIrRate
);
2301 for(a
= 0;a
< hData
->mAzCount
[e
];a
++)
2303 j
= (hData
->mEvOffset
[e
] + a
) * step
;
2304 CalcAzIndices(hData
, oi
, a
* 2.0 * M_PI
/ hData
->mAzCount
[e
], &j0
, &j1
, &jf
);
2311 for(i
= 0;i
< n
;i
++)
2313 s0
= hData
->mHrirs
[i
];
2314 s1
= Lerp(hData
->mHrirs
[j0
+i
], hData
->mHrirs
[j1
+i
], jf
);
2315 s0
= Lerp(s0
, s1
, of
);
2316 lp
[0] = Lerp(s0
, lp
[0], b
);
2317 lp
[1] = Lerp(lp
[0], lp
[1], b
);
2318 lp
[2] = Lerp(lp
[1], lp
[2], b
);
2319 lp
[3] = Lerp(lp
[2], lp
[3], b
);
2320 hData
->mHrirs
[j
+i
] = lp
[3];
2324 b
= 3.5e-6 * hData
->mIrRate
;
2329 for(i
= 0;i
< n
;i
++)
2331 s0
= hData
->mHrirs
[i
];
2332 lp
[0] = Lerp(s0
, lp
[0], b
);
2333 lp
[1] = Lerp(lp
[0], lp
[1], b
);
2334 lp
[2] = Lerp(lp
[1], lp
[2], b
);
2335 lp
[3] = Lerp(lp
[2], lp
[3], b
);
2336 hData
->mHrirs
[i
] = lp
[3];
2338 hData
->mEvStart
= 0;
2341 // The following routines assume a full set of HRIRs for all elevations.
2343 // Normalize the HRIR set and slightly attenuate the result.
2344 static void NormalizeHrirs (const HrirDataT
*hData
)
2346 uint step
, end
, n
, j
, i
;
2349 step
= hData
->mIrSize
;
2350 end
= hData
->mIrCount
* step
;
2351 n
= hData
->mIrPoints
;
2353 for(j
= 0;j
< end
;j
+= step
)
2355 for(i
= 0;i
< n
;i
++)
2356 maxLevel
= fmax(fabs(hData
->mHrirs
[j
+i
]), maxLevel
);
2358 maxLevel
= 1.01 * maxLevel
;
2359 for(j
= 0;j
< end
;j
+= step
)
2361 for(i
= 0;i
< n
;i
++)
2362 hData
->mHrirs
[j
+i
] /= maxLevel
;
2366 // Calculate the left-ear time delay using a spherical head model.
2367 static double CalcLTD(const double ev
, const double az
, const double rad
, const double dist
)
2369 double azp
, dlp
, l
, al
;
2371 azp
= asin(cos(ev
) * sin(az
));
2372 dlp
= sqrt((dist
*dist
) + (rad
*rad
) + (2.0*dist
*rad
*sin(azp
)));
2373 l
= sqrt((dist
*dist
) - (rad
*rad
));
2374 al
= (0.5 * M_PI
) + azp
;
2376 dlp
= l
+ (rad
* (al
- acos(rad
/ dist
)));
2377 return (dlp
/ 343.3);
2380 // Calculate the effective head-related time delays for each minimum-phase
2382 static void CalculateHrtds (const HeadModelT model
, const double radius
, HrirDataT
*hData
)
2384 double minHrtd
, maxHrtd
;
2390 for(e
= 0;e
< hData
->mEvCount
;e
++)
2392 for(a
= 0;a
< hData
->mAzCount
[e
];a
++)
2394 j
= hData
->mEvOffset
[e
] + a
;
2395 if(model
== HM_DATASET
)
2396 t
= hData
->mHrtds
[j
] * radius
/ hData
->mRadius
;
2398 t
= CalcLTD((-90.0 + (e
* 180.0 / (hData
->mEvCount
- 1))) * M_PI
/ 180.0,
2399 (a
* 360.0 / hData
->mAzCount
[e
]) * M_PI
/ 180.0,
2400 radius
, hData
->mDistance
);
2401 hData
->mHrtds
[j
] = t
;
2402 maxHrtd
= fmax(t
, maxHrtd
);
2403 minHrtd
= fmin(t
, minHrtd
);
2407 for(j
= 0;j
< hData
->mIrCount
;j
++)
2408 hData
->mHrtds
[j
] -= minHrtd
;
2409 hData
->mMaxHrtd
= maxHrtd
;
2413 // Process the data set definition to read and validate the data set metrics.
2414 static int ProcessMetrics(TokenReaderT
*tr
, const uint fftSize
, const uint truncSize
, HrirDataT
*hData
)
2416 int hasRate
= 0, hasPoints
= 0, hasAzimuths
= 0;
2417 int hasRadius
= 0, hasDistance
= 0;
2418 char ident
[MAX_IDENT_LEN
+1];
2424 while(!(hasRate
&& hasPoints
&& hasAzimuths
&& hasRadius
&& hasDistance
))
2426 TrIndication(tr
, & line
, & col
);
2427 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2429 if(strcasecmp(ident
, "rate") == 0)
2433 TrErrorAt(tr
, line
, col
, "Redefinition of 'rate'.\n");
2436 if(!TrReadOperator(tr
, "="))
2438 if(!TrReadInt(tr
, MIN_RATE
, MAX_RATE
, &intVal
))
2440 hData
->mIrRate
= (uint
)intVal
;
2443 else if(strcasecmp(ident
, "points") == 0)
2446 TrErrorAt(tr
, line
, col
, "Redefinition of 'points'.\n");
2449 if(!TrReadOperator(tr
, "="))
2451 TrIndication(tr
, &line
, &col
);
2452 if(!TrReadInt(tr
, MIN_POINTS
, MAX_POINTS
, &intVal
))
2454 points
= (uint
)intVal
;
2455 if(fftSize
> 0 && points
> fftSize
)
2457 TrErrorAt(tr
, line
, col
, "Value exceeds the overridden FFT size.\n");
2460 if(points
< truncSize
)
2462 TrErrorAt(tr
, line
, col
, "Value is below the truncation size.\n");
2465 hData
->mIrPoints
= points
;
2468 hData
->mFftSize
= DEFAULT_FFTSIZE
;
2469 hData
->mIrSize
= 1 + (DEFAULT_FFTSIZE
/ 2);
2473 hData
->mFftSize
= fftSize
;
2474 hData
->mIrSize
= 1 + (fftSize
/ 2);
2475 if(points
> hData
->mIrSize
)
2476 hData
->mIrSize
= points
;
2480 else if(strcasecmp(ident
, "azimuths") == 0)
2484 TrErrorAt(tr
, line
, col
, "Redefinition of 'azimuths'.\n");
2487 if(!TrReadOperator(tr
, "="))
2489 hData
->mIrCount
= 0;
2490 hData
->mEvCount
= 0;
2491 hData
->mEvOffset
[0] = 0;
2494 if(!TrReadInt(tr
, MIN_AZ_COUNT
, MAX_AZ_COUNT
, &intVal
))
2496 hData
->mAzCount
[hData
->mEvCount
] = (uint
)intVal
;
2497 hData
->mIrCount
+= (uint
)intVal
;
2499 if(!TrIsOperator(tr
, ","))
2501 if(hData
->mEvCount
>= MAX_EV_COUNT
)
2503 TrError(tr
, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT
);
2506 hData
->mEvOffset
[hData
->mEvCount
] = hData
->mEvOffset
[hData
->mEvCount
- 1] + ((uint
)intVal
);
2507 TrReadOperator(tr
, ",");
2509 if(hData
->mEvCount
< MIN_EV_COUNT
)
2511 TrErrorAt(tr
, line
, col
, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT
);
2516 else if(strcasecmp(ident
, "radius") == 0)
2520 TrErrorAt(tr
, line
, col
, "Redefinition of 'radius'.\n");
2523 if(!TrReadOperator(tr
, "="))
2525 if(!TrReadFloat(tr
, MIN_RADIUS
, MAX_RADIUS
, &fpVal
))
2527 hData
->mRadius
= fpVal
;
2530 else if(strcasecmp(ident
, "distance") == 0)
2534 TrErrorAt(tr
, line
, col
, "Redefinition of 'distance'.\n");
2537 if(!TrReadOperator(tr
, "="))
2539 if(!TrReadFloat(tr
, MIN_DISTANCE
, MAX_DISTANCE
, & fpVal
))
2541 hData
->mDistance
= fpVal
;
2546 TrErrorAt(tr
, line
, col
, "Expected a metric name.\n");
2549 TrSkipWhitespace (tr
);
2554 // Parse an index pair from the data set definition.
2555 static int ReadIndexPair(TokenReaderT
*tr
, const HrirDataT
*hData
, uint
*ei
, uint
*ai
)
2558 if(!TrReadInt(tr
, 0, (int)hData
->mEvCount
, &intVal
))
2561 if(!TrReadOperator(tr
, ","))
2563 if(!TrReadInt(tr
, 0, (int)hData
->mAzCount
[*ei
], &intVal
))
2569 // Match the source format from a given identifier.
2570 static SourceFormatT
MatchSourceFormat(const char *ident
)
2572 if(strcasecmp(ident
, "wave") == 0)
2574 if(strcasecmp(ident
, "bin_le") == 0)
2576 if(strcasecmp(ident
, "bin_be") == 0)
2578 if(strcasecmp(ident
, "ascii") == 0)
2583 // Match the source element type from a given identifier.
2584 static ElementTypeT
MatchElementType(const char *ident
)
2586 if(strcasecmp(ident
, "int") == 0)
2588 if(strcasecmp(ident
, "fp") == 0)
2593 // Parse and validate a source reference from the data set definition.
2594 static int ReadSourceRef(TokenReaderT
*tr
, SourceRefT
*src
)
2596 char ident
[MAX_IDENT_LEN
+1];
2600 TrIndication(tr
, &line
, &col
);
2601 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2603 src
->mFormat
= MatchSourceFormat(ident
);
2604 if(src
->mFormat
== SF_NONE
)
2606 TrErrorAt(tr
, line
, col
, "Expected a source format.\n");
2609 if(!TrReadOperator(tr
, "("))
2611 if(src
->mFormat
== SF_WAVE
)
2613 if(!TrReadInt(tr
, 0, MAX_WAVE_CHANNELS
, &intVal
))
2615 src
->mType
= ET_NONE
;
2618 src
->mChannel
= (uint
)intVal
;
2623 TrIndication(tr
, &line
, &col
);
2624 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2626 src
->mType
= MatchElementType(ident
);
2627 if(src
->mType
== ET_NONE
)
2629 TrErrorAt(tr
, line
, col
, "Expected a source element type.\n");
2632 if(src
->mFormat
== SF_BIN_LE
|| src
->mFormat
== SF_BIN_BE
)
2634 if(!TrReadOperator(tr
, ","))
2636 if(src
->mType
== ET_INT
)
2638 if(!TrReadInt(tr
, MIN_BIN_SIZE
, MAX_BIN_SIZE
, &intVal
))
2640 src
->mSize
= (uint
)intVal
;
2641 if(!TrIsOperator(tr
, ","))
2642 src
->mBits
= (int)(8*src
->mSize
);
2645 TrReadOperator(tr
, ",");
2646 TrIndication(tr
, &line
, &col
);
2647 if(!TrReadInt(tr
, -2147483647-1, 2147483647, &intVal
))
2649 if(abs(intVal
) < MIN_BIN_BITS
|| ((uint
)abs(intVal
)) > (8*src
->mSize
))
2651 TrErrorAt(tr
, line
, col
, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS
, 8*src
->mSize
);
2654 src
->mBits
= intVal
;
2659 TrIndication(tr
, &line
, &col
);
2660 if(!TrReadInt(tr
, -2147483647-1, 2147483647, &intVal
))
2662 if(intVal
!= 4 && intVal
!= 8)
2664 TrErrorAt(tr
, line
, col
, "Expected a value of 4 or 8.\n");
2667 src
->mSize
= (uint
)intVal
;
2671 else if(src
->mFormat
== SF_ASCII
&& src
->mType
== ET_INT
)
2673 if(!TrReadOperator(tr
, ","))
2675 if(!TrReadInt(tr
, MIN_ASCII_BITS
, MAX_ASCII_BITS
, &intVal
))
2678 src
->mBits
= intVal
;
2686 if(!TrIsOperator(tr
, ";"))
2690 TrReadOperator(tr
, ";");
2691 if(!TrReadInt (tr
, 0, 0x7FFFFFFF, &intVal
))
2693 src
->mSkip
= (uint
)intVal
;
2696 if(!TrReadOperator(tr
, ")"))
2698 if(TrIsOperator(tr
, "@"))
2700 TrReadOperator(tr
, "@");
2701 if(!TrReadInt(tr
, 0, 0x7FFFFFFF, &intVal
))
2703 src
->mOffset
= (uint
)intVal
;
2707 if(!TrReadOperator(tr
, ":"))
2709 if(!TrReadString(tr
, MAX_PATH_LEN
, src
->mPath
))
2714 // Process the list of sources in the data set definition.
2715 static int ProcessSources(const HeadModelT model
, TokenReaderT
*tr
, HrirDataT
*hData
)
2717 uint
*setCount
, *setFlag
;
2718 uint line
, col
, ei
, ai
;
2724 printf("Loading sources...");
2728 setCount
= (uint
*)calloc(hData
->mEvCount
, sizeof(uint
));
2729 setFlag
= (uint
*)calloc(hData
->mIrCount
, sizeof(uint
));
2730 hrir
= CreateArray(hData
->mIrPoints
);
2731 while(TrIsOperator(tr
, "["))
2733 TrIndication(tr
, & line
, & col
);
2734 TrReadOperator(tr
, "[");
2735 if(!ReadIndexPair(tr
, hData
, &ei
, &ai
))
2737 if(!TrReadOperator(tr
, "]"))
2739 if(setFlag
[hData
->mEvOffset
[ei
] + ai
])
2741 TrErrorAt(tr
, line
, col
, "Redefinition of source.\n");
2744 if(!TrReadOperator(tr
, "="))
2750 if(!ReadSourceRef(tr
, &src
))
2753 // TODO: Would be nice to display 'x of y files', but that would
2754 // require preparing the source refs first to get a total count
2755 // before loading them.
2757 printf("\rLoading sources... %d file%s", count
, (count
==1)?"":"s");
2760 if(!LoadSource(&src
, hData
->mIrRate
, hData
->mIrPoints
, hrir
))
2763 if(model
== HM_DATASET
)
2764 AverageHrirOnset(hrir
, 1.0 / factor
, ei
, ai
, hData
);
2765 AverageHrirMagnitude(hrir
, 1.0 / factor
, ei
, ai
, hData
);
2767 if(!TrIsOperator(tr
, "+"))
2769 TrReadOperator(tr
, "+");
2771 setFlag
[hData
->mEvOffset
[ei
] + ai
] = 1;
2777 while(ei
< hData
->mEvCount
&& setCount
[ei
] < 1)
2779 if(ei
< hData
->mEvCount
)
2781 hData
->mEvStart
= ei
;
2782 while(ei
< hData
->mEvCount
&& setCount
[ei
] == hData
->mAzCount
[ei
])
2784 if(ei
>= hData
->mEvCount
)
2793 TrError(tr
, "Errant data at end of source list.\n");
2796 TrError(tr
, "Missing sources for elevation index %d.\n", ei
);
2799 TrError(tr
, "Missing source references.\n");
2808 /* Parse the data set definition and process the source data, storing the
2809 * resulting data set as desired. If the input name is NULL it will read
2810 * from standard input.
2812 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 int experimental
, const char *outName
)
2814 char rateStr
[8+1], expName
[MAX_PATH_LEN
];
2821 hData
.mSampleType
= ST_S24
;
2822 hData
.mChannelType
= CT_LEFTONLY
;
2823 hData
.mIrPoints
= 0;
2829 hData
.mDistance
= 0;
2830 fprintf(stdout
, "Reading HRIR definition...\n");
2833 fp
= fopen(inName
, "r");
2836 fprintf(stderr
, "Error: Could not open definition file '%s'\n", inName
);
2839 TrSetup(fp
, inName
, &tr
);
2844 TrSetup(fp
, "<stdin>", &tr
);
2846 if(!ProcessMetrics(&tr
, fftSize
, truncSize
, &hData
))
2852 hData
.mHrirs
= CreateArray(hData
.mIrCount
* hData
.mIrSize
);
2853 hData
.mHrtds
= CreateArray(hData
.mIrCount
);
2854 if(!ProcessSources(model
, &tr
, &hData
))
2856 DestroyArray(hData
.mHrtds
);
2857 DestroyArray(hData
.mHrirs
);
2866 dfa
= CreateArray(1 + (hData
.mFftSize
/2));
2867 fprintf(stdout
, "Calculating diffuse-field average...\n");
2868 CalculateDiffuseFieldAverage(&hData
, surface
, limit
, dfa
);
2869 fprintf(stdout
, "Performing diffuse-field equalization...\n");
2870 DiffuseFieldEqualize(dfa
, &hData
);
2873 fprintf(stdout
, "Performing minimum phase reconstruction...\n");
2874 ReconstructHrirs(&hData
);
2875 if(outRate
!= 0 && outRate
!= hData
.mIrRate
)
2877 fprintf(stdout
, "Resampling HRIRs...\n");
2878 ResampleHrirs(outRate
, &hData
);
2880 fprintf(stdout
, "Truncating minimum-phase HRIRs...\n");
2881 hData
.mIrPoints
= truncSize
;
2882 fprintf(stdout
, "Synthesizing missing elevations...\n");
2883 if(model
== HM_DATASET
)
2884 SynthesizeOnsets(&hData
);
2885 SynthesizeHrirs(&hData
);
2886 fprintf(stdout
, "Normalizing final HRIRs...\n");
2887 NormalizeHrirs(&hData
);
2888 fprintf(stdout
, "Calculating impulse delays...\n");
2889 CalculateHrtds(model
, (radius
> DEFAULT_CUSTOM_RADIUS
) ? radius
: hData
.mRadius
, &hData
);
2890 snprintf(rateStr
, 8, "%u", hData
.mIrRate
);
2891 StrSubst(outName
, "%r", rateStr
, MAX_PATH_LEN
, expName
);
2895 fprintf(stdout
, "Creating MHR data set file...\n");
2896 if(!StoreMhr(&hData
, experimental
, expName
))
2898 DestroyArray(hData
.mHrtds
);
2899 DestroyArray(hData
.mHrirs
);
2906 DestroyArray(hData
.mHrtds
);
2907 DestroyArray(hData
.mHrirs
);
2911 static void PrintHelp(const char *argv0
, FILE *ofile
)
2913 fprintf(ofile
, "Usage: %s <command> [<option>...]\n\n", argv0
);
2914 fprintf(ofile
, "Commands:\n");
2915 fprintf(ofile
, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n");
2916 fprintf(ofile
, " Defaults output to: ./oalsoft_hrtf_%%r.mhr\n");
2917 fprintf(ofile
, " -h, --help Displays this help information.\n\n");
2918 fprintf(ofile
, "Options:\n");
2919 fprintf(ofile
, " -r <rate> Change the data set sample rate to the specified value and\n");
2920 fprintf(ofile
, " resample the HRIRs accordingly.\n");
2921 fprintf(ofile
, " -f <points> Override the FFT window size (default: %u).\n", DEFAULT_FFTSIZE
);
2922 fprintf(ofile
, " -e {on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE
? "on" : "off"));
2923 fprintf(ofile
, " -s {on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE
? "on" : "off"));
2924 fprintf(ofile
, " -l {<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
2925 fprintf(ofile
, " average (default: %.2f).\n", DEFAULT_LIMIT
);
2926 fprintf(ofile
, " -w <points> Specify the size of the truncation window that's applied\n");
2927 fprintf(ofile
, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE
);
2928 fprintf(ofile
, " -d {dataset| Specify the model used for calculating the head-delay timing\n");
2929 fprintf(ofile
, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL
== HM_DATASET
) ? "dataset" : "sphere"));
2930 fprintf(ofile
, " -c <size> Use a customized head radius measured ear-to-ear in meters.\n");
2931 fprintf(ofile
, " -i <filename> Specify an HRIR definition file to use (defaults to stdin).\n");
2932 fprintf(ofile
, " -o <filename> Specify an output file. Overrides command-selected default.\n");
2933 fprintf(ofile
, " Use of '%%r' will be substituted with the data set sample rate.\n");
2937 #define main my_main
2938 int main(int argc
, char *argv
[]);
2940 static char **arglist
;
2941 static void cleanup_arglist(void)
2944 for(i
= 0;arglist
[i
];i
++)
2949 int wmain(int argc
, const wchar_t *wargv
[])
2953 atexit(cleanup_arglist
);
2954 arglist
= calloc(sizeof(*arglist
), argc
);
2955 for(i
= 0;i
< argc
;i
++)
2956 arglist
[i
] = ToUTF8(wargv
[i
]);
2958 return main(argc
, arglist
);
2962 // Standard command line dispatch.
2963 int main(int argc
, char *argv
[])
2965 const char *inName
= NULL
, *outName
= NULL
;
2966 OutputFormatT outFormat
;
2967 uint outRate
, fftSize
;
2968 int equalize
, surface
;
2977 if(argc
< 2 || strcmp(argv
[1], "--help") == 0 || strcmp(argv
[1], "-h") == 0)
2979 fprintf(stdout
, "HRTF Processing and Composition Utility\n\n");
2980 PrintHelp(argv
[0], stdout
);
2984 if(strcmp(argv
[1], "--make-mhr") == 0 || strcmp(argv
[1], "-m") == 0)
2986 outName
= "./oalsoft_hrtf_%r.mhr";
2991 fprintf(stderr
, "Error: Invalid command '%s'.\n\n", argv
[1]);
2992 PrintHelp(argv
[0], stderr
);
2998 equalize
= DEFAULT_EQUALIZE
;
2999 surface
= DEFAULT_SURFACE
;
3000 limit
= DEFAULT_LIMIT
;
3001 truncSize
= DEFAULT_TRUNCSIZE
;
3002 model
= DEFAULT_HEAD_MODEL
;
3003 radius
= DEFAULT_CUSTOM_RADIUS
;
3007 while((opt
=getopt(argc
, argv
, "r:f:e:s:k:w:d:c:e:i:o:xh")) != -1)
3012 outRate
= strtoul(optarg
, &end
, 10);
3013 if(end
[0] != '\0' || outRate
< MIN_RATE
|| outRate
> MAX_RATE
)
3015 fprintf(stderr
, "Error: Expected a value from %u to %u for '-r'.\n", MIN_RATE
, MAX_RATE
);
3021 fftSize
= strtoul(optarg
, &end
, 10);
3022 if(end
[0] != '\0' || (fftSize
&(fftSize
-1)) || fftSize
< MIN_FFTSIZE
|| fftSize
> MAX_FFTSIZE
)
3024 fprintf(stderr
, "Error: Expected a power-of-two value from %u to %u for '-f'.\n", MIN_FFTSIZE
, MAX_FFTSIZE
);
3030 if(strcmp(optarg
, "on") == 0)
3032 else if(strcmp(optarg
, "off") == 0)
3036 fprintf(stderr
, "Error: Expected 'on' or 'off' for '-e'.\n");
3042 if(strcmp(optarg
, "on") == 0)
3044 else if(strcmp(optarg
, "off") == 0)
3048 fprintf(stderr
, "Error: Expected 'on' or 'off' for '-s'.\n");
3054 if(strcmp(optarg
, "none") == 0)
3058 limit
= strtod(optarg
, &end
);
3059 if(end
[0] != '\0' || limit
< MIN_LIMIT
|| limit
> MAX_LIMIT
)
3061 fprintf(stderr
, "Error: Expected 'none' or a value from %.2f to %.2f for '-l'.\n", MIN_LIMIT
, MAX_LIMIT
);
3068 truncSize
= strtoul(optarg
, &end
, 10);
3069 if(end
[0] != '\0' || truncSize
< MIN_TRUNCSIZE
|| truncSize
> MAX_TRUNCSIZE
|| (truncSize
%MOD_TRUNCSIZE
))
3071 fprintf(stderr
, "Error: Expected a value from %u to %u in multiples of %u for '-w'.\n", MIN_TRUNCSIZE
, MAX_TRUNCSIZE
, MOD_TRUNCSIZE
);
3077 if(strcmp(optarg
, "dataset") == 0)
3079 else if(strcmp(optarg
, "sphere") == 0)
3083 fprintf(stderr
, "Error: Expected 'dataset' or 'sphere' for '-d'.\n");
3089 radius
= strtod(optarg
, &end
);
3090 if(end
[0] != '\0' || radius
< MIN_CUSTOM_RADIUS
|| radius
> MAX_CUSTOM_RADIUS
)
3092 fprintf(stderr
, "Error: Expected a value from %.2f to %.2f for '-c'.\n", MIN_CUSTOM_RADIUS
, MAX_CUSTOM_RADIUS
);
3110 PrintHelp(argv
[0], stdout
);
3114 PrintHelp(argv
[0], stderr
);
3119 if(!ProcessDefinition(inName
, outRate
, fftSize
, equalize
, surface
, limit
,
3120 truncSize
, model
, radius
, outFormat
, experimental
,
3123 fprintf(stdout
, "Operation completed.\n");
3125 return EXIT_SUCCESS
;