Avoid including AL headers in makehrtf
[openal-soft.git] / utils / makehrtf.c
blob516318e8620be7af9215a711fab07996eb67f0e3
1 /*
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
36 * measurement.
37 * 3. Take the diffuse-field average and limit its magnitude range.
38 * 4. Equalize the responses by using the inverse of the diffuse-field
39 * average.
40 * 5. Reconstruct the minimum-phase responses.
41 * 5. Zero the DC component.
42 * 6. IFFT the result and truncate to the desired-length minimum-phase FIR.
44 * The spherical head algorithm for calculating propagation delay was adapted
45 * from the paper:
47 * Modeling Interaural Time Difference Assuming a Spherical Head
48 * Joel David Miller
49 * Music 150, Musical Acoustics, Stanford University
50 * December 2, 2001
52 * The formulae for calculating the Kaiser window metrics are from the
53 * the textbook:
55 * Discrete-Time Signal Processing
56 * Alan V. Oppenheim and Ronald W. Schafer
57 * Prentice-Hall Signal Processing Series
58 * 1999
61 #include "config.h"
63 #define _UNICODE
64 #include <stdio.h>
65 #include <stdlib.h>
66 #include <stdarg.h>
67 #include <stddef.h>
68 #include <string.h>
69 #include <limits.h>
70 #include <ctype.h>
71 #include <math.h>
72 #ifdef HAVE_STRINGS_H
73 #include <strings.h>
74 #endif
75 #ifdef HAVE_GETOPT
76 #include <unistd.h>
77 #else
78 #include "getopt.h"
79 #endif
81 #include "win_main_utf8.h"
83 /* Define int64_t and uint64_t types */
84 #if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 199901L
85 #include <inttypes.h>
86 #elif defined(_WIN32) && defined(__GNUC__)
87 #include <stdint.h>
88 #elif defined(_WIN32)
89 typedef __int64 int64_t;
90 typedef unsigned __int64 uint64_t;
91 #else
92 /* Fallback if nothing above works */
93 #include <inttypes.h>
94 #endif
96 #ifndef M_PI
97 #define M_PI (3.14159265358979323846)
98 #endif
100 #ifndef HUGE_VAL
101 #define HUGE_VAL (1.0 / 0.0)
102 #endif
105 // The epsilon used to maintain signal stability.
106 #define EPSILON (1e-9)
108 // Constants for accessing the token reader's ring buffer.
109 #define TR_RING_BITS (16)
110 #define TR_RING_SIZE (1 << TR_RING_BITS)
111 #define TR_RING_MASK (TR_RING_SIZE - 1)
113 // The token reader's load interval in bytes.
114 #define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
116 // The maximum identifier length used when processing the data set
117 // definition.
118 #define MAX_IDENT_LEN (16)
120 // The maximum path length used when processing filenames.
121 #define MAX_PATH_LEN (256)
123 // The limits for the sample 'rate' metric in the data set definition and for
124 // resampling.
125 #define MIN_RATE (32000)
126 #define MAX_RATE (96000)
128 // The limits for the HRIR 'points' metric in the data set definition.
129 #define MIN_POINTS (16)
130 #define MAX_POINTS (8192)
132 // The limits to the number of 'azimuths' listed in the data set definition.
133 #define MIN_EV_COUNT (5)
134 #define MAX_EV_COUNT (128)
136 // The limits for each of the 'azimuths' listed in the data set definition.
137 #define MIN_AZ_COUNT (1)
138 #define MAX_AZ_COUNT (128)
140 // The limits for the listener's head 'radius' in the data set definition.
141 #define MIN_RADIUS (0.05)
142 #define MAX_RADIUS (0.15)
144 // The limits for the 'distance' from source to listener in the definition
145 // file.
146 #define MIN_DISTANCE (0.5)
147 #define MAX_DISTANCE (2.5)
149 // The maximum number of channels that can be addressed for a WAVE file
150 // source listed in the data set definition.
151 #define MAX_WAVE_CHANNELS (65535)
153 // The limits to the byte size for a binary source listed in the definition
154 // file.
155 #define MIN_BIN_SIZE (2)
156 #define MAX_BIN_SIZE (4)
158 // The minimum number of significant bits for binary sources listed in the
159 // data set definition. The maximum is calculated from the byte size.
160 #define MIN_BIN_BITS (16)
162 // The limits to the number of significant bits for an ASCII source listed in
163 // the data set definition.
164 #define MIN_ASCII_BITS (16)
165 #define MAX_ASCII_BITS (32)
167 // The limits to the FFT window size override on the command line.
168 #define MIN_FFTSIZE (65536)
169 #define MAX_FFTSIZE (131072)
171 // The limits to the equalization range limit on the command line.
172 #define MIN_LIMIT (2.0)
173 #define MAX_LIMIT (120.0)
175 // The limits to the truncation window size on the command line.
176 #define MIN_TRUNCSIZE (16)
177 #define MAX_TRUNCSIZE (512)
179 // The limits to the custom head radius on the command line.
180 #define MIN_CUSTOM_RADIUS (0.05)
181 #define MAX_CUSTOM_RADIUS (0.15)
183 // The truncation window size must be a multiple of the below value to allow
184 // for vectorized convolution.
185 #define MOD_TRUNCSIZE (8)
187 // The defaults for the command line options.
188 #define DEFAULT_FFTSIZE (65536)
189 #define DEFAULT_EQUALIZE (1)
190 #define DEFAULT_SURFACE (1)
191 #define DEFAULT_LIMIT (24.0)
192 #define DEFAULT_TRUNCSIZE (32)
193 #define DEFAULT_HEAD_MODEL (HM_DATASET)
194 #define DEFAULT_CUSTOM_RADIUS (0.0)
196 // The four-character-codes for RIFF/RIFX WAVE file chunks.
197 #define FOURCC_RIFF (0x46464952) // 'RIFF'
198 #define FOURCC_RIFX (0x58464952) // 'RIFX'
199 #define FOURCC_WAVE (0x45564157) // 'WAVE'
200 #define FOURCC_FMT (0x20746D66) // 'fmt '
201 #define FOURCC_DATA (0x61746164) // 'data'
202 #define FOURCC_LIST (0x5453494C) // 'LIST'
203 #define FOURCC_WAVL (0x6C766177) // 'wavl'
204 #define FOURCC_SLNT (0x746E6C73) // 'slnt'
206 // The supported wave formats.
207 #define WAVE_FORMAT_PCM (0x0001)
208 #define WAVE_FORMAT_IEEE_FLOAT (0x0003)
209 #define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
211 // The maximum propagation delay value supported by OpenAL Soft.
212 #define MAX_HRTD (63.0)
214 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
215 // response protocol 01.
216 #define MHR_FORMAT ("MinPHR01")
218 #define MHR_FORMAT_EXPERIMENTAL ("MinPHRTEMPDONOTUSE")
220 // Sample and channel type enum values
221 typedef enum SampleTypeT {
222 ST_S16 = 0,
223 ST_S24 = 1
224 } SampleTypeT;
226 typedef enum ChannelTypeT {
227 CT_LEFTONLY = 0,
228 CT_LEFTRIGHT = 1
229 } ChannelTypeT;
231 // Byte order for the serialization routines.
232 typedef enum ByteOrderT {
233 BO_NONE,
234 BO_LITTLE,
235 BO_BIG
236 } ByteOrderT;
238 // Source format for the references listed in the data set definition.
239 typedef enum SourceFormatT {
240 SF_NONE,
241 SF_WAVE, // RIFF/RIFX WAVE file.
242 SF_BIN_LE, // Little-endian binary file.
243 SF_BIN_BE, // Big-endian binary file.
244 SF_ASCII // ASCII text file.
245 } SourceFormatT;
247 // Element types for the references listed in the data set definition.
248 typedef enum ElementTypeT {
249 ET_NONE,
250 ET_INT, // Integer elements.
251 ET_FP // Floating-point elements.
252 } ElementTypeT;
254 // Head model used for calculating the impulse delays.
255 typedef enum HeadModelT {
256 HM_NONE,
257 HM_DATASET, // Measure the onset from the dataset.
258 HM_SPHERE // Calculate the onset using a spherical head model.
259 } HeadModelT;
261 // Unsigned integer type.
262 typedef unsigned int uint;
264 // Serialization types. The trailing digit indicates the number of bits.
265 typedef unsigned char uint8;
266 typedef int int32;
267 typedef unsigned int uint32;
268 typedef uint64_t uint64;
270 // Token reader state for parsing the data set definition.
271 typedef struct TokenReaderT {
272 FILE *mFile;
273 const char *mName;
274 uint mLine;
275 uint mColumn;
276 char mRing[TR_RING_SIZE];
277 size_t mIn;
278 size_t mOut;
279 } TokenReaderT;
281 // Source reference state used when loading sources.
282 typedef struct SourceRefT {
283 SourceFormatT mFormat;
284 ElementTypeT mType;
285 uint mSize;
286 int mBits;
287 uint mChannel;
288 uint mSkip;
289 uint mOffset;
290 char mPath[MAX_PATH_LEN+1];
291 } SourceRefT;
293 // The HRIR metrics and data set used when loading, processing, and storing
294 // the resulting HRTF.
295 typedef struct HrirDataT {
296 uint mIrRate;
297 SampleTypeT mSampleType;
298 ChannelTypeT mChannelType;
299 uint mIrCount;
300 uint mIrSize;
301 uint mIrPoints;
302 uint mFftSize;
303 uint mEvCount;
304 uint mEvStart;
305 uint mAzCount[MAX_EV_COUNT];
306 uint mEvOffset[MAX_EV_COUNT];
307 double mRadius;
308 double mDistance;
309 double *mHrirs;
310 double *mHrtds;
311 double mMaxHrtd;
312 } HrirDataT;
314 // The resampler metrics and FIR filter.
315 typedef struct ResamplerT {
316 uint mP, mQ, mM, mL;
317 double *mF;
318 } ResamplerT;
321 /****************************************
322 *** Complex number type and routines ***
323 ****************************************/
325 typedef struct {
326 double Real, Imag;
327 } Complex;
329 static Complex MakeComplex(double r, double i)
331 Complex c = { r, i };
332 return c;
335 static Complex c_add(Complex a, Complex b)
337 Complex r;
338 r.Real = a.Real + b.Real;
339 r.Imag = a.Imag + b.Imag;
340 return r;
343 static Complex c_sub(Complex a, Complex b)
345 Complex r;
346 r.Real = a.Real - b.Real;
347 r.Imag = a.Imag - b.Imag;
348 return r;
351 static Complex c_mul(Complex a, Complex b)
353 Complex r;
354 r.Real = a.Real*b.Real - a.Imag*b.Imag;
355 r.Imag = a.Imag*b.Real + a.Real*b.Imag;
356 return r;
359 static Complex c_muls(Complex a, double s)
361 Complex r;
362 r.Real = a.Real * s;
363 r.Imag = a.Imag * s;
364 return r;
367 static double c_abs(Complex a)
369 return sqrt(a.Real*a.Real + a.Imag*a.Imag);
372 static Complex c_exp(Complex a)
374 Complex r;
375 double e = exp(a.Real);
376 r.Real = e * cos(a.Imag);
377 r.Imag = e * sin(a.Imag);
378 return r;
381 /*****************************
382 *** Token reader routines ***
383 *****************************/
385 /* Whitespace is not significant. It can process tokens as identifiers, numbers
386 * (integer and floating-point), strings, and operators. Strings must be
387 * encapsulated by double-quotes and cannot span multiple lines.
390 // Setup the reader on the given file. The filename can be NULL if no error
391 // output is desired.
392 static void TrSetup(FILE *fp, const char *filename, TokenReaderT *tr)
394 const char *name = NULL;
396 if(filename)
398 const char *slash = strrchr(filename, '/');
399 if(slash)
401 const char *bslash = strrchr(slash+1, '\\');
402 if(bslash) name = bslash+1;
403 else name = slash+1;
405 else
407 const char *bslash = strrchr(filename, '\\');
408 if(bslash) name = bslash+1;
409 else name = filename;
413 tr->mFile = fp;
414 tr->mName = name;
415 tr->mLine = 1;
416 tr->mColumn = 1;
417 tr->mIn = 0;
418 tr->mOut = 0;
421 // Prime the reader's ring buffer, and return a result indicating that there
422 // is text to process.
423 static int TrLoad(TokenReaderT *tr)
425 size_t toLoad, in, count;
427 toLoad = TR_RING_SIZE - (tr->mIn - tr->mOut);
428 if(toLoad >= TR_LOAD_SIZE && !feof(tr->mFile))
430 // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
431 toLoad = TR_LOAD_SIZE;
432 in = tr->mIn&TR_RING_MASK;
433 count = TR_RING_SIZE - in;
434 if(count < toLoad)
436 tr->mIn += fread(&tr->mRing[in], 1, count, tr->mFile);
437 tr->mIn += fread(&tr->mRing[0], 1, toLoad-count, tr->mFile);
439 else
440 tr->mIn += fread(&tr->mRing[in], 1, toLoad, tr->mFile);
442 if(tr->mOut >= TR_RING_SIZE)
444 tr->mOut -= TR_RING_SIZE;
445 tr->mIn -= TR_RING_SIZE;
448 if(tr->mIn > tr->mOut)
449 return 1;
450 return 0;
453 // Error display routine. Only displays when the base name is not NULL.
454 static void TrErrorVA(const TokenReaderT *tr, uint line, uint column, const char *format, va_list argPtr)
456 if(!tr->mName)
457 return;
458 fprintf(stderr, "Error (%s:%u:%u): ", tr->mName, line, column);
459 vfprintf(stderr, format, argPtr);
462 // Used to display an error at a saved line/column.
463 static void TrErrorAt(const TokenReaderT *tr, uint line, uint column, const char *format, ...)
465 va_list argPtr;
467 va_start(argPtr, format);
468 TrErrorVA(tr, line, column, format, argPtr);
469 va_end(argPtr);
472 // Used to display an error at the current line/column.
473 static void TrError(const TokenReaderT *tr, const char *format, ...)
475 va_list argPtr;
477 va_start(argPtr, format);
478 TrErrorVA(tr, tr->mLine, tr->mColumn, format, argPtr);
479 va_end(argPtr);
482 // Skips to the next line.
483 static void TrSkipLine(TokenReaderT *tr)
485 char ch;
487 while(TrLoad(tr))
489 ch = tr->mRing[tr->mOut&TR_RING_MASK];
490 tr->mOut++;
491 if(ch == '\n')
493 tr->mLine++;
494 tr->mColumn = 1;
495 break;
497 tr->mColumn ++;
501 // Skips to the next token.
502 static int TrSkipWhitespace(TokenReaderT *tr)
504 char ch;
506 while(TrLoad(tr))
508 ch = tr->mRing[tr->mOut&TR_RING_MASK];
509 if(isspace(ch))
511 tr->mOut++;
512 if(ch == '\n')
514 tr->mLine++;
515 tr->mColumn = 1;
517 else
518 tr->mColumn++;
520 else if(ch == '#')
521 TrSkipLine(tr);
522 else
523 return 1;
525 return 0;
528 // Get the line and/or column of the next token (or the end of input).
529 static void TrIndication(TokenReaderT *tr, uint *line, uint *column)
531 TrSkipWhitespace(tr);
532 if(line) *line = tr->mLine;
533 if(column) *column = tr->mColumn;
536 // Checks to see if a token is the given operator. It does not display any
537 // errors and will not proceed to the next token.
538 static int TrIsOperator(TokenReaderT *tr, const char *op)
540 size_t out, len;
541 char ch;
543 if(!TrSkipWhitespace(tr))
544 return 0;
545 out = tr->mOut;
546 len = 0;
547 while(op[len] != '\0' && out < tr->mIn)
549 ch = tr->mRing[out&TR_RING_MASK];
550 if(ch != op[len]) break;
551 len++;
552 out++;
554 if(op[len] == '\0')
555 return 1;
556 return 0;
559 /* The TrRead*() routines obtain the value of a matching token type. They
560 * display type, form, and boundary errors and will proceed to the next
561 * token.
564 // Reads and validates an identifier token.
565 static int TrReadIdent(TokenReaderT *tr, const uint maxLen, char *ident)
567 uint col, len;
568 char ch;
570 col = tr->mColumn;
571 if(TrSkipWhitespace(tr))
573 col = tr->mColumn;
574 ch = tr->mRing[tr->mOut&TR_RING_MASK];
575 if(ch == '_' || isalpha(ch))
577 len = 0;
578 do {
579 if(len < maxLen)
580 ident[len] = ch;
581 len++;
582 tr->mOut++;
583 if(!TrLoad(tr))
584 break;
585 ch = tr->mRing[tr->mOut&TR_RING_MASK];
586 } while(ch == '_' || isdigit(ch) || isalpha(ch));
588 tr->mColumn += len;
589 if(len < maxLen)
591 ident[len] = '\0';
592 return 1;
594 TrErrorAt(tr, tr->mLine, col, "Identifier is too long.\n");
595 return 0;
598 TrErrorAt(tr, tr->mLine, col, "Expected an identifier.\n");
599 return 0;
602 // Reads and validates (including bounds) an integer token.
603 static int TrReadInt(TokenReaderT *tr, const int loBound, const int hiBound, int *value)
605 uint col, digis, len;
606 char ch, temp[64+1];
608 col = tr->mColumn;
609 if(TrSkipWhitespace(tr))
611 col = tr->mColumn;
612 len = 0;
613 ch = tr->mRing[tr->mOut&TR_RING_MASK];
614 if(ch == '+' || ch == '-')
616 temp[len] = ch;
617 len++;
618 tr->mOut++;
620 digis = 0;
621 while(TrLoad(tr))
623 ch = tr->mRing[tr->mOut&TR_RING_MASK];
624 if(!isdigit(ch)) break;
625 if(len < 64)
626 temp[len] = ch;
627 len++;
628 digis++;
629 tr->mOut++;
631 tr->mColumn += len;
632 if(digis > 0 && ch != '.' && !isalpha(ch))
634 if(len > 64)
636 TrErrorAt(tr, tr->mLine, col, "Integer is too long.");
637 return 0;
639 temp[len] = '\0';
640 *value = strtol(temp, NULL, 10);
641 if(*value < loBound || *value > hiBound)
643 TrErrorAt(tr, tr->mLine, col, "Expected a value from %d to %d.\n", loBound, hiBound);
644 return (0);
646 return (1);
649 TrErrorAt(tr, tr->mLine, col, "Expected an integer.\n");
650 return 0;
653 // Reads and validates (including bounds) a float token.
654 static int TrReadFloat(TokenReaderT *tr, const double loBound, const double hiBound, double *value)
656 uint col, digis, len;
657 char ch, temp[64+1];
659 col = tr->mColumn;
660 if(TrSkipWhitespace(tr))
662 col = tr->mColumn;
663 len = 0;
664 ch = tr->mRing[tr->mOut&TR_RING_MASK];
665 if(ch == '+' || ch == '-')
667 temp[len] = ch;
668 len++;
669 tr->mOut++;
672 digis = 0;
673 while(TrLoad(tr))
675 ch = tr->mRing[tr->mOut&TR_RING_MASK];
676 if(!isdigit(ch)) break;
677 if(len < 64)
678 temp[len] = ch;
679 len++;
680 digis++;
681 tr->mOut++;
683 if(ch == '.')
685 if(len < 64)
686 temp[len] = ch;
687 len++;
688 tr->mOut++;
690 while(TrLoad(tr))
692 ch = tr->mRing[tr->mOut&TR_RING_MASK];
693 if(!isdigit(ch)) break;
694 if(len < 64)
695 temp[len] = ch;
696 len++;
697 digis++;
698 tr->mOut++;
700 if(digis > 0)
702 if(ch == 'E' || ch == 'e')
704 if(len < 64)
705 temp[len] = ch;
706 len++;
707 digis = 0;
708 tr->mOut++;
709 if(ch == '+' || ch == '-')
711 if(len < 64)
712 temp[len] = ch;
713 len++;
714 tr->mOut++;
716 while(TrLoad(tr))
718 ch = tr->mRing[tr->mOut&TR_RING_MASK];
719 if(!isdigit(ch)) break;
720 if(len < 64)
721 temp[len] = ch;
722 len++;
723 digis++;
724 tr->mOut++;
727 tr->mColumn += len;
728 if(digis > 0 && ch != '.' && !isalpha(ch))
730 if(len > 64)
732 TrErrorAt(tr, tr->mLine, col, "Float is too long.");
733 return 0;
735 temp[len] = '\0';
736 *value = strtod(temp, NULL);
737 if(*value < loBound || *value > hiBound)
739 TrErrorAt (tr, tr->mLine, col, "Expected a value from %f to %f.\n", loBound, hiBound);
740 return 0;
742 return 1;
745 else
746 tr->mColumn += len;
748 TrErrorAt(tr, tr->mLine, col, "Expected a float.\n");
749 return 0;
752 // Reads and validates a string token.
753 static int TrReadString(TokenReaderT *tr, const uint maxLen, char *text)
755 uint col, len;
756 char ch;
758 col = tr->mColumn;
759 if(TrSkipWhitespace(tr))
761 col = tr->mColumn;
762 ch = tr->mRing[tr->mOut&TR_RING_MASK];
763 if(ch == '\"')
765 tr->mOut++;
766 len = 0;
767 while(TrLoad(tr))
769 ch = tr->mRing[tr->mOut&TR_RING_MASK];
770 tr->mOut++;
771 if(ch == '\"')
772 break;
773 if(ch == '\n')
775 TrErrorAt (tr, tr->mLine, col, "Unterminated string at end of line.\n");
776 return 0;
778 if(len < maxLen)
779 text[len] = ch;
780 len++;
782 if(ch != '\"')
784 tr->mColumn += 1 + len;
785 TrErrorAt(tr, tr->mLine, col, "Unterminated string at end of input.\n");
786 return 0;
788 tr->mColumn += 2 + len;
789 if(len > maxLen)
791 TrErrorAt (tr, tr->mLine, col, "String is too long.\n");
792 return 0;
794 text[len] = '\0';
795 return 1;
798 TrErrorAt(tr, tr->mLine, col, "Expected a string.\n");
799 return 0;
802 // Reads and validates the given operator.
803 static int TrReadOperator(TokenReaderT *tr, const char *op)
805 uint col, len;
806 char ch;
808 col = tr->mColumn;
809 if(TrSkipWhitespace(tr))
811 col = tr->mColumn;
812 len = 0;
813 while(op[len] != '\0' && TrLoad(tr))
815 ch = tr->mRing[tr->mOut&TR_RING_MASK];
816 if(ch != op[len]) break;
817 len++;
818 tr->mOut++;
820 tr->mColumn += len;
821 if(op[len] == '\0')
822 return 1;
824 TrErrorAt(tr, tr->mLine, col, "Expected '%s' operator.\n", op);
825 return 0;
828 /* Performs a string substitution. Any case-insensitive occurrences of the
829 * pattern string are replaced with the replacement string. The result is
830 * truncated if necessary.
832 static int StrSubst(const char *in, const char *pat, const char *rep, const size_t maxLen, char *out)
834 size_t inLen, patLen, repLen;
835 size_t si, di;
836 int truncated;
838 inLen = strlen(in);
839 patLen = strlen(pat);
840 repLen = strlen(rep);
841 si = 0;
842 di = 0;
843 truncated = 0;
844 while(si < inLen && di < maxLen)
846 if(patLen <= inLen-si)
848 if(strncasecmp(&in[si], pat, patLen) == 0)
850 if(repLen > maxLen-di)
852 repLen = maxLen - di;
853 truncated = 1;
855 strncpy(&out[di], rep, repLen);
856 si += patLen;
857 di += repLen;
860 out[di] = in[si];
861 si++;
862 di++;
864 if(si < inLen)
865 truncated = 1;
866 out[di] = '\0';
867 return !truncated;
871 /*********************
872 *** Math routines ***
873 *********************/
875 // Provide missing math routines for MSVC versions < 1800 (Visual Studio 2013).
876 #if defined(_MSC_VER) && _MSC_VER < 1800
877 static double round(double val)
879 if(val < 0.0)
880 return ceil(val-0.5);
881 return floor(val+0.5);
884 static double fmin(double a, double b)
886 return (a<b) ? a : b;
889 static double fmax(double a, double b)
891 return (a>b) ? a : b;
893 #endif
895 // Simple clamp routine.
896 static double Clamp(const double val, const double lower, const double upper)
898 return fmin(fmax(val, lower), upper);
901 // Performs linear interpolation.
902 static double Lerp(const double a, const double b, const double f)
904 return a + (f * (b - a));
907 static inline uint dither_rng(uint *seed)
909 *seed = (*seed * 96314165) + 907633515;
910 return *seed;
913 // Performs a triangular probability density function dither. It assumes the
914 // input sample is already scaled.
915 static inline double TpdfDither(const double in, uint *seed)
917 static const double PRNG_SCALE = 1.0 / UINT_MAX;
918 uint prn0, prn1;
920 prn0 = dither_rng(seed);
921 prn1 = dither_rng(seed);
922 return round(in + (prn0*PRNG_SCALE - prn1*PRNG_SCALE));
925 // Allocates an array of doubles.
926 static double *CreateArray(size_t n)
928 double *a;
930 if(n == 0) n = 1;
931 a = calloc(n, sizeof(double));
932 if(a == NULL)
934 fprintf(stderr, "Error: Out of memory.\n");
935 exit(-1);
937 return a;
940 // Frees an array of doubles.
941 static void DestroyArray(double *a)
942 { free(a); }
944 /* Fast Fourier transform routines. The number of points must be a power of
945 * two. In-place operation is possible only if both the real and imaginary
946 * parts are in-place together.
949 // Performs bit-reversal ordering.
950 static void FftArrange(const uint n, const Complex *in, Complex *out)
952 uint rk, k, m;
954 if(in == out)
956 // Handle in-place arrangement.
957 rk = 0;
958 for(k = 0;k < n;k++)
960 if(rk > k)
962 Complex temp = in[rk];
963 out[rk] = in[k];
964 out[k] = temp;
966 m = n;
967 while(rk&(m >>= 1))
968 rk &= ~m;
969 rk |= m;
972 else
974 // Handle copy arrangement.
975 rk = 0;
976 for(k = 0;k < n;k++)
978 out[rk] = in[k];
979 m = n;
980 while(rk&(m >>= 1))
981 rk &= ~m;
982 rk |= m;
987 // Performs the summation.
988 static void FftSummation(const int n, const double s, Complex *cplx)
990 double pi;
991 int m, m2;
992 int i, k, mk;
994 pi = s * M_PI;
995 for(m = 1, m2 = 2;m < n; m <<= 1, m2 <<= 1)
997 // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
998 double sm = sin(0.5 * pi / m);
999 Complex v = MakeComplex(-2.0*sm*sm, -sin(pi / m));
1000 Complex w = MakeComplex(1.0, 0.0);
1001 for(i = 0;i < m;i++)
1003 for(k = i;k < n;k += m2)
1005 Complex t;
1006 mk = k + m;
1007 t = c_mul(w, cplx[mk]);
1008 cplx[mk] = c_sub(cplx[k], t);
1009 cplx[k] = c_add(cplx[k], t);
1011 w = c_add(w, c_mul(v, w));
1016 // Performs a forward FFT.
1017 static void FftForward(const uint n, const Complex *in, Complex *out)
1019 FftArrange(n, in, out);
1020 FftSummation(n, 1.0, out);
1023 // Performs an inverse FFT.
1024 static void FftInverse(const uint n, const Complex *in, Complex *out)
1026 double f;
1027 uint i;
1029 FftArrange(n, in, out);
1030 FftSummation(n, -1.0, out);
1031 f = 1.0 / n;
1032 for(i = 0;i < n;i++)
1033 out[i] = c_muls(out[i], f);
1036 /* Calculate the complex helical sequence (or discrete-time analytical signal)
1037 * of the given input using the Hilbert transform. Given the natural logarithm
1038 * of a signal's magnitude response, the imaginary components can be used as
1039 * the angles for minimum-phase reconstruction.
1041 static void Hilbert(const uint n, const Complex *in, Complex *out)
1043 uint i;
1045 if(in == out)
1047 // Handle in-place operation.
1048 for(i = 0;i < n;i++)
1049 out[i].Imag = 0.0;
1051 else
1053 // Handle copy operation.
1054 for(i = 0;i < n;i++)
1055 out[i] = MakeComplex(in[i].Real, 0.0);
1057 FftInverse(n, out, out);
1058 for(i = 1;i < (n+1)/2;i++)
1059 out[i] = c_muls(out[i], 2.0);
1060 /* Increment i if n is even. */
1061 i += (n&1)^1;
1062 for(;i < n;i++)
1063 out[i] = MakeComplex(0.0, 0.0);
1064 FftForward(n, out, out);
1067 /* Calculate the magnitude response of the given input. This is used in
1068 * place of phase decomposition, since the phase residuals are discarded for
1069 * minimum phase reconstruction. The mirrored half of the response is also
1070 * discarded.
1072 static void MagnitudeResponse(const uint n, const Complex *in, double *out)
1074 const uint m = 1 + (n / 2);
1075 uint i;
1076 for(i = 0;i < m;i++)
1077 out[i] = fmax(c_abs(in[i]), EPSILON);
1080 /* Apply a range limit (in dB) to the given magnitude response. This is used
1081 * to adjust the effects of the diffuse-field average on the equalization
1082 * process.
1084 static void LimitMagnitudeResponse(const uint n, const double limit, const double *in, double *out)
1086 const uint m = 1 + (n / 2);
1087 double halfLim;
1088 uint i, lower, upper;
1089 double ave;
1091 halfLim = limit / 2.0;
1092 // Convert the response to dB.
1093 for(i = 0;i < m;i++)
1094 out[i] = 20.0 * log10(in[i]);
1095 // Use six octaves to calculate the average magnitude of the signal.
1096 lower = ((uint)ceil(n / pow(2.0, 8.0))) - 1;
1097 upper = ((uint)floor(n / pow(2.0, 2.0))) - 1;
1098 ave = 0.0;
1099 for(i = lower;i <= upper;i++)
1100 ave += out[i];
1101 ave /= upper - lower + 1;
1102 // Keep the response within range of the average magnitude.
1103 for(i = 0;i < m;i++)
1104 out[i] = Clamp(out[i], ave - halfLim, ave + halfLim);
1105 // Convert the response back to linear magnitude.
1106 for(i = 0;i < m;i++)
1107 out[i] = pow(10.0, out[i] / 20.0);
1110 /* Reconstructs the minimum-phase component for the given magnitude response
1111 * of a signal. This is equivalent to phase recomposition, sans the missing
1112 * residuals (which were discarded). The mirrored half of the response is
1113 * reconstructed.
1115 static void MinimumPhase(const uint n, const double *in, Complex *out)
1117 const uint m = 1 + (n / 2);
1118 double *mags;
1119 uint i;
1121 mags = CreateArray(n);
1122 for(i = 0;i < m;i++)
1124 mags[i] = fmax(EPSILON, in[i]);
1125 out[i] = MakeComplex(log(mags[i]), 0.0);
1127 for(;i < n;i++)
1129 mags[i] = mags[n - i];
1130 out[i] = out[n - i];
1132 Hilbert(n, out, out);
1133 // Remove any DC offset the filter has.
1134 mags[0] = EPSILON;
1135 for(i = 0;i < n;i++)
1137 Complex a = c_exp(MakeComplex(0.0, out[i].Imag));
1138 out[i] = c_mul(MakeComplex(mags[i], 0.0), a);
1140 DestroyArray(mags);
1144 /***************************
1145 *** Resampler functions ***
1146 ***************************/
1148 /* This is the normalized cardinal sine (sinc) function.
1150 * sinc(x) = { 1, x = 0
1151 * { sin(pi x) / (pi x), otherwise.
1153 static double Sinc(const double x)
1155 if(fabs(x) < EPSILON)
1156 return 1.0;
1157 return sin(M_PI * x) / (M_PI * x);
1160 /* The zero-order modified Bessel function of the first kind, used for the
1161 * Kaiser window.
1163 * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
1164 * = sum_{k=0}^inf ((x / 2)^k / k!)^2
1166 static double BesselI_0(const double x)
1168 double term, sum, x2, y, last_sum;
1169 int k;
1171 // Start at k=1 since k=0 is trivial.
1172 term = 1.0;
1173 sum = 1.0;
1174 x2 = x/2.0;
1175 k = 1;
1177 // Let the integration converge until the term of the sum is no longer
1178 // significant.
1179 do {
1180 y = x2 / k;
1181 k++;
1182 last_sum = sum;
1183 term *= y * y;
1184 sum += term;
1185 } while(sum != last_sum);
1186 return sum;
1189 /* Calculate a Kaiser window from the given beta value and a normalized k
1190 * [-1, 1].
1192 * w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
1193 * { 0, elsewhere.
1195 * Where k can be calculated as:
1197 * k = i / l, where -l <= i <= l.
1199 * or:
1201 * k = 2 i / M - 1, where 0 <= i <= M.
1203 static double Kaiser(const double b, const double k)
1205 if(!(k >= -1.0 && k <= 1.0))
1206 return 0.0;
1207 return BesselI_0(b * sqrt(1.0 - k*k)) / BesselI_0(b);
1210 // Calculates the greatest common divisor of a and b.
1211 static uint Gcd(uint x, uint y)
1213 while(y > 0)
1215 uint z = y;
1216 y = x % y;
1217 x = z;
1219 return x;
1222 /* Calculates the size (order) of the Kaiser window. Rejection is in dB and
1223 * the transition width is normalized frequency (0.5 is nyquist).
1225 * M = { ceil((r - 7.95) / (2.285 2 pi f_t)), r > 21
1226 * { ceil(5.79 / 2 pi f_t), r <= 21.
1229 static uint CalcKaiserOrder(const double rejection, const double transition)
1231 double w_t = 2.0 * M_PI * transition;
1232 if(rejection > 21.0)
1233 return (uint)ceil((rejection - 7.95) / (2.285 * w_t));
1234 return (uint)ceil(5.79 / w_t);
1237 // Calculates the beta value of the Kaiser window. Rejection is in dB.
1238 static double CalcKaiserBeta(const double rejection)
1240 if(rejection > 50.0)
1241 return 0.1102 * (rejection - 8.7);
1242 if(rejection >= 21.0)
1243 return (0.5842 * pow(rejection - 21.0, 0.4)) +
1244 (0.07886 * (rejection - 21.0));
1245 return 0.0;
1248 /* Calculates a point on the Kaiser-windowed sinc filter for the given half-
1249 * width, beta, gain, and cutoff. The point is specified in non-normalized
1250 * samples, from 0 to M, where M = (2 l + 1).
1252 * w(k) 2 p f_t sinc(2 f_t x)
1254 * x -- centered sample index (i - l)
1255 * k -- normalized and centered window index (x / l)
1256 * w(k) -- window function (Kaiser)
1257 * p -- gain compensation factor when sampling
1258 * f_t -- normalized center frequency (or cutoff; 0.5 is nyquist)
1260 static double SincFilter(const int l, const double b, const double gain, const double cutoff, const int i)
1262 return Kaiser(b, (double)(i - l) / l) * 2.0 * gain * cutoff * Sinc(2.0 * cutoff * (i - l));
1265 /* This is a polyphase sinc-filtered resampler.
1267 * Upsample Downsample
1269 * p/q = 3/2 p/q = 3/5
1271 * M-+-+-+-> M-+-+-+->
1272 * -------------------+ ---------------------+
1273 * p s * f f f f|f| | p s * f f f f f |
1274 * | 0 * 0 0 0|0|0 | | 0 * 0 0 0 0|0| |
1275 * v 0 * 0 0|0|0 0 | v 0 * 0 0 0|0|0 |
1276 * s * f|f|f f f | s * f f|f|f f |
1277 * 0 * |0|0 0 0 0 | 0 * 0|0|0 0 0 |
1278 * --------+=+--------+ 0 * |0|0 0 0 0 |
1279 * d . d .|d|. d . d ----------+=+--------+
1280 * d . . . .|d|. . . .
1281 * q->
1282 * q-+-+-+->
1284 * P_f(i,j) = q i mod p + pj
1285 * P_s(i,j) = floor(q i / p) - j
1286 * d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
1287 * { f[P_f(i,j)] s[P_s(i,j)], P_f(i,j) < M
1288 * { 0, P_f(i,j) >= M. }
1291 // Calculate the resampling metrics and build the Kaiser-windowed sinc filter
1292 // that's used to cut frequencies above the destination nyquist.
1293 static void ResamplerSetup(ResamplerT *rs, const uint srcRate, const uint dstRate)
1295 double cutoff, width, beta;
1296 uint gcd, l;
1297 int i;
1299 gcd = Gcd(srcRate, dstRate);
1300 rs->mP = dstRate / gcd;
1301 rs->mQ = srcRate / gcd;
1302 /* The cutoff is adjusted by half the transition width, so the transition
1303 * ends before the nyquist (0.5). Both are scaled by the downsampling
1304 * factor.
1306 if(rs->mP > rs->mQ)
1308 cutoff = 0.475 / rs->mP;
1309 width = 0.05 / rs->mP;
1311 else
1313 cutoff = 0.475 / rs->mQ;
1314 width = 0.05 / rs->mQ;
1316 // A rejection of -180 dB is used for the stop band. Round up when
1317 // calculating the left offset to avoid increasing the transition width.
1318 l = (CalcKaiserOrder(180.0, width)+1) / 2;
1319 beta = CalcKaiserBeta(180.0);
1320 rs->mM = l*2 + 1;
1321 rs->mL = l;
1322 rs->mF = CreateArray(rs->mM);
1323 for(i = 0;i < ((int)rs->mM);i++)
1324 rs->mF[i] = SincFilter((int)l, beta, rs->mP, cutoff, i);
1327 // Clean up after the resampler.
1328 static void ResamplerClear(ResamplerT *rs)
1330 DestroyArray(rs->mF);
1331 rs->mF = NULL;
1334 // Perform the upsample-filter-downsample resampling operation using a
1335 // polyphase filter implementation.
1336 static void ResamplerRun(ResamplerT *rs, const uint inN, const double *in, const uint outN, double *out)
1338 const uint p = rs->mP, q = rs->mQ, m = rs->mM, l = rs->mL;
1339 const double *f = rs->mF;
1340 uint j_f, j_s;
1341 double *work;
1342 uint i;
1344 if(outN == 0)
1345 return;
1347 // Handle in-place operation.
1348 if(in == out)
1349 work = CreateArray(outN);
1350 else
1351 work = out;
1352 // Resample the input.
1353 for(i = 0;i < outN;i++)
1355 double r = 0.0;
1356 // Input starts at l to compensate for the filter delay. This will
1357 // drop any build-up from the first half of the filter.
1358 j_f = (l + (q * i)) % p;
1359 j_s = (l + (q * i)) / p;
1360 while(j_f < m)
1362 // Only take input when 0 <= j_s < inN. This single unsigned
1363 // comparison catches both cases.
1364 if(j_s < inN)
1365 r += f[j_f] * in[j_s];
1366 j_f += p;
1367 j_s--;
1369 work[i] = r;
1371 // Clean up after in-place operation.
1372 if(work != out)
1374 for(i = 0;i < outN;i++)
1375 out[i] = work[i];
1376 DestroyArray(work);
1380 /*************************
1381 *** File source input ***
1382 *************************/
1384 // Read a binary value of the specified byte order and byte size from a file,
1385 // storing it as a 32-bit unsigned integer.
1386 static int ReadBin4(FILE *fp, const char *filename, const ByteOrderT order, const uint bytes, uint32 *out)
1388 uint8 in[4];
1389 uint32 accum;
1390 uint i;
1392 if(fread(in, 1, bytes, fp) != bytes)
1394 fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
1395 return 0;
1397 accum = 0;
1398 switch(order)
1400 case BO_LITTLE:
1401 for(i = 0;i < bytes;i++)
1402 accum = (accum<<8) | in[bytes - i - 1];
1403 break;
1404 case BO_BIG:
1405 for(i = 0;i < bytes;i++)
1406 accum = (accum<<8) | in[i];
1407 break;
1408 default:
1409 break;
1411 *out = accum;
1412 return 1;
1415 // Read a binary value of the specified byte order from a file, storing it as
1416 // a 64-bit unsigned integer.
1417 static int ReadBin8(FILE *fp, const char *filename, const ByteOrderT order, uint64 *out)
1419 uint8 in [8];
1420 uint64 accum;
1421 uint i;
1423 if(fread(in, 1, 8, fp) != 8)
1425 fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
1426 return 0;
1428 accum = 0ULL;
1429 switch(order)
1431 case BO_LITTLE:
1432 for(i = 0;i < 8;i++)
1433 accum = (accum<<8) | in[8 - i - 1];
1434 break;
1435 case BO_BIG:
1436 for(i = 0;i < 8;i++)
1437 accum = (accum<<8) | in[i];
1438 break;
1439 default:
1440 break;
1442 *out = accum;
1443 return 1;
1446 /* Read a binary value of the specified type, byte order, and byte size from
1447 * a file, converting it to a double. For integer types, the significant
1448 * bits are used to normalize the result. The sign of bits determines
1449 * whether they are padded toward the MSB (negative) or LSB (positive).
1450 * Floating-point types are not normalized.
1452 static int ReadBinAsDouble(FILE *fp, const char *filename, const ByteOrderT order, const ElementTypeT type, const uint bytes, const int bits, double *out)
1454 union {
1455 uint32 ui;
1456 int32 i;
1457 float f;
1458 } v4;
1459 union {
1460 uint64 ui;
1461 double f;
1462 } v8;
1464 *out = 0.0;
1465 if(bytes > 4)
1467 if(!ReadBin8(fp, filename, order, &v8.ui))
1468 return 0;
1469 if(type == ET_FP)
1470 *out = v8.f;
1472 else
1474 if(!ReadBin4(fp, filename, order, bytes, &v4.ui))
1475 return 0;
1476 if(type == ET_FP)
1477 *out = v4.f;
1478 else
1480 if(bits > 0)
1481 v4.ui >>= (8*bytes) - ((uint)bits);
1482 else
1483 v4.ui &= (0xFFFFFFFF >> (32+bits));
1485 if(v4.ui&(uint)(1<<(abs(bits)-1)))
1486 v4.ui |= (0xFFFFFFFF << abs (bits));
1487 *out = v4.i / (double)(1<<(abs(bits)-1));
1490 return 1;
1493 /* Read an ascii value of the specified type from a file, converting it to a
1494 * double. For integer types, the significant bits are used to normalize the
1495 * result. The sign of the bits should always be positive. This also skips
1496 * up to one separator character before the element itself.
1498 static int ReadAsciiAsDouble(TokenReaderT *tr, const char *filename, const ElementTypeT type, const uint bits, double *out)
1500 if(TrIsOperator(tr, ","))
1501 TrReadOperator(tr, ",");
1502 else if(TrIsOperator(tr, ":"))
1503 TrReadOperator(tr, ":");
1504 else if(TrIsOperator(tr, ";"))
1505 TrReadOperator(tr, ";");
1506 else if(TrIsOperator(tr, "|"))
1507 TrReadOperator(tr, "|");
1509 if(type == ET_FP)
1511 if(!TrReadFloat(tr, -HUGE_VAL, HUGE_VAL, out))
1513 fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
1514 return 0;
1517 else
1519 int v;
1520 if(!TrReadInt(tr, -(1<<(bits-1)), (1<<(bits-1))-1, &v))
1522 fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
1523 return 0;
1525 *out = v / (double)((1<<(bits-1))-1);
1527 return 1;
1530 // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
1531 // the source parameters and data set metrics.
1532 static int ReadWaveFormat(FILE *fp, const ByteOrderT order, const uint hrirRate, SourceRefT *src)
1534 uint32 fourCC, chunkSize;
1535 uint32 format, channels, rate, dummy, block, size, bits;
1537 chunkSize = 0;
1538 do {
1539 if (chunkSize > 0)
1540 fseek (fp, (long) chunkSize, SEEK_CUR);
1541 if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
1542 !ReadBin4(fp, src->mPath, order, 4, &chunkSize))
1543 return 0;
1544 } while(fourCC != FOURCC_FMT);
1545 if(!ReadBin4(fp, src->mPath, order, 2, & format) ||
1546 !ReadBin4(fp, src->mPath, order, 2, & channels) ||
1547 !ReadBin4(fp, src->mPath, order, 4, & rate) ||
1548 !ReadBin4(fp, src->mPath, order, 4, & dummy) ||
1549 !ReadBin4(fp, src->mPath, order, 2, & block))
1550 return (0);
1551 block /= channels;
1552 if(chunkSize > 14)
1554 if(!ReadBin4(fp, src->mPath, order, 2, &size))
1555 return 0;
1556 size /= 8;
1557 if(block > size)
1558 size = block;
1560 else
1561 size = block;
1562 if(format == WAVE_FORMAT_EXTENSIBLE)
1564 fseek(fp, 2, SEEK_CUR);
1565 if(!ReadBin4(fp, src->mPath, order, 2, &bits))
1566 return 0;
1567 if(bits == 0)
1568 bits = 8 * size;
1569 fseek(fp, 4, SEEK_CUR);
1570 if(!ReadBin4(fp, src->mPath, order, 2, &format))
1571 return 0;
1572 fseek(fp, (long)(chunkSize - 26), SEEK_CUR);
1574 else
1576 bits = 8 * size;
1577 if(chunkSize > 14)
1578 fseek(fp, (long)(chunkSize - 16), SEEK_CUR);
1579 else
1580 fseek(fp, (long)(chunkSize - 14), SEEK_CUR);
1582 if(format != WAVE_FORMAT_PCM && format != WAVE_FORMAT_IEEE_FLOAT)
1584 fprintf(stderr, "Error: Unsupported WAVE format in file '%s'.\n", src->mPath);
1585 return 0;
1587 if(src->mChannel >= channels)
1589 fprintf(stderr, "Error: Missing source channel in WAVE file '%s'.\n", src->mPath);
1590 return 0;
1592 if(rate != hrirRate)
1594 fprintf(stderr, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src->mPath);
1595 return 0;
1597 if(format == WAVE_FORMAT_PCM)
1599 if(size < 2 || size > 4)
1601 fprintf(stderr, "Error: Unsupported sample size in WAVE file '%s'.\n", src->mPath);
1602 return 0;
1604 if(bits < 16 || bits > (8*size))
1606 fprintf (stderr, "Error: Bad significant bits in WAVE file '%s'.\n", src->mPath);
1607 return 0;
1609 src->mType = ET_INT;
1611 else
1613 if(size != 4 && size != 8)
1615 fprintf(stderr, "Error: Unsupported sample size in WAVE file '%s'.\n", src->mPath);
1616 return 0;
1618 src->mType = ET_FP;
1620 src->mSize = size;
1621 src->mBits = (int)bits;
1622 src->mSkip = channels;
1623 return 1;
1626 // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
1627 static int ReadWaveData(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
1629 int pre, post, skip;
1630 uint i;
1632 pre = (int)(src->mSize * src->mChannel);
1633 post = (int)(src->mSize * (src->mSkip - src->mChannel - 1));
1634 skip = 0;
1635 for(i = 0;i < n;i++)
1637 skip += pre;
1638 if(skip > 0)
1639 fseek(fp, skip, SEEK_CUR);
1640 if(!ReadBinAsDouble(fp, src->mPath, order, src->mType, src->mSize, src->mBits, &hrir[i]))
1641 return 0;
1642 skip = post;
1644 if(skip > 0)
1645 fseek(fp, skip, SEEK_CUR);
1646 return 1;
1649 // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
1650 // doubles.
1651 static int ReadWaveList(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
1653 uint32 fourCC, chunkSize, listSize, count;
1654 uint block, skip, offset, i;
1655 double lastSample;
1657 for (;;) {
1658 if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, & fourCC) ||
1659 !ReadBin4(fp, src->mPath, order, 4, & chunkSize))
1660 return (0);
1662 if(fourCC == FOURCC_DATA)
1664 block = src->mSize * src->mSkip;
1665 count = chunkSize / block;
1666 if(count < (src->mOffset + n))
1668 fprintf(stderr, "Error: Bad read from file '%s'.\n", src->mPath);
1669 return 0;
1671 fseek(fp, (long)(src->mOffset * block), SEEK_CUR);
1672 if(!ReadWaveData(fp, src, order, n, &hrir[0]))
1673 return 0;
1674 return 1;
1676 else if(fourCC == FOURCC_LIST)
1678 if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC))
1679 return 0;
1680 chunkSize -= 4;
1681 if(fourCC == FOURCC_WAVL)
1682 break;
1684 if(chunkSize > 0)
1685 fseek(fp, (long)chunkSize, SEEK_CUR);
1687 listSize = chunkSize;
1688 block = src->mSize * src->mSkip;
1689 skip = src->mOffset;
1690 offset = 0;
1691 lastSample = 0.0;
1692 while(offset < n && listSize > 8)
1694 if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
1695 !ReadBin4(fp, src->mPath, order, 4, &chunkSize))
1696 return 0;
1697 listSize -= 8 + chunkSize;
1698 if(fourCC == FOURCC_DATA)
1700 count = chunkSize / block;
1701 if(count > skip)
1703 fseek(fp, (long)(skip * block), SEEK_CUR);
1704 chunkSize -= skip * block;
1705 count -= skip;
1706 skip = 0;
1707 if(count > (n - offset))
1708 count = n - offset;
1709 if(!ReadWaveData(fp, src, order, count, &hrir[offset]))
1710 return 0;
1711 chunkSize -= count * block;
1712 offset += count;
1713 lastSample = hrir [offset - 1];
1715 else
1717 skip -= count;
1718 count = 0;
1721 else if(fourCC == FOURCC_SLNT)
1723 if(!ReadBin4(fp, src->mPath, order, 4, &count))
1724 return 0;
1725 chunkSize -= 4;
1726 if(count > skip)
1728 count -= skip;
1729 skip = 0;
1730 if(count > (n - offset))
1731 count = n - offset;
1732 for(i = 0; i < count; i ++)
1733 hrir[offset + i] = lastSample;
1734 offset += count;
1736 else
1738 skip -= count;
1739 count = 0;
1742 if(chunkSize > 0)
1743 fseek(fp, (long)chunkSize, SEEK_CUR);
1745 if(offset < n)
1747 fprintf(stderr, "Error: Bad read from file '%s'.\n", src->mPath);
1748 return 0;
1750 return 1;
1753 // Load a source HRIR from a RIFF/RIFX WAVE file.
1754 static int LoadWaveSource(FILE *fp, SourceRefT *src, const uint hrirRate, const uint n, double *hrir)
1756 uint32 fourCC, dummy;
1757 ByteOrderT order;
1759 if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
1760 !ReadBin4(fp, src->mPath, BO_LITTLE, 4, &dummy))
1761 return 0;
1762 if(fourCC == FOURCC_RIFF)
1763 order = BO_LITTLE;
1764 else if(fourCC == FOURCC_RIFX)
1765 order = BO_BIG;
1766 else
1768 fprintf(stderr, "Error: No RIFF/RIFX chunk in file '%s'.\n", src->mPath);
1769 return 0;
1772 if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC))
1773 return 0;
1774 if(fourCC != FOURCC_WAVE)
1776 fprintf(stderr, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src->mPath);
1777 return 0;
1779 if(!ReadWaveFormat(fp, order, hrirRate, src))
1780 return 0;
1781 if(!ReadWaveList(fp, src, order, n, hrir))
1782 return 0;
1783 return 1;
1786 // Load a source HRIR from a binary file.
1787 static int LoadBinarySource(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
1789 uint i;
1791 fseek(fp, (long)src->mOffset, SEEK_SET);
1792 for(i = 0;i < n;i++)
1794 if(!ReadBinAsDouble(fp, src->mPath, order, src->mType, src->mSize, src->mBits, &hrir[i]))
1795 return 0;
1796 if(src->mSkip > 0)
1797 fseek(fp, (long)src->mSkip, SEEK_CUR);
1799 return 1;
1802 // Load a source HRIR from an ASCII text file containing a list of elements
1803 // separated by whitespace or common list operators (',', ';', ':', '|').
1804 static int LoadAsciiSource(FILE *fp, const SourceRefT *src, const uint n, double *hrir)
1806 TokenReaderT tr;
1807 uint i, j;
1808 double dummy;
1810 TrSetup(fp, NULL, &tr);
1811 for(i = 0;i < src->mOffset;i++)
1813 if(!ReadAsciiAsDouble(&tr, src->mPath, src->mType, (uint)src->mBits, &dummy))
1814 return (0);
1816 for(i = 0;i < n;i++)
1818 if(!ReadAsciiAsDouble(&tr, src->mPath, src->mType, (uint)src->mBits, &hrir[i]))
1819 return 0;
1820 for(j = 0;j < src->mSkip;j++)
1822 if(!ReadAsciiAsDouble(&tr, src->mPath, src->mType, (uint)src->mBits, &dummy))
1823 return 0;
1826 return 1;
1829 // Load a source HRIR from a supported file type.
1830 static int LoadSource(SourceRefT *src, const uint hrirRate, const uint n, double *hrir)
1832 int result;
1833 FILE *fp;
1835 if (src->mFormat == SF_ASCII)
1836 fp = fopen(src->mPath, "r");
1837 else
1838 fp = fopen(src->mPath, "rb");
1839 if(fp == NULL)
1841 fprintf(stderr, "Error: Could not open source file '%s'.\n", src->mPath);
1842 return 0;
1844 if(src->mFormat == SF_WAVE)
1845 result = LoadWaveSource(fp, src, hrirRate, n, hrir);
1846 else if(src->mFormat == SF_BIN_LE)
1847 result = LoadBinarySource(fp, src, BO_LITTLE, n, hrir);
1848 else if(src->mFormat == SF_BIN_BE)
1849 result = LoadBinarySource(fp, src, BO_BIG, n, hrir);
1850 else
1851 result = LoadAsciiSource(fp, src, n, hrir);
1852 fclose(fp);
1853 return result;
1857 /***************************
1858 *** File storage output ***
1859 ***************************/
1861 // Write an ASCII string to a file.
1862 static int WriteAscii(const char *out, FILE *fp, const char *filename)
1864 size_t len;
1866 len = strlen(out);
1867 if(fwrite(out, 1, len, fp) != len)
1869 fclose(fp);
1870 fprintf(stderr, "Error: Bad write to file '%s'.\n", filename);
1871 return 0;
1873 return 1;
1876 // Write a binary value of the given byte order and byte size to a file,
1877 // loading it from a 32-bit unsigned integer.
1878 static int WriteBin4(const ByteOrderT order, const uint bytes, const uint32 in, FILE *fp, const char *filename)
1880 uint8 out[4];
1881 uint i;
1883 switch(order)
1885 case BO_LITTLE:
1886 for(i = 0;i < bytes;i++)
1887 out[i] = (in>>(i*8)) & 0x000000FF;
1888 break;
1889 case BO_BIG:
1890 for(i = 0;i < bytes;i++)
1891 out[bytes - i - 1] = (in>>(i*8)) & 0x000000FF;
1892 break;
1893 default:
1894 break;
1896 if(fwrite(out, 1, bytes, fp) != bytes)
1898 fprintf(stderr, "Error: Bad write to file '%s'.\n", filename);
1899 return 0;
1901 return 1;
1904 // Store the OpenAL Soft HRTF data set.
1905 static int StoreMhr(const HrirDataT *hData, const int experimental, const char *filename)
1907 uint e, step, end, n, j, i;
1908 uint dither_seed;
1909 FILE *fp;
1910 int v;
1912 if((fp=fopen(filename, "wb")) == NULL)
1914 fprintf(stderr, "Error: Could not open MHR file '%s'.\n", filename);
1915 return 0;
1917 if(!WriteAscii(experimental ? MHR_FORMAT_EXPERIMENTAL : MHR_FORMAT, fp, filename))
1918 return 0;
1919 if(!WriteBin4(BO_LITTLE, 4, (uint32)hData->mIrRate, fp, filename))
1920 return 0;
1921 if(experimental)
1923 if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mSampleType, fp, filename))
1924 return 0;
1925 if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mChannelType, fp, filename))
1926 return 0;
1928 if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mIrPoints, fp, filename))
1929 return 0;
1930 if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mEvCount, fp, filename))
1931 return 0;
1932 for(e = 0;e < hData->mEvCount;e++)
1934 if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mAzCount[e], fp, filename))
1935 return 0;
1937 step = hData->mIrSize;
1938 end = hData->mIrCount * step;
1939 n = hData->mIrPoints;
1940 dither_seed = 22222;
1941 for(j = 0;j < end;j += step)
1943 const double scale = (!experimental || hData->mSampleType == ST_S16) ? 32767.0 :
1944 ((hData->mSampleType == ST_S24) ? 8388607.0 : 0.0);
1945 const int bps = (!experimental || hData->mSampleType == ST_S16) ? 2 :
1946 ((hData->mSampleType == ST_S24) ? 3 : 0);
1947 double out[MAX_TRUNCSIZE];
1948 for(i = 0;i < n;i++)
1949 out[i] = TpdfDither(scale * hData->mHrirs[j+i], &dither_seed);
1950 for(i = 0;i < n;i++)
1952 v = (int)Clamp(out[i], -scale-1.0, scale);
1953 if(!WriteBin4(BO_LITTLE, bps, (uint32)v, fp, filename))
1954 return 0;
1957 for(j = 0;j < hData->mIrCount;j++)
1959 v = (int)fmin(round(hData->mIrRate * hData->mHrtds[j]), MAX_HRTD);
1960 if(!WriteBin4(BO_LITTLE, 1, (uint32)v, fp, filename))
1961 return 0;
1963 fclose(fp);
1964 return 1;
1968 /***********************
1969 *** HRTF processing ***
1970 ***********************/
1972 // Calculate the onset time of an HRIR and average it with any existing
1973 // timing for its elevation and azimuth.
1974 static void AverageHrirOnset(const double *hrir, const double f, const uint ei, const uint ai, const HrirDataT *hData)
1976 double mag;
1977 uint n, i, j;
1979 mag = 0.0;
1980 n = hData->mIrPoints;
1981 for(i = 0;i < n;i++)
1982 mag = fmax(fabs(hrir[i]), mag);
1983 mag *= 0.15;
1984 for(i = 0;i < n;i++)
1986 if(fabs(hrir[i]) >= mag)
1987 break;
1989 j = hData->mEvOffset[ei] + ai;
1990 hData->mHrtds[j] = Lerp(hData->mHrtds[j], ((double)i) / hData->mIrRate, f);
1993 // Calculate the magnitude response of an HRIR and average it with any
1994 // existing responses for its elevation and azimuth.
1995 static void AverageHrirMagnitude(const double *hrir, const double f, const uint ei, const uint ai, const HrirDataT *hData)
1997 uint n, m, i, j;
1998 Complex *cplx;
1999 double *mags;
2001 n = hData->mFftSize;
2002 cplx = calloc(sizeof(*cplx), n);
2003 mags = calloc(sizeof(*mags), n);
2004 for(i = 0;i < hData->mIrPoints;i++)
2005 cplx[i] = MakeComplex(hrir[i], 0.0);
2006 for(;i < n;i++)
2007 cplx[i] = MakeComplex(0.0, 0.0);
2008 FftForward(n, cplx, cplx);
2009 MagnitudeResponse(n, cplx, mags);
2010 m = 1 + (n / 2);
2011 j = (hData->mEvOffset[ei] + ai) * hData->mIrSize;
2012 for(i = 0;i < m;i++)
2013 hData->mHrirs[j+i] = Lerp(hData->mHrirs[j+i], mags[i], f);
2014 free(mags);
2015 free(cplx);
2018 /* Calculate the contribution of each HRIR to the diffuse-field average based
2019 * on the area of its surface patch. All patches are centered at the HRIR
2020 * coordinates on the unit sphere and are measured by solid angle.
2022 static void CalculateDfWeights(const HrirDataT *hData, double *weights)
2024 double evs, sum, ev, up_ev, down_ev, solidAngle;
2025 uint ei;
2027 evs = 90.0 / (hData->mEvCount - 1);
2028 sum = 0.0;
2029 for(ei = hData->mEvStart;ei < hData->mEvCount;ei++)
2031 // For each elevation, calculate the upper and lower limits of the
2032 // patch band.
2033 ev = -90.0 + (ei * 2.0 * evs);
2034 if(ei < (hData->mEvCount - 1))
2035 up_ev = (ev + evs) * M_PI / 180.0;
2036 else
2037 up_ev = M_PI / 2.0;
2038 if(ei > 0)
2039 down_ev = (ev - evs) * M_PI / 180.0;
2040 else
2041 down_ev = -M_PI / 2.0;
2042 // Calculate the area of the patch band.
2043 solidAngle = 2.0 * M_PI * (sin(up_ev) - sin(down_ev));
2044 // Each weight is the area of one patch.
2045 weights[ei] = solidAngle / hData->mAzCount [ei];
2046 // Sum the total surface area covered by the HRIRs.
2047 sum += solidAngle;
2049 // Normalize the weights given the total surface coverage.
2050 for(ei = hData->mEvStart;ei < hData->mEvCount;ei++)
2051 weights[ei] /= sum;
2054 /* Calculate the diffuse-field average from the given magnitude responses of
2055 * the HRIR set. Weighting can be applied to compensate for the varying
2056 * surface area covered by each HRIR. The final average can then be limited
2057 * by the specified magnitude range (in positive dB; 0.0 to skip).
2059 static void CalculateDiffuseFieldAverage(const HrirDataT *hData, const int weighted, const double limit, double *dfa)
2061 uint ei, ai, count, step, start, end, m, j, i;
2062 double *weights;
2064 weights = CreateArray(hData->mEvCount);
2065 if(weighted)
2067 // Use coverage weighting to calculate the average.
2068 CalculateDfWeights(hData, weights);
2070 else
2072 // If coverage weighting is not used, the weights still need to be
2073 // averaged by the number of HRIRs.
2074 count = 0;
2075 for(ei = hData->mEvStart;ei < hData->mEvCount;ei++)
2076 count += hData->mAzCount [ei];
2077 for(ei = hData->mEvStart;ei < hData->mEvCount;ei++)
2078 weights[ei] = 1.0 / count;
2080 ei = hData->mEvStart;
2081 ai = 0;
2082 step = hData->mIrSize;
2083 start = hData->mEvOffset[ei] * step;
2084 end = hData->mIrCount * step;
2085 m = 1 + (hData->mFftSize / 2);
2086 for(i = 0;i < m;i++)
2087 dfa[i] = 0.0;
2088 for(j = start;j < end;j += step)
2090 // Get the weight for this HRIR's contribution.
2091 double weight = weights[ei];
2092 // Add this HRIR's weighted power average to the total.
2093 for(i = 0;i < m;i++)
2094 dfa[i] += weight * hData->mHrirs[j+i] * hData->mHrirs[j+i];
2095 // Determine the next weight to use.
2096 ai++;
2097 if(ai >= hData->mAzCount[ei])
2099 ei++;
2100 ai = 0;
2103 // Finish the average calculation and keep it from being too small.
2104 for(i = 0;i < m;i++)
2105 dfa[i] = fmax(sqrt(dfa[i]), EPSILON);
2106 // Apply a limit to the magnitude range of the diffuse-field average if
2107 // desired.
2108 if(limit > 0.0)
2109 LimitMagnitudeResponse(hData->mFftSize, limit, dfa, dfa);
2110 DestroyArray(weights);
2113 // Perform diffuse-field equalization on the magnitude responses of the HRIR
2114 // set using the given average response.
2115 static void DiffuseFieldEqualize(const double *dfa, const HrirDataT *hData)
2117 uint step, start, end, m, j, i;
2119 step = hData->mIrSize;
2120 start = hData->mEvOffset[hData->mEvStart] * step;
2121 end = hData->mIrCount * step;
2122 m = 1 + (hData->mFftSize / 2);
2123 for(j = start;j < end;j += step)
2125 for(i = 0;i < m;i++)
2126 hData->mHrirs[j+i] /= dfa[i];
2130 // Perform minimum-phase reconstruction using the magnitude responses of the
2131 // HRIR set.
2132 static void ReconstructHrirs(const HrirDataT *hData)
2134 uint step, start, end, n, j, i;
2135 uint pcdone, lastpc;
2136 Complex *cplx;
2138 pcdone = lastpc = 0;
2139 printf("%3d%% done.", pcdone);
2140 fflush(stdout);
2142 step = hData->mIrSize;
2143 start = hData->mEvOffset[hData->mEvStart] * step;
2144 end = hData->mIrCount * step;
2145 n = hData->mFftSize;
2146 cplx = calloc(sizeof(*cplx), n);
2147 for(j = start;j < end;j += step)
2149 MinimumPhase(n, &hData->mHrirs[j], cplx);
2150 FftInverse(n, cplx, cplx);
2151 for(i = 0;i < hData->mIrPoints;i++)
2152 hData->mHrirs[j+i] = cplx[i].Real;
2153 pcdone = (j+step-start) * 100 / (end-start);
2154 if(pcdone != lastpc)
2156 lastpc = pcdone;
2157 printf("\r%3d%% done.", pcdone);
2158 fflush(stdout);
2161 free(cplx);
2162 printf("\n");
2165 // Resamples the HRIRs for use at the given sampling rate.
2166 static void ResampleHrirs(const uint rate, HrirDataT *hData)
2168 uint n, step, start, end, j;
2169 ResamplerT rs;
2171 ResamplerSetup(&rs, hData->mIrRate, rate);
2172 n = hData->mIrPoints;
2173 step = hData->mIrSize;
2174 start = hData->mEvOffset[hData->mEvStart] * step;
2175 end = hData->mIrCount * step;
2176 for(j = start;j < end;j += step)
2177 ResamplerRun(&rs, n, &hData->mHrirs[j], n, &hData->mHrirs[j]);
2178 ResamplerClear(&rs);
2179 hData->mIrRate = rate;
2182 /* Given an elevation index and an azimuth, calculate the indices of the two
2183 * HRIRs that bound the coordinate along with a factor for calculating the
2184 * continous HRIR using interpolation.
2186 static void CalcAzIndices(const HrirDataT *hData, const uint ei, const double az, uint *j0, uint *j1, double *jf)
2188 double af;
2189 uint ai;
2191 af = ((2.0*M_PI) + az) * hData->mAzCount[ei] / (2.0*M_PI);
2192 ai = ((uint)af) % hData->mAzCount[ei];
2193 af -= floor(af);
2195 *j0 = hData->mEvOffset[ei] + ai;
2196 *j1 = hData->mEvOffset[ei] + ((ai+1) % hData->mAzCount [ei]);
2197 *jf = af;
2200 // Synthesize any missing onset timings at the bottom elevations. This just
2201 // blends between slightly exaggerated known onsets. Not an accurate model.
2202 static void SynthesizeOnsets(HrirDataT *hData)
2204 uint oi, e, a, j0, j1;
2205 double t, of, jf;
2207 oi = hData->mEvStart;
2208 t = 0.0;
2209 for(a = 0;a < hData->mAzCount[oi];a++)
2210 t += hData->mHrtds[hData->mEvOffset[oi] + a];
2211 hData->mHrtds[0] = 1.32e-4 + (t / hData->mAzCount[oi]);
2212 for(e = 1;e < hData->mEvStart;e++)
2214 of = ((double)e) / hData->mEvStart;
2215 for(a = 0;a < hData->mAzCount[e];a++)
2217 CalcAzIndices(hData, oi, a * 2.0 * M_PI / hData->mAzCount[e], &j0, &j1, &jf);
2218 hData->mHrtds[hData->mEvOffset[e] + a] = Lerp(hData->mHrtds[0], Lerp(hData->mHrtds[j0], hData->mHrtds[j1], jf), of);
2223 /* Attempt to synthesize any missing HRIRs at the bottom elevations. Right
2224 * now this just blends the lowest elevation HRIRs together and applies some
2225 * attenuation and high frequency damping. It is a simple, if inaccurate
2226 * model.
2228 static void SynthesizeHrirs (HrirDataT *hData)
2230 uint oi, a, e, step, n, i, j;
2231 double lp[4], s0, s1;
2232 double of, b;
2233 uint j0, j1;
2234 double jf;
2236 if(hData->mEvStart <= 0)
2237 return;
2238 step = hData->mIrSize;
2239 oi = hData->mEvStart;
2240 n = hData->mIrPoints;
2241 for(i = 0;i < n;i++)
2242 hData->mHrirs[i] = 0.0;
2243 for(a = 0;a < hData->mAzCount[oi];a++)
2245 j = (hData->mEvOffset[oi] + a) * step;
2246 for(i = 0;i < n;i++)
2247 hData->mHrirs[i] += hData->mHrirs[j+i] / hData->mAzCount[oi];
2249 for(e = 1;e < hData->mEvStart;e++)
2251 of = ((double)e) / hData->mEvStart;
2252 b = (1.0 - of) * (3.5e-6 * hData->mIrRate);
2253 for(a = 0;a < hData->mAzCount[e];a++)
2255 j = (hData->mEvOffset[e] + a) * step;
2256 CalcAzIndices(hData, oi, a * 2.0 * M_PI / hData->mAzCount[e], &j0, &j1, &jf);
2257 j0 *= step;
2258 j1 *= step;
2259 lp[0] = 0.0;
2260 lp[1] = 0.0;
2261 lp[2] = 0.0;
2262 lp[3] = 0.0;
2263 for(i = 0;i < n;i++)
2265 s0 = hData->mHrirs[i];
2266 s1 = Lerp(hData->mHrirs[j0+i], hData->mHrirs[j1+i], jf);
2267 s0 = Lerp(s0, s1, of);
2268 lp[0] = Lerp(s0, lp[0], b);
2269 lp[1] = Lerp(lp[0], lp[1], b);
2270 lp[2] = Lerp(lp[1], lp[2], b);
2271 lp[3] = Lerp(lp[2], lp[3], b);
2272 hData->mHrirs[j+i] = lp[3];
2276 b = 3.5e-6 * hData->mIrRate;
2277 lp[0] = 0.0;
2278 lp[1] = 0.0;
2279 lp[2] = 0.0;
2280 lp[3] = 0.0;
2281 for(i = 0;i < n;i++)
2283 s0 = hData->mHrirs[i];
2284 lp[0] = Lerp(s0, lp[0], b);
2285 lp[1] = Lerp(lp[0], lp[1], b);
2286 lp[2] = Lerp(lp[1], lp[2], b);
2287 lp[3] = Lerp(lp[2], lp[3], b);
2288 hData->mHrirs[i] = lp[3];
2290 hData->mEvStart = 0;
2293 // The following routines assume a full set of HRIRs for all elevations.
2295 // Normalize the HRIR set and slightly attenuate the result.
2296 static void NormalizeHrirs (const HrirDataT *hData)
2298 uint step, end, n, j, i;
2299 double maxLevel;
2301 step = hData->mIrSize;
2302 end = hData->mIrCount * step;
2303 n = hData->mIrPoints;
2304 maxLevel = 0.0;
2305 for(j = 0;j < end;j += step)
2307 for(i = 0;i < n;i++)
2308 maxLevel = fmax(fabs(hData->mHrirs[j+i]), maxLevel);
2310 maxLevel = 1.01 * maxLevel;
2311 for(j = 0;j < end;j += step)
2313 for(i = 0;i < n;i++)
2314 hData->mHrirs[j+i] /= maxLevel;
2318 // Calculate the left-ear time delay using a spherical head model.
2319 static double CalcLTD(const double ev, const double az, const double rad, const double dist)
2321 double azp, dlp, l, al;
2323 azp = asin(cos(ev) * sin(az));
2324 dlp = sqrt((dist*dist) + (rad*rad) + (2.0*dist*rad*sin(azp)));
2325 l = sqrt((dist*dist) - (rad*rad));
2326 al = (0.5 * M_PI) + azp;
2327 if(dlp > l)
2328 dlp = l + (rad * (al - acos(rad / dist)));
2329 return (dlp / 343.3);
2332 // Calculate the effective head-related time delays for each minimum-phase
2333 // HRIR.
2334 static void CalculateHrtds (const HeadModelT model, const double radius, HrirDataT *hData)
2336 double minHrtd, maxHrtd;
2337 uint e, a, j;
2338 double t;
2340 minHrtd = 1000.0;
2341 maxHrtd = -1000.0;
2342 for(e = 0;e < hData->mEvCount;e++)
2344 for(a = 0;a < hData->mAzCount[e];a++)
2346 j = hData->mEvOffset[e] + a;
2347 if(model == HM_DATASET)
2348 t = hData->mHrtds[j] * radius / hData->mRadius;
2349 else
2350 t = CalcLTD((-90.0 + (e * 180.0 / (hData->mEvCount - 1))) * M_PI / 180.0,
2351 (a * 360.0 / hData->mAzCount [e]) * M_PI / 180.0,
2352 radius, hData->mDistance);
2353 hData->mHrtds[j] = t;
2354 maxHrtd = fmax(t, maxHrtd);
2355 minHrtd = fmin(t, minHrtd);
2358 maxHrtd -= minHrtd;
2359 for(j = 0;j < hData->mIrCount;j++)
2360 hData->mHrtds[j] -= minHrtd;
2361 hData->mMaxHrtd = maxHrtd;
2365 // Process the data set definition to read and validate the data set metrics.
2366 static int ProcessMetrics(TokenReaderT *tr, const uint fftSize, const uint truncSize, HrirDataT *hData)
2368 int hasRate = 0, hasPoints = 0, hasAzimuths = 0;
2369 int hasRadius = 0, hasDistance = 0;
2370 char ident[MAX_IDENT_LEN+1];
2371 uint line, col;
2372 double fpVal;
2373 uint points;
2374 int intVal;
2376 while(!(hasRate && hasPoints && hasAzimuths && hasRadius && hasDistance))
2378 TrIndication(tr, & line, & col);
2379 if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
2380 return 0;
2381 if(strcasecmp(ident, "rate") == 0)
2383 if(hasRate)
2385 TrErrorAt(tr, line, col, "Redefinition of 'rate'.\n");
2386 return 0;
2388 if(!TrReadOperator(tr, "="))
2389 return 0;
2390 if(!TrReadInt(tr, MIN_RATE, MAX_RATE, &intVal))
2391 return 0;
2392 hData->mIrRate = (uint)intVal;
2393 hasRate = 1;
2395 else if(strcasecmp(ident, "points") == 0)
2397 if (hasPoints) {
2398 TrErrorAt(tr, line, col, "Redefinition of 'points'.\n");
2399 return 0;
2401 if(!TrReadOperator(tr, "="))
2402 return 0;
2403 TrIndication(tr, &line, &col);
2404 if(!TrReadInt(tr, MIN_POINTS, MAX_POINTS, &intVal))
2405 return 0;
2406 points = (uint)intVal;
2407 if(fftSize > 0 && points > fftSize)
2409 TrErrorAt(tr, line, col, "Value exceeds the overridden FFT size.\n");
2410 return 0;
2412 if(points < truncSize)
2414 TrErrorAt(tr, line, col, "Value is below the truncation size.\n");
2415 return 0;
2417 hData->mIrPoints = points;
2418 if(fftSize <= 0)
2420 hData->mFftSize = DEFAULT_FFTSIZE;
2421 hData->mIrSize = 1 + (DEFAULT_FFTSIZE / 2);
2423 else
2425 hData->mFftSize = fftSize;
2426 hData->mIrSize = 1 + (fftSize / 2);
2427 if(points > hData->mIrSize)
2428 hData->mIrSize = points;
2430 hasPoints = 1;
2432 else if(strcasecmp(ident, "azimuths") == 0)
2434 if(hasAzimuths)
2436 TrErrorAt(tr, line, col, "Redefinition of 'azimuths'.\n");
2437 return 0;
2439 if(!TrReadOperator(tr, "="))
2440 return 0;
2441 hData->mIrCount = 0;
2442 hData->mEvCount = 0;
2443 hData->mEvOffset[0] = 0;
2444 for(;;)
2446 if(!TrReadInt(tr, MIN_AZ_COUNT, MAX_AZ_COUNT, &intVal))
2447 return 0;
2448 hData->mAzCount[hData->mEvCount] = (uint)intVal;
2449 hData->mIrCount += (uint)intVal;
2450 hData->mEvCount ++;
2451 if(!TrIsOperator(tr, ","))
2452 break;
2453 if(hData->mEvCount >= MAX_EV_COUNT)
2455 TrError(tr, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT);
2456 return 0;
2458 hData->mEvOffset[hData->mEvCount] = hData->mEvOffset[hData->mEvCount - 1] + ((uint)intVal);
2459 TrReadOperator(tr, ",");
2461 if(hData->mEvCount < MIN_EV_COUNT)
2463 TrErrorAt(tr, line, col, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT);
2464 return 0;
2466 hasAzimuths = 1;
2468 else if(strcasecmp(ident, "radius") == 0)
2470 if(hasRadius)
2472 TrErrorAt(tr, line, col, "Redefinition of 'radius'.\n");
2473 return 0;
2475 if(!TrReadOperator(tr, "="))
2476 return 0;
2477 if(!TrReadFloat(tr, MIN_RADIUS, MAX_RADIUS, &fpVal))
2478 return 0;
2479 hData->mRadius = fpVal;
2480 hasRadius = 1;
2482 else if(strcasecmp(ident, "distance") == 0)
2484 if(hasDistance)
2486 TrErrorAt(tr, line, col, "Redefinition of 'distance'.\n");
2487 return 0;
2489 if(!TrReadOperator(tr, "="))
2490 return 0;
2491 if(!TrReadFloat(tr, MIN_DISTANCE, MAX_DISTANCE, & fpVal))
2492 return 0;
2493 hData->mDistance = fpVal;
2494 hasDistance = 1;
2496 else
2498 TrErrorAt(tr, line, col, "Expected a metric name.\n");
2499 return 0;
2501 TrSkipWhitespace (tr);
2503 return 1;
2506 // Parse an index pair from the data set definition.
2507 static int ReadIndexPair(TokenReaderT *tr, const HrirDataT *hData, uint *ei, uint *ai)
2509 int intVal;
2510 if(!TrReadInt(tr, 0, (int)hData->mEvCount, &intVal))
2511 return 0;
2512 *ei = (uint)intVal;
2513 if(!TrReadOperator(tr, ","))
2514 return 0;
2515 if(!TrReadInt(tr, 0, (int)hData->mAzCount[*ei], &intVal))
2516 return 0;
2517 *ai = (uint)intVal;
2518 return 1;
2521 // Match the source format from a given identifier.
2522 static SourceFormatT MatchSourceFormat(const char *ident)
2524 if(strcasecmp(ident, "wave") == 0)
2525 return SF_WAVE;
2526 if(strcasecmp(ident, "bin_le") == 0)
2527 return SF_BIN_LE;
2528 if(strcasecmp(ident, "bin_be") == 0)
2529 return SF_BIN_BE;
2530 if(strcasecmp(ident, "ascii") == 0)
2531 return SF_ASCII;
2532 return SF_NONE;
2535 // Match the source element type from a given identifier.
2536 static ElementTypeT MatchElementType(const char *ident)
2538 if(strcasecmp(ident, "int") == 0)
2539 return ET_INT;
2540 if(strcasecmp(ident, "fp") == 0)
2541 return ET_FP;
2542 return ET_NONE;
2545 // Parse and validate a source reference from the data set definition.
2546 static int ReadSourceRef(TokenReaderT *tr, SourceRefT *src)
2548 char ident[MAX_IDENT_LEN+1];
2549 uint line, col;
2550 int intVal;
2552 TrIndication(tr, &line, &col);
2553 if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
2554 return 0;
2555 src->mFormat = MatchSourceFormat(ident);
2556 if(src->mFormat == SF_NONE)
2558 TrErrorAt(tr, line, col, "Expected a source format.\n");
2559 return 0;
2561 if(!TrReadOperator(tr, "("))
2562 return 0;
2563 if(src->mFormat == SF_WAVE)
2565 if(!TrReadInt(tr, 0, MAX_WAVE_CHANNELS, &intVal))
2566 return 0;
2567 src->mType = ET_NONE;
2568 src->mSize = 0;
2569 src->mBits = 0;
2570 src->mChannel = (uint)intVal;
2571 src->mSkip = 0;
2573 else
2575 TrIndication(tr, &line, &col);
2576 if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
2577 return 0;
2578 src->mType = MatchElementType(ident);
2579 if(src->mType == ET_NONE)
2581 TrErrorAt(tr, line, col, "Expected a source element type.\n");
2582 return 0;
2584 if(src->mFormat == SF_BIN_LE || src->mFormat == SF_BIN_BE)
2586 if(!TrReadOperator(tr, ","))
2587 return 0;
2588 if(src->mType == ET_INT)
2590 if(!TrReadInt(tr, MIN_BIN_SIZE, MAX_BIN_SIZE, &intVal))
2591 return 0;
2592 src->mSize = (uint)intVal;
2593 if(!TrIsOperator(tr, ","))
2594 src->mBits = (int)(8*src->mSize);
2595 else
2597 TrReadOperator(tr, ",");
2598 TrIndication(tr, &line, &col);
2599 if(!TrReadInt(tr, -2147483647-1, 2147483647, &intVal))
2600 return 0;
2601 if(abs(intVal) < MIN_BIN_BITS || ((uint)abs(intVal)) > (8*src->mSize))
2603 TrErrorAt(tr, line, col, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS, 8*src->mSize);
2604 return 0;
2606 src->mBits = intVal;
2609 else
2611 TrIndication(tr, &line, &col);
2612 if(!TrReadInt(tr, -2147483647-1, 2147483647, &intVal))
2613 return 0;
2614 if(intVal != 4 && intVal != 8)
2616 TrErrorAt(tr, line, col, "Expected a value of 4 or 8.\n");
2617 return 0;
2619 src->mSize = (uint)intVal;
2620 src->mBits = 0;
2623 else if(src->mFormat == SF_ASCII && src->mType == ET_INT)
2625 if(!TrReadOperator(tr, ","))
2626 return 0;
2627 if(!TrReadInt(tr, MIN_ASCII_BITS, MAX_ASCII_BITS, &intVal))
2628 return 0;
2629 src->mSize = 0;
2630 src->mBits = intVal;
2632 else
2634 src->mSize = 0;
2635 src->mBits = 0;
2638 if(!TrIsOperator(tr, ";"))
2639 src->mSkip = 0;
2640 else
2642 TrReadOperator(tr, ";");
2643 if(!TrReadInt (tr, 0, 0x7FFFFFFF, &intVal))
2644 return 0;
2645 src->mSkip = (uint)intVal;
2648 if(!TrReadOperator(tr, ")"))
2649 return 0;
2650 if(TrIsOperator(tr, "@"))
2652 TrReadOperator(tr, "@");
2653 if(!TrReadInt(tr, 0, 0x7FFFFFFF, &intVal))
2654 return 0;
2655 src->mOffset = (uint)intVal;
2657 else
2658 src->mOffset = 0;
2659 if(!TrReadOperator(tr, ":"))
2660 return 0;
2661 if(!TrReadString(tr, MAX_PATH_LEN, src->mPath))
2662 return 0;
2663 return 1;
2666 // Process the list of sources in the data set definition.
2667 static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *hData)
2669 uint *setCount, *setFlag;
2670 uint line, col, ei, ai;
2671 SourceRefT src;
2672 double factor;
2673 double *hrir;
2674 int count;
2676 printf("Loading sources...");
2677 fflush(stdout);
2679 count = 0;
2680 setCount = (uint*)calloc(hData->mEvCount, sizeof(uint));
2681 setFlag = (uint*)calloc(hData->mIrCount, sizeof(uint));
2682 hrir = CreateArray(hData->mIrPoints);
2683 while(TrIsOperator(tr, "["))
2685 TrIndication(tr, & line, & col);
2686 TrReadOperator(tr, "[");
2687 if(!ReadIndexPair(tr, hData, &ei, &ai))
2688 goto error;
2689 if(!TrReadOperator(tr, "]"))
2690 goto error;
2691 if(setFlag[hData->mEvOffset[ei] + ai])
2693 TrErrorAt(tr, line, col, "Redefinition of source.\n");
2694 goto error;
2696 if(!TrReadOperator(tr, "="))
2697 goto error;
2699 factor = 1.0;
2700 for(;;)
2702 if(!ReadSourceRef(tr, &src))
2703 goto error;
2705 // TODO: Would be nice to display 'x of y files', but that would
2706 // require preparing the source refs first to get a total count
2707 // before loading them.
2708 ++count;
2709 printf("\rLoading sources... %d file%s", count, (count==1)?"":"s");
2710 fflush(stdout);
2712 if(!LoadSource(&src, hData->mIrRate, hData->mIrPoints, hrir))
2713 goto error;
2715 if(model == HM_DATASET)
2716 AverageHrirOnset(hrir, 1.0 / factor, ei, ai, hData);
2717 AverageHrirMagnitude(hrir, 1.0 / factor, ei, ai, hData);
2718 factor += 1.0;
2719 if(!TrIsOperator(tr, "+"))
2720 break;
2721 TrReadOperator(tr, "+");
2723 setFlag[hData->mEvOffset[ei] + ai] = 1;
2724 setCount[ei]++;
2726 printf("\n");
2728 ei = 0;
2729 while(ei < hData->mEvCount && setCount[ei] < 1)
2730 ei++;
2731 if(ei < hData->mEvCount)
2733 hData->mEvStart = ei;
2734 while(ei < hData->mEvCount && setCount[ei] == hData->mAzCount[ei])
2735 ei++;
2736 if(ei >= hData->mEvCount)
2738 if(!TrLoad(tr))
2740 DestroyArray(hrir);
2741 free(setFlag);
2742 free(setCount);
2743 return 1;
2745 TrError(tr, "Errant data at end of source list.\n");
2747 else
2748 TrError(tr, "Missing sources for elevation index %d.\n", ei);
2750 else
2751 TrError(tr, "Missing source references.\n");
2753 error:
2754 DestroyArray(hrir);
2755 free(setFlag);
2756 free(setCount);
2757 return 0;
2760 /* Parse the data set definition and process the source data, storing the
2761 * resulting data set as desired. If the input name is NULL it will read
2762 * from standard input.
2764 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 int experimental, const char *outName)
2766 char rateStr[8+1], expName[MAX_PATH_LEN];
2767 TokenReaderT tr;
2768 HrirDataT hData;
2769 FILE *fp;
2770 int ret;
2772 hData.mIrRate = 0;
2773 hData.mSampleType = ST_S24;
2774 hData.mChannelType = CT_LEFTONLY;
2775 hData.mIrPoints = 0;
2776 hData.mFftSize = 0;
2777 hData.mIrSize = 0;
2778 hData.mIrCount = 0;
2779 hData.mEvCount = 0;
2780 hData.mRadius = 0;
2781 hData.mDistance = 0;
2782 fprintf(stdout, "Reading HRIR definition from %s...\n", inName?inName:"stdin");
2783 if(inName != NULL)
2785 fp = fopen(inName, "r");
2786 if(fp == NULL)
2788 fprintf(stderr, "Error: Could not open definition file '%s'\n", inName);
2789 return 0;
2791 TrSetup(fp, inName, &tr);
2793 else
2795 fp = stdin;
2796 TrSetup(fp, "<stdin>", &tr);
2798 if(!ProcessMetrics(&tr, fftSize, truncSize, &hData))
2800 if(inName != NULL)
2801 fclose(fp);
2802 return 0;
2804 hData.mHrirs = CreateArray(hData.mIrCount * hData.mIrSize);
2805 hData.mHrtds = CreateArray(hData.mIrCount);
2806 if(!ProcessSources(model, &tr, &hData))
2808 DestroyArray(hData.mHrtds);
2809 DestroyArray(hData.mHrirs);
2810 if(inName != NULL)
2811 fclose(fp);
2812 return 0;
2814 if(fp != stdin)
2815 fclose(fp);
2816 if(equalize)
2818 double *dfa = CreateArray(1 + (hData.mFftSize/2));
2819 fprintf(stdout, "Calculating diffuse-field average...\n");
2820 CalculateDiffuseFieldAverage(&hData, surface, limit, dfa);
2821 fprintf(stdout, "Performing diffuse-field equalization...\n");
2822 DiffuseFieldEqualize(dfa, &hData);
2823 DestroyArray(dfa);
2825 fprintf(stdout, "Performing minimum phase reconstruction...\n");
2826 ReconstructHrirs(&hData);
2827 if(outRate != 0 && outRate != hData.mIrRate)
2829 fprintf(stdout, "Resampling HRIRs...\n");
2830 ResampleHrirs(outRate, &hData);
2832 fprintf(stdout, "Truncating minimum-phase HRIRs...\n");
2833 hData.mIrPoints = truncSize;
2834 fprintf(stdout, "Synthesizing missing elevations...\n");
2835 if(model == HM_DATASET)
2836 SynthesizeOnsets(&hData);
2837 SynthesizeHrirs(&hData);
2838 fprintf(stdout, "Normalizing final HRIRs...\n");
2839 NormalizeHrirs(&hData);
2840 fprintf(stdout, "Calculating impulse delays...\n");
2841 CalculateHrtds(model, (radius > DEFAULT_CUSTOM_RADIUS) ? radius : hData.mRadius, &hData);
2842 snprintf(rateStr, 8, "%u", hData.mIrRate);
2843 StrSubst(outName, "%r", rateStr, MAX_PATH_LEN, expName);
2844 fprintf(stdout, "Creating MHR data set %s...\n", expName);
2845 ret = StoreMhr(&hData, experimental, expName);
2847 DestroyArray(hData.mHrtds);
2848 DestroyArray(hData.mHrirs);
2849 return ret;
2852 static void PrintHelp(const char *argv0, FILE *ofile)
2854 fprintf(ofile, "Usage: %s <command> [<option>...]\n\n", argv0);
2855 fprintf(ofile, "Options:\n");
2856 fprintf(ofile, " -m Ignored for compatibility.\n");
2857 fprintf(ofile, " -r <rate> Change the data set sample rate to the specified value and\n");
2858 fprintf(ofile, " resample the HRIRs accordingly.\n");
2859 fprintf(ofile, " -f <points> Override the FFT window size (default: %u).\n", DEFAULT_FFTSIZE);
2860 fprintf(ofile, " -e {on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE ? "on" : "off"));
2861 fprintf(ofile, " -s {on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE ? "on" : "off"));
2862 fprintf(ofile, " -l {<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
2863 fprintf(ofile, " average (default: %.2f).\n", DEFAULT_LIMIT);
2864 fprintf(ofile, " -w <points> Specify the size of the truncation window that's applied\n");
2865 fprintf(ofile, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE);
2866 fprintf(ofile, " -d {dataset| Specify the model used for calculating the head-delay timing\n");
2867 fprintf(ofile, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL == HM_DATASET) ? "dataset" : "sphere"));
2868 fprintf(ofile, " -c <size> Use a customized head radius measured ear-to-ear in meters.\n");
2869 fprintf(ofile, " -i <filename> Specify an HRIR definition file to use (defaults to stdin).\n");
2870 fprintf(ofile, " -o <filename> Specify an output file. Overrides command-selected default.\n");
2871 fprintf(ofile, " Use of '%%r' will be substituted with the data set sample rate.\n");
2874 // Standard command line dispatch.
2875 int main(int argc, char *argv[])
2877 const char *inName = NULL, *outName = NULL;
2878 uint outRate, fftSize;
2879 int equalize, surface;
2880 int experimental;
2881 char *end = NULL;
2882 HeadModelT model;
2883 uint truncSize;
2884 double radius;
2885 double limit;
2886 int opt;
2888 if(argc < 2)
2890 fprintf(stdout, "HRTF Processing and Composition Utility\n\n");
2891 PrintHelp(argv[0], stdout);
2892 exit(EXIT_SUCCESS);
2895 outName = "./oalsoft_hrtf_%r.mhr";
2896 outRate = 0;
2897 fftSize = 0;
2898 equalize = DEFAULT_EQUALIZE;
2899 surface = DEFAULT_SURFACE;
2900 limit = DEFAULT_LIMIT;
2901 truncSize = DEFAULT_TRUNCSIZE;
2902 model = DEFAULT_HEAD_MODEL;
2903 radius = DEFAULT_CUSTOM_RADIUS;
2904 experimental = 0;
2906 while((opt=getopt(argc, argv, "mr:f:e:s:l:w:d:c:e:i:o:xh")) != -1)
2908 switch(opt)
2910 case 'm':
2911 fprintf(stderr, "Ignoring unused command '-m'.\n");
2912 break;
2914 case 'r':
2915 outRate = strtoul(optarg, &end, 10);
2916 if(end[0] != '\0' || outRate < MIN_RATE || outRate > MAX_RATE)
2918 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected between %u to %u.\n", optarg, opt, MIN_RATE, MAX_RATE);
2919 exit(EXIT_FAILURE);
2921 break;
2923 case 'f':
2924 fftSize = strtoul(optarg, &end, 10);
2925 if(end[0] != '\0' || (fftSize&(fftSize-1)) || fftSize < MIN_FFTSIZE || fftSize > MAX_FFTSIZE)
2927 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected a power-of-two between %u to %u.\n", optarg, opt, MIN_FFTSIZE, MAX_FFTSIZE);
2928 exit(EXIT_FAILURE);
2930 break;
2932 case 'e':
2933 if(strcmp(optarg, "on") == 0)
2934 equalize = 1;
2935 else if(strcmp(optarg, "off") == 0)
2936 equalize = 0;
2937 else
2939 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg, opt);
2940 exit(EXIT_FAILURE);
2942 break;
2944 case 's':
2945 if(strcmp(optarg, "on") == 0)
2946 surface = 1;
2947 else if(strcmp(optarg, "off") == 0)
2948 surface = 0;
2949 else
2951 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg, opt);
2952 exit(EXIT_FAILURE);
2954 break;
2956 case 'l':
2957 if(strcmp(optarg, "none") == 0)
2958 limit = 0.0;
2959 else
2961 limit = strtod(optarg, &end);
2962 if(end[0] != '\0' || limit < MIN_LIMIT || limit > MAX_LIMIT)
2964 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected between %.0f to %.0f.\n", optarg, opt, MIN_LIMIT, MAX_LIMIT);
2965 exit(EXIT_FAILURE);
2968 break;
2970 case 'w':
2971 truncSize = strtoul(optarg, &end, 10);
2972 if(end[0] != '\0' || truncSize < MIN_TRUNCSIZE || truncSize > MAX_TRUNCSIZE || (truncSize%MOD_TRUNCSIZE))
2974 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected multiple of %u between %u to %u.\n", optarg, opt, MOD_TRUNCSIZE, MIN_TRUNCSIZE, MAX_TRUNCSIZE);
2975 exit(EXIT_FAILURE);
2977 break;
2979 case 'd':
2980 if(strcmp(optarg, "dataset") == 0)
2981 model = HM_DATASET;
2982 else if(strcmp(optarg, "sphere") == 0)
2983 model = HM_SPHERE;
2984 else
2986 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected dataset or sphere.\n", optarg, opt);
2987 exit(EXIT_FAILURE);
2989 break;
2991 case 'c':
2992 radius = strtod(optarg, &end);
2993 if(end[0] != '\0' || radius < MIN_CUSTOM_RADIUS || radius > MAX_CUSTOM_RADIUS)
2995 fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected between %.2f to %.2f.\n", optarg, opt, MIN_CUSTOM_RADIUS, MAX_CUSTOM_RADIUS);
2996 exit(EXIT_FAILURE);
2998 break;
3000 case 'i':
3001 inName = optarg;
3002 break;
3004 case 'o':
3005 outName = optarg;
3006 break;
3008 case 'x':
3009 experimental = 1;
3010 break;
3012 case 'h':
3013 PrintHelp(argv[0], stdout);
3014 exit(EXIT_SUCCESS);
3016 default: /* '?' */
3017 PrintHelp(argv[0], stderr);
3018 exit(EXIT_FAILURE);
3022 if(!ProcessDefinition(inName, outRate, fftSize, equalize, surface, limit,
3023 truncSize, model, radius, experimental, outName))
3024 return -1;
3025 fprintf(stdout, "Operation completed.\n");
3027 return EXIT_SUCCESS;