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
88 #include "win_main_utf8.h"
91 #define M_PI (3.14159265358979323846)
95 // The epsilon used to maintain signal stability.
96 #define EPSILON (1e-9)
98 // Constants for accessing the token reader's ring buffer.
99 #define TR_RING_BITS (16)
100 #define TR_RING_SIZE (1 << TR_RING_BITS)
101 #define TR_RING_MASK (TR_RING_SIZE - 1)
103 // The token reader's load interval in bytes.
104 #define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
106 // The maximum identifier length used when processing the data set
108 #define MAX_IDENT_LEN (16)
110 // The maximum path length used when processing filenames.
111 #define MAX_PATH_LEN (256)
113 // The limits for the sample 'rate' metric in the data set definition and for
115 #define MIN_RATE (32000)
116 #define MAX_RATE (96000)
118 // The limits for the HRIR 'points' metric in the data set definition.
119 #define MIN_POINTS (16)
120 #define MAX_POINTS (8192)
122 // The limit to the number of 'distances' listed in the data set definition.
123 #define MAX_FD_COUNT (16)
125 // The limits to the number of 'azimuths' listed in the data set definition.
126 #define MIN_EV_COUNT (5)
127 #define MAX_EV_COUNT (128)
129 // The limits for each of the 'azimuths' listed in the data set definition.
130 #define MIN_AZ_COUNT (1)
131 #define MAX_AZ_COUNT (128)
133 // The limits for the listener's head 'radius' in the data set definition.
134 #define MIN_RADIUS (0.05)
135 #define MAX_RADIUS (0.15)
137 // The limits for the 'distance' from source to listener for each field in
138 // the definition file.
139 #define MIN_DISTANCE (0.05)
140 #define MAX_DISTANCE (2.50)
142 // The maximum number of channels that can be addressed for a WAVE file
143 // source listed in the data set definition.
144 #define MAX_WAVE_CHANNELS (65535)
146 // The limits to the byte size for a binary source listed in the definition
148 #define MIN_BIN_SIZE (2)
149 #define MAX_BIN_SIZE (4)
151 // The minimum number of significant bits for binary sources listed in the
152 // data set definition. The maximum is calculated from the byte size.
153 #define MIN_BIN_BITS (16)
155 // The limits to the number of significant bits for an ASCII source listed in
156 // the data set definition.
157 #define MIN_ASCII_BITS (16)
158 #define MAX_ASCII_BITS (32)
160 // The limits to the FFT window size override on the command line.
161 #define MIN_FFTSIZE (65536)
162 #define MAX_FFTSIZE (131072)
164 // The limits to the equalization range limit on the command line.
165 #define MIN_LIMIT (2.0)
166 #define MAX_LIMIT (120.0)
168 // The limits to the truncation window size on the command line.
169 #define MIN_TRUNCSIZE (16)
170 #define MAX_TRUNCSIZE (512)
172 // The limits to the custom head radius on the command line.
173 #define MIN_CUSTOM_RADIUS (0.05)
174 #define MAX_CUSTOM_RADIUS (0.15)
176 // The truncation window size must be a multiple of the below value to allow
177 // for vectorized convolution.
178 #define MOD_TRUNCSIZE (8)
180 // The defaults for the command line options.
181 #define DEFAULT_FFTSIZE (65536)
182 #define DEFAULT_EQUALIZE (1)
183 #define DEFAULT_SURFACE (1)
184 #define DEFAULT_LIMIT (24.0)
185 #define DEFAULT_TRUNCSIZE (32)
186 #define DEFAULT_HEAD_MODEL (HM_DATASET)
187 #define DEFAULT_CUSTOM_RADIUS (0.0)
189 // The four-character-codes for RIFF/RIFX WAVE file chunks.
190 #define FOURCC_RIFF (0x46464952) // 'RIFF'
191 #define FOURCC_RIFX (0x58464952) // 'RIFX'
192 #define FOURCC_WAVE (0x45564157) // 'WAVE'
193 #define FOURCC_FMT (0x20746D66) // 'fmt '
194 #define FOURCC_DATA (0x61746164) // 'data'
195 #define FOURCC_LIST (0x5453494C) // 'LIST'
196 #define FOURCC_WAVL (0x6C766177) // 'wavl'
197 #define FOURCC_SLNT (0x746E6C73) // 'slnt'
199 // The supported wave formats.
200 #define WAVE_FORMAT_PCM (0x0001)
201 #define WAVE_FORMAT_IEEE_FLOAT (0x0003)
202 #define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
204 // The maximum propagation delay value supported by OpenAL Soft.
205 #define MAX_HRTD (63.0)
207 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
208 // response protocol 02.
209 #define MHR_FORMAT ("MinPHR02")
211 // Sample and channel type enum values.
217 // Certain iterations rely on these integer enum values.
224 // Byte order for the serialization routines.
231 // Source format for the references listed in the data set definition.
234 SF_WAVE
, // RIFF/RIFX WAVE file.
235 SF_BIN_LE
, // Little-endian binary file.
236 SF_BIN_BE
, // Big-endian binary file.
237 SF_ASCII
// ASCII text file.
240 // Element types for the references listed in the data set definition.
243 ET_INT
, // Integer elements.
244 ET_FP
// Floating-point elements.
247 // Head model used for calculating the impulse delays.
250 HM_DATASET
, // Measure the onset from the dataset.
251 HM_SPHERE
// Calculate the onset using a spherical head model.
254 /* Unsigned integer type. */
255 using uint
= unsigned int;
257 /* Complex double type. */
258 using complex_d
= std::complex<double>;
261 // Token reader state for parsing the data set definition.
262 struct TokenReaderT
{
267 char mRing
[TR_RING_SIZE
];
272 // Source reference state used when loading sources.
274 SourceFormatT mFormat
;
281 char mPath
[MAX_PATH_LEN
+1];
284 // Structured HRIR storage for stereo azimuth pairs, elevations, and fields.
286 double mAzimuth
{0.0};
288 double mDelays
[2]{0.0, 0.0};
289 double *mIrs
[2]{nullptr, nullptr};
293 double mElevation
{0.0};
296 HrirAzT
*mAzs
{nullptr};
300 double mDistance
{0.0};
304 HrirEvT
*mEvs
{nullptr};
307 // The HRIR metrics and data set used when loading, processing, and storing
308 // the resulting HRTF.
311 SampleTypeT mSampleType
{ST_S24
};
312 ChannelTypeT mChannelType
{CT_NONE
};
320 std::vector
<double> mHrirsBase
;
321 std::vector
<HrirEvT
> mEvsBase
;
322 std::vector
<HrirAzT
> mAzsBase
;
324 std::vector
<HrirFdT
> mFds
;
327 // The resampler metrics and FIR filter.
330 std::vector
<double> mF
;
334 /*****************************
335 *** Token reader routines ***
336 *****************************/
338 /* Whitespace is not significant. It can process tokens as identifiers, numbers
339 * (integer and floating-point), strings, and operators. Strings must be
340 * encapsulated by double-quotes and cannot span multiple lines.
343 // Setup the reader on the given file. The filename can be NULL if no error
344 // output is desired.
345 static void TrSetup(FILE *fp
, const char *filename
, TokenReaderT
*tr
)
347 const char *name
= nullptr;
351 const char *slash
= strrchr(filename
, '/');
354 const char *bslash
= strrchr(slash
+1, '\\');
355 if(bslash
) name
= bslash
+1;
360 const char *bslash
= strrchr(filename
, '\\');
361 if(bslash
) name
= bslash
+1;
362 else name
= filename
;
374 // Prime the reader's ring buffer, and return a result indicating that there
375 // is text to process.
376 static int TrLoad(TokenReaderT
*tr
)
378 size_t toLoad
, in
, count
;
380 toLoad
= TR_RING_SIZE
- (tr
->mIn
- tr
->mOut
);
381 if(toLoad
>= TR_LOAD_SIZE
&& !feof(tr
->mFile
))
383 // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
384 toLoad
= TR_LOAD_SIZE
;
385 in
= tr
->mIn
&TR_RING_MASK
;
386 count
= TR_RING_SIZE
- in
;
389 tr
->mIn
+= fread(&tr
->mRing
[in
], 1, count
, tr
->mFile
);
390 tr
->mIn
+= fread(&tr
->mRing
[0], 1, toLoad
-count
, tr
->mFile
);
393 tr
->mIn
+= fread(&tr
->mRing
[in
], 1, toLoad
, tr
->mFile
);
395 if(tr
->mOut
>= TR_RING_SIZE
)
397 tr
->mOut
-= TR_RING_SIZE
;
398 tr
->mIn
-= TR_RING_SIZE
;
401 if(tr
->mIn
> tr
->mOut
)
406 // Error display routine. Only displays when the base name is not NULL.
407 static void TrErrorVA(const TokenReaderT
*tr
, uint line
, uint column
, const char *format
, va_list argPtr
)
411 fprintf(stderr
, "Error (%s:%u:%u): ", tr
->mName
, line
, column
);
412 vfprintf(stderr
, format
, argPtr
);
415 // Used to display an error at a saved line/column.
416 static void TrErrorAt(const TokenReaderT
*tr
, uint line
, uint column
, const char *format
, ...)
420 va_start(argPtr
, format
);
421 TrErrorVA(tr
, line
, column
, format
, argPtr
);
425 // Used to display an error at the current line/column.
426 static void TrError(const TokenReaderT
*tr
, const char *format
, ...)
430 va_start(argPtr
, format
);
431 TrErrorVA(tr
, tr
->mLine
, tr
->mColumn
, format
, argPtr
);
435 // Skips to the next line.
436 static void TrSkipLine(TokenReaderT
*tr
)
442 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
454 // Skips to the next token.
455 static int TrSkipWhitespace(TokenReaderT
*tr
)
459 char ch
{tr
->mRing
[tr
->mOut
&TR_RING_MASK
]};
479 // Get the line and/or column of the next token (or the end of input).
480 static void TrIndication(TokenReaderT
*tr
, uint
*line
, uint
*column
)
482 TrSkipWhitespace(tr
);
483 if(line
) *line
= tr
->mLine
;
484 if(column
) *column
= tr
->mColumn
;
487 // Checks to see if a token is (likely to be) an identifier. It does not
488 // display any errors and will not proceed to the next token.
489 static int TrIsIdent(TokenReaderT
*tr
)
491 if(!TrSkipWhitespace(tr
))
493 char ch
{tr
->mRing
[tr
->mOut
&TR_RING_MASK
]};
494 return ch
== '_' || isalpha(ch
);
498 // Checks to see if a token is the given operator. It does not display any
499 // errors and will not proceed to the next token.
500 static int TrIsOperator(TokenReaderT
*tr
, const char *op
)
505 if(!TrSkipWhitespace(tr
))
509 while(op
[len
] != '\0' && out
< tr
->mIn
)
511 ch
= tr
->mRing
[out
&TR_RING_MASK
];
512 if(ch
!= op
[len
]) break;
521 /* The TrRead*() routines obtain the value of a matching token type. They
522 * display type, form, and boundary errors and will proceed to the next
526 // Reads and validates an identifier token.
527 static int TrReadIdent(TokenReaderT
*tr
, const uint maxLen
, char *ident
)
533 if(TrSkipWhitespace(tr
))
536 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
537 if(ch
== '_' || isalpha(ch
))
547 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
548 } while(ch
== '_' || isdigit(ch
) || isalpha(ch
));
556 TrErrorAt(tr
, tr
->mLine
, col
, "Identifier is too long.\n");
560 TrErrorAt(tr
, tr
->mLine
, col
, "Expected an identifier.\n");
564 // Reads and validates (including bounds) an integer token.
565 static int TrReadInt(TokenReaderT
*tr
, const int loBound
, const int hiBound
, int *value
)
567 uint col
, digis
, len
;
571 if(TrSkipWhitespace(tr
))
575 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
576 if(ch
== '+' || ch
== '-')
585 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
586 if(!isdigit(ch
)) break;
594 if(digis
> 0 && ch
!= '.' && !isalpha(ch
))
598 TrErrorAt(tr
, tr
->mLine
, col
, "Integer is too long.");
602 *value
= strtol(temp
, nullptr, 10);
603 if(*value
< loBound
|| *value
> hiBound
)
605 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a value from %d to %d.\n", loBound
, hiBound
);
611 TrErrorAt(tr
, tr
->mLine
, col
, "Expected an integer.\n");
615 // Reads and validates (including bounds) a float token.
616 static int TrReadFloat(TokenReaderT
*tr
, const double loBound
, const double hiBound
, double *value
)
618 uint col
, digis
, len
;
622 if(TrSkipWhitespace(tr
))
626 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
627 if(ch
== '+' || ch
== '-')
637 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
638 if(!isdigit(ch
)) break;
654 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
655 if(!isdigit(ch
)) break;
664 if(ch
== 'E' || ch
== 'e')
671 if(ch
== '+' || ch
== '-')
680 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
681 if(!isdigit(ch
)) break;
690 if(digis
> 0 && ch
!= '.' && !isalpha(ch
))
694 TrErrorAt(tr
, tr
->mLine
, col
, "Float is too long.");
698 *value
= strtod(temp
, nullptr);
699 if(*value
< loBound
|| *value
> hiBound
)
701 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a value from %f to %f.\n", loBound
, hiBound
);
710 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a float.\n");
714 // Reads and validates a string token.
715 static int TrReadString(TokenReaderT
*tr
, const uint maxLen
, char *text
)
721 if(TrSkipWhitespace(tr
))
724 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
731 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
737 TrErrorAt(tr
, tr
->mLine
, col
, "Unterminated string at end of line.\n");
746 tr
->mColumn
+= 1 + len
;
747 TrErrorAt(tr
, tr
->mLine
, col
, "Unterminated string at end of input.\n");
750 tr
->mColumn
+= 2 + len
;
753 TrErrorAt(tr
, tr
->mLine
, col
, "String is too long.\n");
760 TrErrorAt(tr
, tr
->mLine
, col
, "Expected a string.\n");
764 // Reads and validates the given operator.
765 static int TrReadOperator(TokenReaderT
*tr
, const char *op
)
771 if(TrSkipWhitespace(tr
))
775 while(op
[len
] != '\0' && TrLoad(tr
))
777 ch
= tr
->mRing
[tr
->mOut
&TR_RING_MASK
];
778 if(ch
!= op
[len
]) break;
786 TrErrorAt(tr
, tr
->mLine
, col
, "Expected '%s' operator.\n", op
);
790 /* Performs a string substitution. Any case-insensitive occurrences of the
791 * pattern string are replaced with the replacement string. The result is
792 * truncated if necessary.
794 static int StrSubst(const char *in
, const char *pat
, const char *rep
, const size_t maxLen
, char *out
)
796 size_t inLen
, patLen
, repLen
;
801 patLen
= strlen(pat
);
802 repLen
= strlen(rep
);
806 while(si
< inLen
&& di
< maxLen
)
808 if(patLen
<= inLen
-si
)
810 if(strncasecmp(&in
[si
], pat
, patLen
) == 0)
812 if(repLen
> maxLen
-di
)
814 repLen
= maxLen
- di
;
817 strncpy(&out
[di
], rep
, repLen
);
833 /*********************
834 *** Math routines ***
835 *********************/
837 // Simple clamp routine.
838 static double Clamp(const double val
, const double lower
, const double upper
)
840 return std::min(std::max(val
, lower
), upper
);
843 // Performs linear interpolation.
844 static double Lerp(const double a
, const double b
, const double f
)
846 return a
+ f
* (b
- a
);
849 static inline uint
dither_rng(uint
*seed
)
851 *seed
= *seed
* 96314165 + 907633515;
855 // Performs a triangular probability density function dither. The input samples
856 // should be normalized (-1 to +1).
857 static void TpdfDither(double *RESTRICT out
, const double *RESTRICT in
, const double scale
,
858 const int count
, const int step
, uint
*seed
)
860 static constexpr double PRNG_SCALE
= 1.0 / UINT_MAX
;
862 for(int i
{0};i
< count
;i
++)
864 uint prn0
{dither_rng(seed
)};
865 uint prn1
{dither_rng(seed
)};
866 out
[i
*step
] = std::round(in
[i
]*scale
+ (prn0
*PRNG_SCALE
- prn1
*PRNG_SCALE
));
870 /* Fast Fourier transform routines. The number of points must be a power of
874 // Performs bit-reversal ordering.
875 static void FftArrange(const uint n
, complex_d
*inout
)
879 // Handle in-place arrangement.
884 std::swap(inout
[rk
], inout
[k
]);
893 // Performs the summation.
894 static void FftSummation(const int n
, const double s
, complex_d
*cplx
)
901 for(m
= 1, m2
= 2;m
< n
; m
<<= 1, m2
<<= 1)
903 // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
904 double sm
= sin(0.5 * pi
/ m
);
905 auto v
= complex_d
{-2.0*sm
*sm
, -sin(pi
/ m
)};
906 auto w
= complex_d
{1.0, 0.0};
909 for(k
= i
;k
< n
;k
+= m2
)
912 auto t
= w
* cplx
[mk
];
913 cplx
[mk
] = cplx
[k
] - t
;
914 cplx
[k
] = cplx
[k
] + t
;
921 // Performs a forward FFT.
922 static void FftForward(const uint n
, complex_d
*inout
)
924 FftArrange(n
, inout
);
925 FftSummation(n
, 1.0, inout
);
928 // Performs an inverse FFT.
929 static void FftInverse(const uint n
, complex_d
*inout
)
931 FftArrange(n
, inout
);
932 FftSummation(n
, -1.0, inout
);
934 for(uint i
{0};i
< n
;i
++)
938 /* Calculate the complex helical sequence (or discrete-time analytical signal)
939 * of the given input using the Hilbert transform. Given the natural logarithm
940 * of a signal's magnitude response, the imaginary components can be used as
941 * the angles for minimum-phase reconstruction.
943 static void Hilbert(const uint n
, complex_d
*inout
)
947 // Handle in-place operation.
951 FftInverse(n
, inout
);
952 for(i
= 1;i
< (n
+1)/2;i
++)
954 /* Increment i if n is even. */
957 inout
[i
] = complex_d
{0.0, 0.0};
958 FftForward(n
, inout
);
961 /* Calculate the magnitude response of the given input. This is used in
962 * place of phase decomposition, since the phase residuals are discarded for
963 * minimum phase reconstruction. The mirrored half of the response is also
966 static void MagnitudeResponse(const uint n
, const complex_d
*in
, double *out
)
968 const uint m
= 1 + (n
/ 2);
971 out
[i
] = std::max(std::abs(in
[i
]), EPSILON
);
974 /* Apply a range limit (in dB) to the given magnitude response. This is used
975 * to adjust the effects of the diffuse-field average on the equalization
978 static void LimitMagnitudeResponse(const uint n
, const uint m
, const double limit
, const double *in
, double *out
)
981 uint i
, lower
, upper
;
984 halfLim
= limit
/ 2.0;
985 // Convert the response to dB.
987 out
[i
] = 20.0 * std::log10(in
[i
]);
988 // Use six octaves to calculate the average magnitude of the signal.
989 lower
= (static_cast<uint
>(std::ceil(n
/ std::pow(2.0, 8.0)))) - 1;
990 upper
= (static_cast<uint
>(std::floor(n
/ std::pow(2.0, 2.0)))) - 1;
992 for(i
= lower
;i
<= upper
;i
++)
994 ave
/= upper
- lower
+ 1;
995 // Keep the response within range of the average magnitude.
997 out
[i
] = Clamp(out
[i
], ave
- halfLim
, ave
+ halfLim
);
998 // Convert the response back to linear magnitude.
1000 out
[i
] = std::pow(10.0, out
[i
] / 20.0);
1003 /* Reconstructs the minimum-phase component for the given magnitude response
1004 * of a signal. This is equivalent to phase recomposition, sans the missing
1005 * residuals (which were discarded). The mirrored half of the response is
1008 static void MinimumPhase(const uint n
, const double *in
, complex_d
*out
)
1010 const uint m
= 1 + (n
/ 2);
1011 std::vector
<double> mags(n
);
1014 for(i
= 0;i
< m
;i
++)
1016 mags
[i
] = std::max(EPSILON
, in
[i
]);
1017 out
[i
] = complex_d
{std::log(mags
[i
]), 0.0};
1021 mags
[i
] = mags
[n
- i
];
1022 out
[i
] = out
[n
- i
];
1025 // Remove any DC offset the filter has.
1027 for(i
= 0;i
< n
;i
++)
1029 auto a
= std::exp(complex_d
{0.0, out
[i
].imag()});
1030 out
[i
] = complex_d
{mags
[i
], 0.0} * a
;
1035 /***************************
1036 *** Resampler functions ***
1037 ***************************/
1039 /* This is the normalized cardinal sine (sinc) function.
1041 * sinc(x) = { 1, x = 0
1042 * { sin(pi x) / (pi x), otherwise.
1044 static double Sinc(const double x
)
1046 if(std::abs(x
) < EPSILON
)
1048 return std::sin(M_PI
* x
) / (M_PI
* x
);
1051 /* The zero-order modified Bessel function of the first kind, used for the
1054 * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
1055 * = sum_{k=0}^inf ((x / 2)^k / k!)^2
1057 static double BesselI_0(const double x
)
1059 double term
, sum
, x2
, y
, last_sum
;
1062 // Start at k=1 since k=0 is trivial.
1068 // Let the integration converge until the term of the sum is no longer
1076 } while(sum
!= last_sum
);
1080 /* Calculate a Kaiser window from the given beta value and a normalized k
1083 * w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
1086 * Where k can be calculated as:
1088 * k = i / l, where -l <= i <= l.
1092 * k = 2 i / M - 1, where 0 <= i <= M.
1094 static double Kaiser(const double b
, const double k
)
1096 if(!(k
>= -1.0 && k
<= 1.0))
1098 return BesselI_0(b
* std::sqrt(1.0 - k
*k
)) / BesselI_0(b
);
1101 // Calculates the greatest common divisor of a and b.
1102 static uint
Gcd(uint x
, uint y
)
1113 /* Calculates the size (order) of the Kaiser window. Rejection is in dB and
1114 * the transition width is normalized frequency (0.5 is nyquist).
1116 * M = { ceil((r - 7.95) / (2.285 2 pi f_t)), r > 21
1117 * { ceil(5.79 / 2 pi f_t), r <= 21.
1120 static uint
CalcKaiserOrder(const double rejection
, const double transition
)
1122 double w_t
= 2.0 * M_PI
* transition
;
1123 if(rejection
> 21.0)
1124 return static_cast<uint
>(std::ceil((rejection
- 7.95) / (2.285 * w_t
)));
1125 return static_cast<uint
>(std::ceil(5.79 / w_t
));
1128 // Calculates the beta value of the Kaiser window. Rejection is in dB.
1129 static double CalcKaiserBeta(const double rejection
)
1131 if(rejection
> 50.0)
1132 return 0.1102 * (rejection
- 8.7);
1133 if(rejection
>= 21.0)
1134 return (0.5842 * std::pow(rejection
- 21.0, 0.4)) +
1135 (0.07886 * (rejection
- 21.0));
1139 /* Calculates a point on the Kaiser-windowed sinc filter for the given half-
1140 * width, beta, gain, and cutoff. The point is specified in non-normalized
1141 * samples, from 0 to M, where M = (2 l + 1).
1143 * w(k) 2 p f_t sinc(2 f_t x)
1145 * x -- centered sample index (i - l)
1146 * k -- normalized and centered window index (x / l)
1147 * w(k) -- window function (Kaiser)
1148 * p -- gain compensation factor when sampling
1149 * f_t -- normalized center frequency (or cutoff; 0.5 is nyquist)
1151 static double SincFilter(const int l
, const double b
, const double gain
, const double cutoff
, const int i
)
1153 return Kaiser(b
, static_cast<double>(i
- l
) / l
) * 2.0 * gain
* cutoff
* Sinc(2.0 * cutoff
* (i
- l
));
1156 /* This is a polyphase sinc-filtered resampler.
1158 * Upsample Downsample
1160 * p/q = 3/2 p/q = 3/5
1162 * M-+-+-+-> M-+-+-+->
1163 * -------------------+ ---------------------+
1164 * p s * f f f f|f| | p s * f f f f f |
1165 * | 0 * 0 0 0|0|0 | | 0 * 0 0 0 0|0| |
1166 * v 0 * 0 0|0|0 0 | v 0 * 0 0 0|0|0 |
1167 * s * f|f|f f f | s * f f|f|f f |
1168 * 0 * |0|0 0 0 0 | 0 * 0|0|0 0 0 |
1169 * --------+=+--------+ 0 * |0|0 0 0 0 |
1170 * d . d .|d|. d . d ----------+=+--------+
1171 * d . . . .|d|. . . .
1175 * P_f(i,j) = q i mod p + pj
1176 * P_s(i,j) = floor(q i / p) - j
1177 * d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
1178 * { f[P_f(i,j)] s[P_s(i,j)], P_f(i,j) < M
1179 * { 0, P_f(i,j) >= M. }
1182 // Calculate the resampling metrics and build the Kaiser-windowed sinc filter
1183 // that's used to cut frequencies above the destination nyquist.
1184 static void ResamplerSetup(ResamplerT
*rs
, const uint srcRate
, const uint dstRate
)
1186 double cutoff
, width
, beta
;
1190 gcd
= Gcd(srcRate
, dstRate
);
1191 rs
->mP
= dstRate
/ gcd
;
1192 rs
->mQ
= srcRate
/ gcd
;
1193 /* The cutoff is adjusted by half the transition width, so the transition
1194 * ends before the nyquist (0.5). Both are scaled by the downsampling
1199 cutoff
= 0.475 / rs
->mP
;
1200 width
= 0.05 / rs
->mP
;
1204 cutoff
= 0.475 / rs
->mQ
;
1205 width
= 0.05 / rs
->mQ
;
1207 // A rejection of -180 dB is used for the stop band. Round up when
1208 // calculating the left offset to avoid increasing the transition width.
1209 l
= (CalcKaiserOrder(180.0, width
)+1) / 2;
1210 beta
= CalcKaiserBeta(180.0);
1213 rs
->mF
.resize(rs
->mM
);
1214 for(i
= 0;i
< (static_cast<int>(rs
->mM
));i
++)
1215 rs
->mF
[i
] = SincFilter(static_cast<int>(l
), beta
, rs
->mP
, cutoff
, i
);
1218 // Perform the upsample-filter-downsample resampling operation using a
1219 // polyphase filter implementation.
1220 static void ResamplerRun(ResamplerT
*rs
, const uint inN
, const double *in
, const uint outN
, double *out
)
1222 const uint p
= rs
->mP
, q
= rs
->mQ
, m
= rs
->mM
, l
= rs
->mL
;
1223 std::vector
<double> workspace
;
1224 const double *f
= rs
->mF
.data();
1232 // Handle in-place operation.
1235 workspace
.resize(outN
);
1236 work
= workspace
.data();
1240 // Resample the input.
1241 for(i
= 0;i
< outN
;i
++)
1244 // Input starts at l to compensate for the filter delay. This will
1245 // drop any build-up from the first half of the filter.
1246 j_f
= (l
+ (q
* i
)) % p
;
1247 j_s
= (l
+ (q
* i
)) / p
;
1250 // Only take input when 0 <= j_s < inN. This single unsigned
1251 // comparison catches both cases.
1253 r
+= f
[j_f
] * in
[j_s
];
1259 // Clean up after in-place operation.
1262 for(i
= 0;i
< outN
;i
++)
1267 /*************************
1268 *** File source input ***
1269 *************************/
1271 // Read a binary value of the specified byte order and byte size from a file,
1272 // storing it as a 32-bit unsigned integer.
1273 static int ReadBin4(FILE *fp
, const char *filename
, const ByteOrderT order
, const uint bytes
, uint32_t *out
)
1279 if(fread(in
, 1, bytes
, fp
) != bytes
)
1281 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1288 for(i
= 0;i
< bytes
;i
++)
1289 accum
= (accum
<<8) | in
[bytes
- i
- 1];
1292 for(i
= 0;i
< bytes
;i
++)
1293 accum
= (accum
<<8) | in
[i
];
1302 // Read a binary value of the specified byte order from a file, storing it as
1303 // a 64-bit unsigned integer.
1304 static int ReadBin8(FILE *fp
, const char *filename
, const ByteOrderT order
, uint64_t *out
)
1310 if(fread(in
, 1, 8, fp
) != 8)
1312 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1319 for(i
= 0;i
< 8;i
++)
1320 accum
= (accum
<<8) | in
[8 - i
- 1];
1323 for(i
= 0;i
< 8;i
++)
1324 accum
= (accum
<<8) | in
[i
];
1333 /* Read a binary value of the specified type, byte order, and byte size from
1334 * a file, converting it to a double. For integer types, the significant
1335 * bits are used to normalize the result. The sign of bits determines
1336 * whether they are padded toward the MSB (negative) or LSB (positive).
1337 * Floating-point types are not normalized.
1339 static int ReadBinAsDouble(FILE *fp
, const char *filename
, const ByteOrderT order
, const ElementTypeT type
, const uint bytes
, const int bits
, double *out
)
1354 if(!ReadBin8(fp
, filename
, order
, &v8
.ui
))
1361 if(!ReadBin4(fp
, filename
, order
, bytes
, &v4
.ui
))
1368 v4
.ui
>>= (8*bytes
) - (static_cast<uint
>(bits
));
1370 v4
.ui
&= (0xFFFFFFFF >> (32+bits
));
1372 if(v4
.ui
&static_cast<uint
>(1<<(std::abs(bits
)-1)))
1373 v4
.ui
|= (0xFFFFFFFF << std::abs(bits
));
1374 *out
= v4
.i
/ static_cast<double>(1<<(std::abs(bits
)-1));
1380 /* Read an ascii value of the specified type from a file, converting it to a
1381 * double. For integer types, the significant bits are used to normalize the
1382 * result. The sign of the bits should always be positive. This also skips
1383 * up to one separator character before the element itself.
1385 static int ReadAsciiAsDouble(TokenReaderT
*tr
, const char *filename
, const ElementTypeT type
, const uint bits
, double *out
)
1387 if(TrIsOperator(tr
, ","))
1388 TrReadOperator(tr
, ",");
1389 else if(TrIsOperator(tr
, ":"))
1390 TrReadOperator(tr
, ":");
1391 else if(TrIsOperator(tr
, ";"))
1392 TrReadOperator(tr
, ";");
1393 else if(TrIsOperator(tr
, "|"))
1394 TrReadOperator(tr
, "|");
1398 if(!TrReadFloat(tr
, -std::numeric_limits
<double>::infinity(),
1399 std::numeric_limits
<double>::infinity(), out
))
1401 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1408 if(!TrReadInt(tr
, -(1<<(bits
-1)), (1<<(bits
-1))-1, &v
))
1410 fprintf(stderr
, "Error: Bad read from file '%s'.\n", filename
);
1413 *out
= v
/ static_cast<double>((1<<(bits
-1))-1);
1418 // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
1419 // the source parameters and data set metrics.
1420 static int ReadWaveFormat(FILE *fp
, const ByteOrderT order
, const uint hrirRate
, SourceRefT
*src
)
1422 uint32_t fourCC
, chunkSize
;
1423 uint32_t format
, channels
, rate
, dummy
, block
, size
, bits
;
1428 fseek(fp
, static_cast<long>(chunkSize
), SEEK_CUR
);
1429 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1430 !ReadBin4(fp
, src
->mPath
, order
, 4, &chunkSize
))
1432 } while(fourCC
!= FOURCC_FMT
);
1433 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &format
) ||
1434 !ReadBin4(fp
, src
->mPath
, order
, 2, &channels
) ||
1435 !ReadBin4(fp
, src
->mPath
, order
, 4, &rate
) ||
1436 !ReadBin4(fp
, src
->mPath
, order
, 4, &dummy
) ||
1437 !ReadBin4(fp
, src
->mPath
, order
, 2, &block
))
1442 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &size
))
1450 if(format
== WAVE_FORMAT_EXTENSIBLE
)
1452 fseek(fp
, 2, SEEK_CUR
);
1453 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &bits
))
1457 fseek(fp
, 4, SEEK_CUR
);
1458 if(!ReadBin4(fp
, src
->mPath
, order
, 2, &format
))
1460 fseek(fp
, static_cast<long>(chunkSize
- 26), SEEK_CUR
);
1466 fseek(fp
, static_cast<long>(chunkSize
- 16), SEEK_CUR
);
1468 fseek(fp
, static_cast<long>(chunkSize
- 14), SEEK_CUR
);
1470 if(format
!= WAVE_FORMAT_PCM
&& format
!= WAVE_FORMAT_IEEE_FLOAT
)
1472 fprintf(stderr
, "Error: Unsupported WAVE format in file '%s'.\n", src
->mPath
);
1475 if(src
->mChannel
>= channels
)
1477 fprintf(stderr
, "Error: Missing source channel in WAVE file '%s'.\n", src
->mPath
);
1480 if(rate
!= hrirRate
)
1482 fprintf(stderr
, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src
->mPath
);
1485 if(format
== WAVE_FORMAT_PCM
)
1487 if(size
< 2 || size
> 4)
1489 fprintf(stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
->mPath
);
1492 if(bits
< 16 || bits
> (8*size
))
1494 fprintf(stderr
, "Error: Bad significant bits in WAVE file '%s'.\n", src
->mPath
);
1497 src
->mType
= ET_INT
;
1501 if(size
!= 4 && size
!= 8)
1503 fprintf(stderr
, "Error: Unsupported sample size in WAVE file '%s'.\n", src
->mPath
);
1509 src
->mBits
= static_cast<int>(bits
);
1510 src
->mSkip
= channels
;
1514 // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
1515 static int ReadWaveData(FILE *fp
, const SourceRefT
*src
, const ByteOrderT order
, const uint n
, double *hrir
)
1517 int pre
, post
, skip
;
1520 pre
= static_cast<int>(src
->mSize
* src
->mChannel
);
1521 post
= static_cast<int>(src
->mSize
* (src
->mSkip
- src
->mChannel
- 1));
1523 for(i
= 0;i
< n
;i
++)
1527 fseek(fp
, skip
, SEEK_CUR
);
1528 if(!ReadBinAsDouble(fp
, src
->mPath
, order
, src
->mType
, src
->mSize
, src
->mBits
, &hrir
[i
]))
1533 fseek(fp
, skip
, SEEK_CUR
);
1537 // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
1539 static int ReadWaveList(FILE *fp
, const SourceRefT
*src
, const ByteOrderT order
, const uint n
, double *hrir
)
1541 uint32_t fourCC
, chunkSize
, listSize
, count
;
1542 uint block
, skip
, offset
, i
;
1547 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1548 !ReadBin4(fp
, src
->mPath
, order
, 4, &chunkSize
))
1551 if(fourCC
== FOURCC_DATA
)
1553 block
= src
->mSize
* src
->mSkip
;
1554 count
= chunkSize
/ block
;
1555 if(count
< (src
->mOffset
+ n
))
1557 fprintf(stderr
, "Error: Bad read from file '%s'.\n", src
->mPath
);
1560 fseek(fp
, static_cast<long>(src
->mOffset
* block
), SEEK_CUR
);
1561 if(!ReadWaveData(fp
, src
, order
, n
, &hrir
[0]))
1565 else if(fourCC
== FOURCC_LIST
)
1567 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
))
1570 if(fourCC
== FOURCC_WAVL
)
1574 fseek(fp
, static_cast<long>(chunkSize
), SEEK_CUR
);
1576 listSize
= chunkSize
;
1577 block
= src
->mSize
* src
->mSkip
;
1578 skip
= src
->mOffset
;
1581 while(offset
< n
&& listSize
> 8)
1583 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1584 !ReadBin4(fp
, src
->mPath
, order
, 4, &chunkSize
))
1586 listSize
-= 8 + chunkSize
;
1587 if(fourCC
== FOURCC_DATA
)
1589 count
= chunkSize
/ block
;
1592 fseek(fp
, static_cast<long>(skip
* block
), SEEK_CUR
);
1593 chunkSize
-= skip
* block
;
1596 if(count
> (n
- offset
))
1598 if(!ReadWaveData(fp
, src
, order
, count
, &hrir
[offset
]))
1600 chunkSize
-= count
* block
;
1602 lastSample
= hrir
[offset
- 1];
1610 else if(fourCC
== FOURCC_SLNT
)
1612 if(!ReadBin4(fp
, src
->mPath
, order
, 4, &count
))
1619 if(count
> (n
- offset
))
1621 for(i
= 0; i
< count
; i
++)
1622 hrir
[offset
+ i
] = lastSample
;
1632 fseek(fp
, static_cast<long>(chunkSize
), SEEK_CUR
);
1636 fprintf(stderr
, "Error: Bad read from file '%s'.\n", src
->mPath
);
1642 // Load a source HRIR from a RIFF/RIFX WAVE file.
1643 static int LoadWaveSource(FILE *fp
, SourceRefT
*src
, const uint hrirRate
, const uint n
, double *hrir
)
1645 uint32_t fourCC
, dummy
;
1648 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
) ||
1649 !ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &dummy
))
1651 if(fourCC
== FOURCC_RIFF
)
1653 else if(fourCC
== FOURCC_RIFX
)
1657 fprintf(stderr
, "Error: No RIFF/RIFX chunk in file '%s'.\n", src
->mPath
);
1661 if(!ReadBin4(fp
, src
->mPath
, BO_LITTLE
, 4, &fourCC
))
1663 if(fourCC
!= FOURCC_WAVE
)
1665 fprintf(stderr
, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src
->mPath
);
1668 if(!ReadWaveFormat(fp
, order
, hrirRate
, src
))
1670 if(!ReadWaveList(fp
, src
, order
, n
, hrir
))
1675 // Load a source HRIR from a binary file.
1676 static int LoadBinarySource(FILE *fp
, const SourceRefT
*src
, const ByteOrderT order
, const uint n
, double *hrir
)
1680 fseek(fp
, static_cast<long>(src
->mOffset
), SEEK_SET
);
1681 for(i
= 0;i
< n
;i
++)
1683 if(!ReadBinAsDouble(fp
, src
->mPath
, order
, src
->mType
, src
->mSize
, src
->mBits
, &hrir
[i
]))
1686 fseek(fp
, static_cast<long>(src
->mSkip
), SEEK_CUR
);
1691 // Load a source HRIR from an ASCII text file containing a list of elements
1692 // separated by whitespace or common list operators (',', ';', ':', '|').
1693 static int LoadAsciiSource(FILE *fp
, const SourceRefT
*src
, const uint n
, double *hrir
)
1699 TrSetup(fp
, nullptr, &tr
);
1700 for(i
= 0;i
< src
->mOffset
;i
++)
1702 if(!ReadAsciiAsDouble(&tr
, src
->mPath
, src
->mType
, static_cast<uint
>(src
->mBits
), &dummy
))
1705 for(i
= 0;i
< n
;i
++)
1707 if(!ReadAsciiAsDouble(&tr
, src
->mPath
, src
->mType
, static_cast<uint
>(src
->mBits
), &hrir
[i
]))
1709 for(j
= 0;j
< src
->mSkip
;j
++)
1711 if(!ReadAsciiAsDouble(&tr
, src
->mPath
, src
->mType
, static_cast<uint
>(src
->mBits
), &dummy
))
1718 // Load a source HRIR from a supported file type.
1719 static int LoadSource(SourceRefT
*src
, const uint hrirRate
, const uint n
, double *hrir
)
1724 if(src
->mFormat
== SF_ASCII
)
1725 fp
= fopen(src
->mPath
, "r");
1727 fp
= fopen(src
->mPath
, "rb");
1730 fprintf(stderr
, "Error: Could not open source file '%s'.\n", src
->mPath
);
1733 if(src
->mFormat
== SF_WAVE
)
1734 result
= LoadWaveSource(fp
, src
, hrirRate
, n
, hrir
);
1735 else if(src
->mFormat
== SF_BIN_LE
)
1736 result
= LoadBinarySource(fp
, src
, BO_LITTLE
, n
, hrir
);
1737 else if(src
->mFormat
== SF_BIN_BE
)
1738 result
= LoadBinarySource(fp
, src
, BO_BIG
, n
, hrir
);
1740 result
= LoadAsciiSource(fp
, src
, n
, hrir
);
1746 /***************************
1747 *** File storage output ***
1748 ***************************/
1750 // Write an ASCII string to a file.
1751 static int WriteAscii(const char *out
, FILE *fp
, const char *filename
)
1756 if(fwrite(out
, 1, len
, fp
) != len
)
1759 fprintf(stderr
, "Error: Bad write to file '%s'.\n", filename
);
1765 // Write a binary value of the given byte order and byte size to a file,
1766 // loading it from a 32-bit unsigned integer.
1767 static int WriteBin4(const ByteOrderT order
, const uint bytes
, const uint32_t in
, FILE *fp
, const char *filename
)
1775 for(i
= 0;i
< bytes
;i
++)
1776 out
[i
] = (in
>>(i
*8)) & 0x000000FF;
1779 for(i
= 0;i
< bytes
;i
++)
1780 out
[bytes
- i
- 1] = (in
>>(i
*8)) & 0x000000FF;
1785 if(fwrite(out
, 1, bytes
, fp
) != bytes
)
1787 fprintf(stderr
, "Error: Bad write to file '%s'.\n", filename
);
1793 // Store the OpenAL Soft HRTF data set.
1794 static int StoreMhr(const HrirDataT
*hData
, const char *filename
)
1796 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
1797 uint n
= hData
->mIrPoints
;
1800 uint dither_seed
= 22222;
1802 if((fp
=fopen(filename
, "wb")) == nullptr)
1804 fprintf(stderr
, "Error: Could not open MHR file '%s'.\n", filename
);
1807 if(!WriteAscii(MHR_FORMAT
, fp
, filename
))
1809 if(!WriteBin4(BO_LITTLE
, 4, static_cast<uint32_t>(hData
->mIrRate
), fp
, filename
))
1811 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(hData
->mSampleType
), fp
, filename
))
1813 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(hData
->mChannelType
), fp
, filename
))
1815 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(hData
->mIrPoints
), fp
, filename
))
1817 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(hData
->mFdCount
), fp
, filename
))
1819 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1821 if(!WriteBin4(BO_LITTLE
, 2, static_cast<uint32_t>(1000.0 * hData
->mFds
[fi
].mDistance
), fp
, filename
))
1823 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(hData
->mFds
[fi
].mEvCount
), fp
, filename
))
1825 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
1827 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(hData
->mFds
[fi
].mEvs
[ei
].mAzCount
), fp
, filename
))
1832 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1834 const double scale
= (hData
->mSampleType
== ST_S16
) ? 32767.0 :
1835 ((hData
->mSampleType
== ST_S24
) ? 8388607.0 : 0.0);
1836 const int bps
= (hData
->mSampleType
== ST_S16
) ? 2 :
1837 ((hData
->mSampleType
== ST_S24
) ? 3 : 0);
1839 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
1841 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
1843 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
1844 double out
[2 * MAX_TRUNCSIZE
];
1846 TpdfDither(out
, azd
->mIrs
[0], scale
, n
, channels
, &dither_seed
);
1847 if(hData
->mChannelType
== CT_STEREO
)
1848 TpdfDither(out
+1, azd
->mIrs
[1], scale
, n
, channels
, &dither_seed
);
1849 for(i
= 0;i
< (channels
* n
);i
++)
1851 int v
= static_cast<int>(Clamp(out
[i
], -scale
-1.0, scale
));
1852 if(!WriteBin4(BO_LITTLE
, bps
, static_cast<uint32_t>(v
), fp
, filename
))
1858 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1860 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
1862 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
1864 const HrirAzT
&azd
= hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
1865 int v
= static_cast<int>(std::min(std::round(hData
->mIrRate
* azd
.mDelays
[0]), MAX_HRTD
));
1867 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(v
), fp
, filename
))
1869 if(hData
->mChannelType
== CT_STEREO
)
1871 v
= static_cast<int>(std::min(std::round(hData
->mIrRate
* azd
.mDelays
[1]), MAX_HRTD
));
1873 if(!WriteBin4(BO_LITTLE
, 1, static_cast<uint32_t>(v
), fp
, filename
))
1884 /***********************
1885 *** HRTF processing ***
1886 ***********************/
1888 // Calculate the onset time of an HRIR and average it with any existing
1889 // timing for its field, elevation, azimuth, and ear.
1890 static double AverageHrirOnset(const uint rate
, const uint n
, const double *hrir
, const double f
, const double onset
)
1895 for(i
= 0;i
< n
;i
++)
1896 mag
= std::max(std::abs(hrir
[i
]), mag
);
1898 for(i
= 0;i
< n
;i
++)
1900 if(std::abs(hrir
[i
]) >= mag
)
1903 return Lerp(onset
, static_cast<double>(i
) / rate
, f
);
1906 // Calculate the magnitude response of an HRIR and average it with any
1907 // existing responses for its field, elevation, azimuth, and ear.
1908 static void AverageHrirMagnitude(const uint points
, const uint n
, const double *hrir
, const double f
, double *mag
)
1910 uint m
= 1 + (n
/ 2), i
;
1911 std::vector
<complex_d
> h(n
);
1912 std::vector
<double> r(n
);
1914 for(i
= 0;i
< points
;i
++)
1915 h
[i
] = complex_d
{hrir
[i
], 0.0};
1917 h
[i
] = complex_d
{0.0, 0.0};
1918 FftForward(n
, h
.data());
1919 MagnitudeResponse(n
, h
.data(), r
.data());
1920 for(i
= 0;i
< m
;i
++)
1921 mag
[i
] = Lerp(mag
[i
], r
[i
], f
);
1924 /* Calculate the contribution of each HRIR to the diffuse-field average based
1925 * on the area of its surface patch. All patches are centered at the HRIR
1926 * coordinates on the unit sphere and are measured by solid angle.
1928 static void CalculateDfWeights(const HrirDataT
*hData
, double *weights
)
1930 double sum
, evs
, ev
, upperEv
, lowerEv
, solidAngle
;
1934 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1936 evs
= M_PI
/ 2.0 / (hData
->mFds
[fi
].mEvCount
- 1);
1937 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
1939 // For each elevation, calculate the upper and lower limits of
1941 ev
= hData
->mFds
[fi
].mEvs
[ei
].mElevation
;
1942 lowerEv
= std::max(-M_PI
/ 2.0, ev
- evs
);
1943 upperEv
= std::min(M_PI
/ 2.0, ev
+ evs
);
1944 // Calculate the area of the patch band.
1945 solidAngle
= 2.0 * M_PI
* (std::sin(upperEv
) - std::sin(lowerEv
));
1946 // Each weight is the area of one patch.
1947 weights
[(fi
* MAX_EV_COUNT
) + ei
] = solidAngle
/ hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;
1948 // Sum the total surface area covered by the HRIRs of all fields.
1952 /* TODO: It may be interesting to experiment with how a volume-based
1953 weighting performs compared to the existing distance-indepenent
1956 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1958 // Normalize the weights given the total surface coverage for all
1960 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
1961 weights
[(fi
* MAX_EV_COUNT
) + ei
] /= sum
;
1965 /* Calculate the diffuse-field average from the given magnitude responses of
1966 * the HRIR set. Weighting can be applied to compensate for the varying
1967 * surface area covered by each HRIR. The final average can then be limited
1968 * by the specified magnitude range (in positive dB; 0.0 to skip).
1970 static void CalculateDiffuseFieldAverage(const HrirDataT
*hData
, const uint channels
, const uint m
, const int weighted
, const double limit
, double *dfa
)
1972 std::vector
<double> weights(hData
->mFdCount
* MAX_EV_COUNT
);
1973 uint count
, ti
, fi
, ei
, i
, ai
;
1977 // Use coverage weighting to calculate the average.
1978 CalculateDfWeights(hData
, weights
.data());
1984 // If coverage weighting is not used, the weights still need to be
1985 // averaged by the number of existing HRIRs.
1986 count
= hData
->mIrCount
;
1987 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1989 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvStart
;ei
++)
1990 count
-= hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;
1992 weight
= 1.0 / count
;
1994 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
1996 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
1997 weights
[(fi
* MAX_EV_COUNT
) + ei
] = weight
;
2000 for(ti
= 0;ti
< channels
;ti
++)
2002 for(i
= 0;i
< m
;i
++)
2003 dfa
[(ti
* m
) + i
] = 0.0;
2004 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2006 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2008 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2010 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2011 // Get the weight for this HRIR's contribution.
2012 double weight
= weights
[(fi
* MAX_EV_COUNT
) + ei
];
2014 // Add this HRIR's weighted power average to the total.
2015 for(i
= 0;i
< m
;i
++)
2016 dfa
[(ti
* m
) + i
] += weight
* azd
->mIrs
[ti
][i
] * azd
->mIrs
[ti
][i
];
2020 // Finish the average calculation and keep it from being too small.
2021 for(i
= 0;i
< m
;i
++)
2022 dfa
[(ti
* m
) + i
] = std::max(sqrt(dfa
[(ti
* m
) + i
]), EPSILON
);
2023 // Apply a limit to the magnitude range of the diffuse-field average
2026 LimitMagnitudeResponse(hData
->mFftSize
, m
, limit
, &dfa
[ti
* m
], &dfa
[ti
* m
]);
2030 // Perform diffuse-field equalization on the magnitude responses of the HRIR
2031 // set using the given average response.
2032 static void DiffuseFieldEqualize(const uint channels
, const uint m
, const double *dfa
, const HrirDataT
*hData
)
2034 uint ti
, fi
, ei
, ai
, i
;
2036 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2038 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2040 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2042 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2044 for(ti
= 0;ti
< channels
;ti
++)
2046 for(i
= 0;i
< m
;i
++)
2047 azd
->mIrs
[ti
][i
] /= dfa
[(ti
* m
) + i
];
2054 // Perform minimum-phase reconstruction using the magnitude responses of the
2056 static void ReconstructHrirs(const HrirDataT
*hData
)
2058 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2059 uint n
= hData
->mFftSize
;
2060 uint ti
, fi
, ei
, ai
, i
;
2061 std::vector
<complex_d
> h(n
);
2062 uint total
, count
, pcdone
, lastpc
;
2064 total
= hData
->mIrCount
;
2065 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2067 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvStart
;ei
++)
2068 total
-= hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;
2071 count
= pcdone
= lastpc
= 0;
2072 printf("%3d%% done.", pcdone
);
2074 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2076 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2078 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2080 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2082 for(ti
= 0;ti
< channels
;ti
++)
2084 MinimumPhase(n
, azd
->mIrs
[ti
], h
.data());
2085 FftInverse(n
, h
.data());
2086 for(i
= 0;i
< hData
->mIrPoints
;i
++)
2087 azd
->mIrs
[ti
][i
] = h
[i
].real();
2088 pcdone
= ++count
* 100 / total
;
2089 if(pcdone
!= lastpc
)
2092 printf("\r%3d%% done.", pcdone
);
2102 // Resamples the HRIRs for use at the given sampling rate.
2103 static void ResampleHrirs(const uint rate
, HrirDataT
*hData
)
2105 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2106 uint n
= hData
->mIrPoints
;
2107 uint ti
, fi
, ei
, ai
;
2110 ResamplerSetup(&rs
, hData
->mIrRate
, rate
);
2111 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2113 for(ei
= hData
->mFds
[fi
].mEvStart
;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2115 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2117 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2118 for(ti
= 0;ti
< channels
;ti
++)
2119 ResamplerRun(&rs
, n
, azd
->mIrs
[ti
], n
, azd
->mIrs
[ti
]);
2123 hData
->mIrRate
= rate
;
2126 /* Given field and elevation indices and an azimuth, calculate the indices of
2127 * the two HRIRs that bound the coordinate along with a factor for
2128 * calculating the continuous HRIR using interpolation.
2130 static void CalcAzIndices(const HrirDataT
*hData
, const uint fi
, const uint ei
, const double az
, uint
*a0
, uint
*a1
, double *af
)
2132 double f
= (2.0*M_PI
+ az
) * hData
->mFds
[fi
].mEvs
[ei
].mAzCount
/ (2.0*M_PI
);
2133 uint i
= static_cast<uint
>(f
) % hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;
2137 *a1
= (i
+ 1) % hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;
2141 // Synthesize any missing onset timings at the bottom elevations of each
2142 // field. This just blends between slightly exaggerated known onsets (not
2143 // an accurate model).
2144 static void SynthesizeOnsets(HrirDataT
*hData
)
2146 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2147 uint ti
, fi
, oi
, ai
, ei
, a0
, a1
;
2150 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2152 if(hData
->mFds
[fi
].mEvStart
<= 0)
2154 oi
= hData
->mFds
[fi
].mEvStart
;
2156 for(ti
= 0;ti
< channels
;ti
++)
2159 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[oi
].mAzCount
;ai
++)
2160 t
+= hData
->mFds
[fi
].mEvs
[oi
].mAzs
[ai
].mDelays
[ti
];
2161 hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mDelays
[ti
] = 1.32e-4 + (t
/ hData
->mFds
[fi
].mEvs
[oi
].mAzCount
);
2162 for(ei
= 1;ei
< hData
->mFds
[fi
].mEvStart
;ei
++)
2164 of
= static_cast<double>(ei
) / hData
->mFds
[fi
].mEvStart
;
2165 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2167 CalcAzIndices(hData
, fi
, oi
, hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mAzimuth
, &a0
, &a1
, &af
);
2168 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mDelays
[ti
] = Lerp(
2169 hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mDelays
[ti
],
2170 Lerp(hData
->mFds
[fi
].mEvs
[oi
].mAzs
[a0
].mDelays
[ti
],
2171 hData
->mFds
[fi
].mEvs
[oi
].mAzs
[a1
].mDelays
[ti
], af
),
2180 /* Attempt to synthesize any missing HRIRs at the bottom elevations of each
2181 * field. Right now this just blends the lowest elevation HRIRs together and
2182 * applies some attenuation and high frequency damping. It is a simple, if
2185 static void SynthesizeHrirs(HrirDataT
*hData
)
2187 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2188 uint n
= hData
->mIrPoints
;
2189 uint ti
, fi
, ai
, ei
, i
;
2190 double lp
[4], s0
, s1
;
2195 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2197 const uint oi
= hData
->mFds
[fi
].mEvStart
;
2198 if(oi
<= 0) continue;
2200 for(ti
= 0;ti
< channels
;ti
++)
2202 for(i
= 0;i
< n
;i
++)
2203 hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mIrs
[ti
][i
] = 0.0;
2204 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[oi
].mAzCount
;ai
++)
2206 for(i
= 0;i
< n
;i
++)
2207 hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mIrs
[ti
][i
] += hData
->mFds
[fi
].mEvs
[oi
].mAzs
[ai
].mIrs
[ti
][i
] /
2208 hData
->mFds
[fi
].mEvs
[oi
].mAzCount
;
2210 for(ei
= 1;ei
< hData
->mFds
[fi
].mEvStart
;ei
++)
2212 of
= static_cast<double>(ei
) / hData
->mFds
[fi
].mEvStart
;
2213 b
= (1.0 - of
) * (3.5e-6 * hData
->mIrRate
);
2214 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2216 CalcAzIndices(hData
, fi
, oi
, hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mAzimuth
, &a0
, &a1
, &af
);
2221 for(i
= 0;i
< n
;i
++)
2223 s0
= hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mIrs
[ti
][i
];
2224 s1
= Lerp(hData
->mFds
[fi
].mEvs
[oi
].mAzs
[a0
].mIrs
[ti
][i
],
2225 hData
->mFds
[fi
].mEvs
[oi
].mAzs
[a1
].mIrs
[ti
][i
], af
);
2226 s0
= Lerp(s0
, s1
, of
);
2227 lp
[0] = Lerp(s0
, lp
[0], b
);
2228 lp
[1] = Lerp(lp
[0], lp
[1], b
);
2229 lp
[2] = Lerp(lp
[1], lp
[2], b
);
2230 lp
[3] = Lerp(lp
[2], lp
[3], b
);
2231 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mIrs
[ti
][i
] = lp
[3];
2235 b
= 3.5e-6 * hData
->mIrRate
;
2240 for(i
= 0;i
< n
;i
++)
2242 s0
= hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mIrs
[ti
][i
];
2243 lp
[0] = Lerp(s0
, lp
[0], b
);
2244 lp
[1] = Lerp(lp
[0], lp
[1], b
);
2245 lp
[2] = Lerp(lp
[1], lp
[2], b
);
2246 lp
[3] = Lerp(lp
[2], lp
[3], b
);
2247 hData
->mFds
[fi
].mEvs
[0].mAzs
[0].mIrs
[ti
][i
] = lp
[3];
2250 hData
->mFds
[fi
].mEvStart
= 0;
2254 // The following routines assume a full set of HRIRs for all elevations.
2256 // Normalize the HRIR set and slightly attenuate the result.
2257 static void NormalizeHrirs(const HrirDataT
*hData
)
2259 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2260 uint n
= hData
->mIrPoints
;
2261 uint ti
, fi
, ei
, ai
, i
;
2262 double maxLevel
= 0.0;
2264 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2266 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2268 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2270 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2271 for(ti
= 0;ti
< channels
;ti
++)
2273 for(i
= 0;i
< n
;i
++)
2274 maxLevel
= std::max(std::abs(azd
->mIrs
[ti
][i
]), maxLevel
);
2279 maxLevel
= 1.01 * maxLevel
;
2280 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2282 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2284 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2286 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2288 for(ti
= 0;ti
< channels
;ti
++)
2290 for(i
= 0;i
< n
;i
++)
2291 azd
->mIrs
[ti
][i
] /= maxLevel
;
2298 // Calculate the left-ear time delay using a spherical head model.
2299 static double CalcLTD(const double ev
, const double az
, const double rad
, const double dist
)
2301 double azp
, dlp
, l
, al
;
2303 azp
= std::asin(std::cos(ev
) * std::sin(az
));
2304 dlp
= std::sqrt((dist
*dist
) + (rad
*rad
) + (2.0*dist
*rad
*sin(azp
)));
2305 l
= std::sqrt((dist
*dist
) - (rad
*rad
));
2306 al
= (0.5 * M_PI
) + azp
;
2308 dlp
= l
+ (rad
* (al
- std::acos(rad
/ dist
)));
2312 // Calculate the effective head-related time delays for each minimum-phase
2314 static void CalculateHrtds(const HeadModelT model
, const double radius
, HrirDataT
*hData
)
2316 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2317 double minHrtd
{std::numeric_limits
<double>::infinity()};
2318 double maxHrtd
{-minHrtd
};
2319 uint ti
, fi
, ei
, ai
;
2322 if(model
== HM_DATASET
)
2324 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2326 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2328 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2330 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2331 for(ti
= 0;ti
< channels
;ti
++)
2333 t
= azd
->mDelays
[ti
] * radius
/ hData
->mRadius
;
2334 azd
->mDelays
[ti
] = t
;
2335 maxHrtd
= std::max(t
, maxHrtd
);
2336 minHrtd
= std::min(t
, minHrtd
);
2344 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2346 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2348 HrirEvT
*evd
= &hData
->mFds
[fi
].mEvs
[ei
];
2349 for(ai
= 0;ai
< evd
->mAzCount
;ai
++)
2351 HrirAzT
*azd
= &evd
->mAzs
[ai
];
2352 for(ti
= 0;ti
< channels
;ti
++)
2354 t
= CalcLTD(evd
->mElevation
, azd
->mAzimuth
, radius
, hData
->mFds
[fi
].mDistance
);
2355 azd
->mDelays
[ti
] = t
;
2356 maxHrtd
= std::max(t
, maxHrtd
);
2357 minHrtd
= std::min(t
, minHrtd
);
2363 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2365 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2367 for(ti
= 0;ti
< channels
;ti
++)
2369 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2370 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mDelays
[ti
] -= minHrtd
;
2376 // Allocate and configure dynamic HRIR structures.
2377 static int PrepareHrirData(const uint fdCount
, const double distances
[MAX_FD_COUNT
], const uint evCounts
[MAX_FD_COUNT
], const uint azCounts
[MAX_FD_COUNT
* MAX_EV_COUNT
], HrirDataT
*hData
)
2379 uint evTotal
= 0, azTotal
= 0, fi
, ei
, ai
;
2381 for(fi
= 0;fi
< fdCount
;fi
++)
2383 evTotal
+= evCounts
[fi
];
2384 for(ei
= 0;ei
< evCounts
[fi
];ei
++)
2385 azTotal
+= azCounts
[(fi
* MAX_EV_COUNT
) + ei
];
2387 if(!fdCount
|| !evTotal
|| !azTotal
)
2390 hData
->mEvsBase
.resize(evTotal
);
2391 hData
->mAzsBase
.resize(azTotal
);
2392 hData
->mFds
.resize(fdCount
);
2393 hData
->mIrCount
= azTotal
;
2394 hData
->mFdCount
= fdCount
;
2397 for(fi
= 0;fi
< fdCount
;fi
++)
2399 hData
->mFds
[fi
].mDistance
= distances
[fi
];
2400 hData
->mFds
[fi
].mEvCount
= evCounts
[fi
];
2401 hData
->mFds
[fi
].mEvStart
= 0;
2402 hData
->mFds
[fi
].mEvs
= &hData
->mEvsBase
[evTotal
];
2403 evTotal
+= evCounts
[fi
];
2404 for(ei
= 0;ei
< evCounts
[fi
];ei
++)
2406 uint azCount
= azCounts
[(fi
* MAX_EV_COUNT
) + ei
];
2408 hData
->mFds
[fi
].mIrCount
+= azCount
;
2409 hData
->mFds
[fi
].mEvs
[ei
].mElevation
= -M_PI
/ 2.0 + M_PI
* ei
/ (evCounts
[fi
] - 1);
2410 hData
->mFds
[fi
].mEvs
[ei
].mIrCount
+= azCount
;
2411 hData
->mFds
[fi
].mEvs
[ei
].mAzCount
= azCount
;
2412 hData
->mFds
[fi
].mEvs
[ei
].mAzs
= &hData
->mAzsBase
[azTotal
];
2413 for(ai
= 0;ai
< azCount
;ai
++)
2415 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mAzimuth
= 2.0 * M_PI
* ai
/ azCount
;
2416 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mIndex
= azTotal
+ ai
;
2417 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mDelays
[0] = 0.0;
2418 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mDelays
[1] = 0.0;
2419 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mIrs
[0] = nullptr;
2420 hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
].mIrs
[1] = nullptr;
2428 // Match the channel type from a given identifier.
2429 static ChannelTypeT
MatchChannelType(const char *ident
)
2431 if(strcasecmp(ident
, "mono") == 0)
2433 if(strcasecmp(ident
, "stereo") == 0)
2438 // Process the data set definition to read and validate the data set metrics.
2439 static int ProcessMetrics(TokenReaderT
*tr
, const uint fftSize
, const uint truncSize
, HrirDataT
*hData
)
2441 int hasRate
= 0, hasType
= 0, hasPoints
= 0, hasRadius
= 0;
2442 int hasDistance
= 0, hasAzimuths
= 0;
2443 char ident
[MAX_IDENT_LEN
+1];
2448 double distances
[MAX_FD_COUNT
];
2450 uint evCounts
[MAX_FD_COUNT
];
2451 std::vector
<uint
> azCounts(MAX_FD_COUNT
* MAX_EV_COUNT
);
2453 TrIndication(tr
, &line
, &col
);
2454 while(TrIsIdent(tr
))
2456 TrIndication(tr
, &line
, &col
);
2457 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2459 if(strcasecmp(ident
, "rate") == 0)
2463 TrErrorAt(tr
, line
, col
, "Redefinition of 'rate'.\n");
2466 if(!TrReadOperator(tr
, "="))
2468 if(!TrReadInt(tr
, MIN_RATE
, MAX_RATE
, &intVal
))
2470 hData
->mIrRate
= static_cast<uint
>(intVal
);
2473 else if(strcasecmp(ident
, "type") == 0)
2475 char type
[MAX_IDENT_LEN
+1];
2479 TrErrorAt(tr
, line
, col
, "Redefinition of 'type'.\n");
2482 if(!TrReadOperator(tr
, "="))
2485 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, type
))
2487 hData
->mChannelType
= MatchChannelType(type
);
2488 if(hData
->mChannelType
== CT_NONE
)
2490 TrErrorAt(tr
, line
, col
, "Expected a channel type.\n");
2495 else if(strcasecmp(ident
, "points") == 0)
2499 TrErrorAt(tr
, line
, col
, "Redefinition of 'points'.\n");
2502 if(!TrReadOperator(tr
, "="))
2504 TrIndication(tr
, &line
, &col
);
2505 if(!TrReadInt(tr
, MIN_POINTS
, MAX_POINTS
, &intVal
))
2507 points
= static_cast<uint
>(intVal
);
2508 if(fftSize
> 0 && points
> fftSize
)
2510 TrErrorAt(tr
, line
, col
, "Value exceeds the overridden FFT size.\n");
2513 if(points
< truncSize
)
2515 TrErrorAt(tr
, line
, col
, "Value is below the truncation size.\n");
2518 hData
->mIrPoints
= points
;
2521 hData
->mFftSize
= DEFAULT_FFTSIZE
;
2522 hData
->mIrSize
= 1 + (DEFAULT_FFTSIZE
/ 2);
2526 hData
->mFftSize
= fftSize
;
2527 hData
->mIrSize
= 1 + (fftSize
/ 2);
2528 if(points
> hData
->mIrSize
)
2529 hData
->mIrSize
= points
;
2533 else if(strcasecmp(ident
, "radius") == 0)
2537 TrErrorAt(tr
, line
, col
, "Redefinition of 'radius'.\n");
2540 if(!TrReadOperator(tr
, "="))
2542 if(!TrReadFloat(tr
, MIN_RADIUS
, MAX_RADIUS
, &fpVal
))
2544 hData
->mRadius
= fpVal
;
2547 else if(strcasecmp(ident
, "distance") == 0)
2553 TrErrorAt(tr
, line
, col
, "Redefinition of 'distance'.\n");
2556 if(!TrReadOperator(tr
, "="))
2561 if(!TrReadFloat(tr
, MIN_DISTANCE
, MAX_DISTANCE
, &fpVal
))
2563 if(count
> 0 && fpVal
<= distances
[count
- 1])
2565 TrError(tr
, "Distances are not ascending.\n");
2568 distances
[count
++] = fpVal
;
2569 if(!TrIsOperator(tr
, ","))
2571 if(count
>= MAX_FD_COUNT
)
2573 TrError(tr
, "Exceeded the maximum of %d fields.\n", MAX_FD_COUNT
);
2576 TrReadOperator(tr
, ",");
2578 if(fdCount
!= 0 && count
!= fdCount
)
2580 TrError(tr
, "Did not match the specified number of %d fields.\n", fdCount
);
2586 else if(strcasecmp(ident
, "azimuths") == 0)
2592 TrErrorAt(tr
, line
, col
, "Redefinition of 'azimuths'.\n");
2595 if(!TrReadOperator(tr
, "="))
2601 if(!TrReadInt(tr
, MIN_AZ_COUNT
, MAX_AZ_COUNT
, &intVal
))
2603 azCounts
[(count
* MAX_EV_COUNT
) + evCounts
[count
]++] = static_cast<uint
>(intVal
);
2604 if(TrIsOperator(tr
, ","))
2606 if(evCounts
[count
] >= MAX_EV_COUNT
)
2608 TrError(tr
, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT
);
2611 TrReadOperator(tr
, ",");
2615 if(evCounts
[count
] < MIN_EV_COUNT
)
2617 TrErrorAt(tr
, line
, col
, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT
);
2620 if(azCounts
[count
* MAX_EV_COUNT
] != 1 || azCounts
[(count
* MAX_EV_COUNT
) + evCounts
[count
] - 1] != 1)
2622 TrError(tr
, "Poles are not singular for field %d.\n", count
- 1);
2626 if(!TrIsOperator(tr
, ";"))
2629 if(count
>= MAX_FD_COUNT
)
2631 TrError(tr
, "Exceeded the maximum number of %d fields.\n", MAX_FD_COUNT
);
2634 evCounts
[count
] = 0;
2635 TrReadOperator(tr
, ";");
2638 if(fdCount
!= 0 && count
!= fdCount
)
2640 TrError(tr
, "Did not match the specified number of %d fields.\n", fdCount
);
2648 TrErrorAt(tr
, line
, col
, "Expected a metric name.\n");
2651 TrSkipWhitespace(tr
);
2653 if(!(hasRate
&& hasPoints
&& hasRadius
&& hasDistance
&& hasAzimuths
))
2655 TrErrorAt(tr
, line
, col
, "Expected a metric name.\n");
2658 if(distances
[0] < hData
->mRadius
)
2660 TrError(tr
, "Distance cannot start below head radius.\n");
2663 if(hData
->mChannelType
== CT_NONE
)
2664 hData
->mChannelType
= CT_MONO
;
2665 if(!PrepareHrirData(fdCount
, distances
, evCounts
, azCounts
.data(), hData
))
2667 fprintf(stderr
, "Error: Out of memory.\n");
2673 // Parse an index triplet from the data set definition.
2674 static int ReadIndexTriplet(TokenReaderT
*tr
, const HrirDataT
*hData
, uint
*fi
, uint
*ei
, uint
*ai
)
2678 if(hData
->mFdCount
> 1)
2680 if(!TrReadInt(tr
, 0, static_cast<int>(hData
->mFdCount
) - 1, &intVal
))
2682 *fi
= static_cast<uint
>(intVal
);
2683 if(!TrReadOperator(tr
, ","))
2690 if(!TrReadInt(tr
, 0, static_cast<int>(hData
->mFds
[*fi
].mEvCount
) - 1, &intVal
))
2692 *ei
= static_cast<uint
>(intVal
);
2693 if(!TrReadOperator(tr
, ","))
2695 if(!TrReadInt(tr
, 0, static_cast<int>(hData
->mFds
[*fi
].mEvs
[*ei
].mAzCount
) - 1, &intVal
))
2697 *ai
= static_cast<uint
>(intVal
);
2701 // Match the source format from a given identifier.
2702 static SourceFormatT
MatchSourceFormat(const char *ident
)
2704 if(strcasecmp(ident
, "wave") == 0)
2706 if(strcasecmp(ident
, "bin_le") == 0)
2708 if(strcasecmp(ident
, "bin_be") == 0)
2710 if(strcasecmp(ident
, "ascii") == 0)
2715 // Match the source element type from a given identifier.
2716 static ElementTypeT
MatchElementType(const char *ident
)
2718 if(strcasecmp(ident
, "int") == 0)
2720 if(strcasecmp(ident
, "fp") == 0)
2725 // Parse and validate a source reference from the data set definition.
2726 static int ReadSourceRef(TokenReaderT
*tr
, SourceRefT
*src
)
2728 char ident
[MAX_IDENT_LEN
+1];
2732 TrIndication(tr
, &line
, &col
);
2733 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2735 src
->mFormat
= MatchSourceFormat(ident
);
2736 if(src
->mFormat
== SF_NONE
)
2738 TrErrorAt(tr
, line
, col
, "Expected a source format.\n");
2741 if(!TrReadOperator(tr
, "("))
2743 if(src
->mFormat
== SF_WAVE
)
2745 if(!TrReadInt(tr
, 0, MAX_WAVE_CHANNELS
, &intVal
))
2747 src
->mType
= ET_NONE
;
2750 src
->mChannel
= static_cast<uint
>(intVal
);
2755 TrIndication(tr
, &line
, &col
);
2756 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2758 src
->mType
= MatchElementType(ident
);
2759 if(src
->mType
== ET_NONE
)
2761 TrErrorAt(tr
, line
, col
, "Expected a source element type.\n");
2764 if(src
->mFormat
== SF_BIN_LE
|| src
->mFormat
== SF_BIN_BE
)
2766 if(!TrReadOperator(tr
, ","))
2768 if(src
->mType
== ET_INT
)
2770 if(!TrReadInt(tr
, MIN_BIN_SIZE
, MAX_BIN_SIZE
, &intVal
))
2772 src
->mSize
= static_cast<uint
>(intVal
);
2773 if(!TrIsOperator(tr
, ","))
2774 src
->mBits
= static_cast<int>(8*src
->mSize
);
2777 TrReadOperator(tr
, ",");
2778 TrIndication(tr
, &line
, &col
);
2779 if(!TrReadInt(tr
, -2147483647-1, 2147483647, &intVal
))
2781 if(std::abs(intVal
) < MIN_BIN_BITS
|| static_cast<uint
>(std::abs(intVal
)) > (8*src
->mSize
))
2783 TrErrorAt(tr
, line
, col
, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS
, 8*src
->mSize
);
2786 src
->mBits
= intVal
;
2791 TrIndication(tr
, &line
, &col
);
2792 if(!TrReadInt(tr
, -2147483647-1, 2147483647, &intVal
))
2794 if(intVal
!= 4 && intVal
!= 8)
2796 TrErrorAt(tr
, line
, col
, "Expected a value of 4 or 8.\n");
2799 src
->mSize
= static_cast<uint
>(intVal
);
2803 else if(src
->mFormat
== SF_ASCII
&& src
->mType
== ET_INT
)
2805 if(!TrReadOperator(tr
, ","))
2807 if(!TrReadInt(tr
, MIN_ASCII_BITS
, MAX_ASCII_BITS
, &intVal
))
2810 src
->mBits
= intVal
;
2818 if(!TrIsOperator(tr
, ";"))
2822 TrReadOperator(tr
, ";");
2823 if(!TrReadInt(tr
, 0, 0x7FFFFFFF, &intVal
))
2825 src
->mSkip
= static_cast<uint
>(intVal
);
2828 if(!TrReadOperator(tr
, ")"))
2830 if(TrIsOperator(tr
, "@"))
2832 TrReadOperator(tr
, "@");
2833 if(!TrReadInt(tr
, 0, 0x7FFFFFFF, &intVal
))
2835 src
->mOffset
= static_cast<uint
>(intVal
);
2839 if(!TrReadOperator(tr
, ":"))
2841 if(!TrReadString(tr
, MAX_PATH_LEN
, src
->mPath
))
2846 // Match the target ear (index) from a given identifier.
2847 static int MatchTargetEar(const char *ident
)
2849 if(strcasecmp(ident
, "left") == 0)
2851 if(strcasecmp(ident
, "right") == 0)
2856 // Process the list of sources in the data set definition.
2857 static int ProcessSources(const HeadModelT model
, TokenReaderT
*tr
, HrirDataT
*hData
)
2859 uint channels
= (hData
->mChannelType
== CT_STEREO
) ? 2 : 1;
2860 hData
->mHrirsBase
.resize(channels
* hData
->mIrCount
* hData
->mIrSize
);
2861 double *hrirs
= hData
->mHrirsBase
.data();
2862 std::vector
<double> hrir(hData
->mIrPoints
);
2863 uint line
, col
, fi
, ei
, ai
, ti
;
2866 printf("Loading sources...");
2869 while(TrIsOperator(tr
, "["))
2871 double factor
[2]{ 1.0, 1.0 };
2873 TrIndication(tr
, &line
, &col
);
2874 TrReadOperator(tr
, "[");
2875 if(!ReadIndexTriplet(tr
, hData
, &fi
, &ei
, &ai
))
2877 if(!TrReadOperator(tr
, "]"))
2879 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2881 if(azd
->mIrs
[0] != nullptr)
2883 TrErrorAt(tr
, line
, col
, "Redefinition of source.\n");
2886 if(!TrReadOperator(tr
, "="))
2894 if(!ReadSourceRef(tr
, &src
))
2897 // TODO: Would be nice to display 'x of y files', but that would
2898 // require preparing the source refs first to get a total count
2899 // before loading them.
2901 printf("\rLoading sources... %d file%s", count
, (count
==1)?"":"s");
2904 if(!LoadSource(&src
, hData
->mIrRate
, hData
->mIrPoints
, hrir
.data()))
2907 if(hData
->mChannelType
== CT_STEREO
)
2909 char ident
[MAX_IDENT_LEN
+1];
2911 if(!TrReadIdent(tr
, MAX_IDENT_LEN
, ident
))
2913 ti
= MatchTargetEar(ident
);
2914 if(static_cast<int>(ti
) < 0)
2916 TrErrorAt(tr
, line
, col
, "Expected a target ear.\n");
2920 azd
->mIrs
[ti
] = &hrirs
[hData
->mIrSize
* (ti
* hData
->mIrCount
+ azd
->mIndex
)];
2921 if(model
== HM_DATASET
)
2922 azd
->mDelays
[ti
] = AverageHrirOnset(hData
->mIrRate
, hData
->mIrPoints
, hrir
.data(), 1.0 / factor
[ti
], azd
->mDelays
[ti
]);
2923 AverageHrirMagnitude(hData
->mIrPoints
, hData
->mFftSize
, hrir
.data(), 1.0 / factor
[ti
], azd
->mIrs
[ti
]);
2925 if(!TrIsOperator(tr
, "+"))
2927 TrReadOperator(tr
, "+");
2929 if(hData
->mChannelType
== CT_STEREO
)
2931 if(azd
->mIrs
[0] == nullptr)
2933 TrErrorAt(tr
, line
, col
, "Missing left ear source reference(s).\n");
2936 else if(azd
->mIrs
[1] == nullptr)
2938 TrErrorAt(tr
, line
, col
, "Missing right ear source reference(s).\n");
2944 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2946 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2948 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2950 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2951 if(azd
->mIrs
[0] != nullptr)
2954 if(ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
)
2957 if(ei
>= hData
->mFds
[fi
].mEvCount
)
2959 TrError(tr
, "Missing source references [ %d, *, * ].\n", fi
);
2962 hData
->mFds
[fi
].mEvStart
= ei
;
2963 for(;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2965 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2967 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2969 if(azd
->mIrs
[0] == nullptr)
2971 TrError(tr
, "Missing source reference [ %d, %d, %d ].\n", fi
, ei
, ai
);
2977 for(ti
= 0;ti
< channels
;ti
++)
2979 for(fi
= 0;fi
< hData
->mFdCount
;fi
++)
2981 for(ei
= 0;ei
< hData
->mFds
[fi
].mEvCount
;ei
++)
2983 for(ai
= 0;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzCount
;ai
++)
2985 HrirAzT
*azd
= &hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
2987 azd
->mIrs
[ti
] = &hrirs
[hData
->mIrSize
* (ti
* hData
->mIrCount
+ azd
->mIndex
)];
2995 TrError(tr
, "Errant data at end of source list.\n");
2999 /* Parse the data set definition and process the source data, storing the
3000 * resulting data set as desired. If the input name is NULL it will read
3001 * from standard input.
3003 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 char *outName
)
3005 char rateStr
[8+1], expName
[MAX_PATH_LEN
];
3011 fprintf(stdout
, "Reading HRIR definition from %s...\n", inName
?inName
:"stdin");
3012 if(inName
!= nullptr)
3014 fp
= fopen(inName
, "r");
3017 fprintf(stderr
, "Error: Could not open definition file '%s'\n", inName
);
3020 TrSetup(fp
, inName
, &tr
);
3025 TrSetup(fp
, "<stdin>", &tr
);
3027 if(!ProcessMetrics(&tr
, fftSize
, truncSize
, &hData
))
3029 if(inName
!= nullptr)
3033 if(!ProcessSources(model
, &tr
, &hData
))
3043 uint c
= (hData
.mChannelType
== CT_STEREO
) ? 2 : 1;
3044 uint m
= 1 + hData
.mFftSize
/ 2;
3045 std::vector
<double> dfa(c
* m
);
3047 fprintf(stdout
, "Calculating diffuse-field average...\n");
3048 CalculateDiffuseFieldAverage(&hData
, c
, m
, surface
, limit
, dfa
.data());
3049 fprintf(stdout
, "Performing diffuse-field equalization...\n");
3050 DiffuseFieldEqualize(c
, m
, dfa
.data(), &hData
);
3052 fprintf(stdout
, "Performing minimum phase reconstruction...\n");
3053 ReconstructHrirs(&hData
);
3054 if(outRate
!= 0 && outRate
!= hData
.mIrRate
)
3056 fprintf(stdout
, "Resampling HRIRs...\n");
3057 ResampleHrirs(outRate
, &hData
);
3059 fprintf(stdout
, "Truncating minimum-phase HRIRs...\n");
3060 hData
.mIrPoints
= truncSize
;
3061 fprintf(stdout
, "Synthesizing missing elevations...\n");
3062 if(model
== HM_DATASET
)
3063 SynthesizeOnsets(&hData
);
3064 SynthesizeHrirs(&hData
);
3065 fprintf(stdout
, "Normalizing final HRIRs...\n");
3066 NormalizeHrirs(&hData
);
3067 fprintf(stdout
, "Calculating impulse delays...\n");
3068 CalculateHrtds(model
, (radius
> DEFAULT_CUSTOM_RADIUS
) ? radius
: hData
.mRadius
, &hData
);
3069 snprintf(rateStr
, 8, "%u", hData
.mIrRate
);
3070 StrSubst(outName
, "%r", rateStr
, MAX_PATH_LEN
, expName
);
3071 fprintf(stdout
, "Creating MHR data set %s...\n", expName
);
3072 ret
= StoreMhr(&hData
, expName
);
3077 static void PrintHelp(const char *argv0
, FILE *ofile
)
3079 fprintf(ofile
, "Usage: %s [<option>...]\n\n", argv0
);
3080 fprintf(ofile
, "Options:\n");
3081 fprintf(ofile
, " -m Ignored for compatibility.\n");
3082 fprintf(ofile
, " -r <rate> Change the data set sample rate to the specified value and\n");
3083 fprintf(ofile
, " resample the HRIRs accordingly.\n");
3084 fprintf(ofile
, " -f <points> Override the FFT window size (default: %u).\n", DEFAULT_FFTSIZE
);
3085 fprintf(ofile
, " -e {on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE
? "on" : "off"));
3086 fprintf(ofile
, " -s {on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE
? "on" : "off"));
3087 fprintf(ofile
, " -l {<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
3088 fprintf(ofile
, " average (default: %.2f).\n", DEFAULT_LIMIT
);
3089 fprintf(ofile
, " -w <points> Specify the size of the truncation window that's applied\n");
3090 fprintf(ofile
, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE
);
3091 fprintf(ofile
, " -d {dataset| Specify the model used for calculating the head-delay timing\n");
3092 fprintf(ofile
, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL
== HM_DATASET
) ? "dataset" : "sphere"));
3093 fprintf(ofile
, " -c <size> Use a customized head radius measured ear-to-ear in meters.\n");
3094 fprintf(ofile
, " -i <filename> Specify an HRIR definition file to use (defaults to stdin).\n");
3095 fprintf(ofile
, " -o <filename> Specify an output file. Use of '%%r' will be substituted with\n");
3096 fprintf(ofile
, " the data set sample rate.\n");
3099 // Standard command line dispatch.
3100 int main(int argc
, char *argv
[])
3102 const char *inName
= nullptr, *outName
= nullptr;
3103 uint outRate
, fftSize
;
3104 int equalize
, surface
;
3105 char *end
= nullptr;
3112 GET_UNICODE_ARGS(&argc
, &argv
);
3116 fprintf(stdout
, "HRTF Processing and Composition Utility\n\n");
3117 PrintHelp(argv
[0], stdout
);
3121 outName
= "./oalsoft_hrtf_%r.mhr";
3124 equalize
= DEFAULT_EQUALIZE
;
3125 surface
= DEFAULT_SURFACE
;
3126 limit
= DEFAULT_LIMIT
;
3127 truncSize
= DEFAULT_TRUNCSIZE
;
3128 model
= DEFAULT_HEAD_MODEL
;
3129 radius
= DEFAULT_CUSTOM_RADIUS
;
3131 while((opt
=getopt(argc
, argv
, "mr:f:e:s:l:w:d:c:e:i:o:h")) != -1)
3136 fprintf(stderr
, "Ignoring unused command '-m'.\n");
3140 outRate
= strtoul(optarg
, &end
, 10);
3141 if(end
[0] != '\0' || outRate
< MIN_RATE
|| outRate
> MAX_RATE
)
3143 fprintf(stderr
, "Error: Got unexpected value \"%s\" for option -%c, expected between %u to %u.\n", optarg
, opt
, MIN_RATE
, MAX_RATE
);
3149 fftSize
= strtoul(optarg
, &end
, 10);
3150 if(end
[0] != '\0' || (fftSize
&(fftSize
-1)) || fftSize
< MIN_FFTSIZE
|| fftSize
> MAX_FFTSIZE
)
3152 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
);
3158 if(strcmp(optarg
, "on") == 0)
3160 else if(strcmp(optarg
, "off") == 0)
3164 fprintf(stderr
, "Error: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg
, opt
);
3170 if(strcmp(optarg
, "on") == 0)
3172 else if(strcmp(optarg
, "off") == 0)
3176 fprintf(stderr
, "Error: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg
, opt
);
3182 if(strcmp(optarg
, "none") == 0)
3186 limit
= strtod(optarg
, &end
);
3187 if(end
[0] != '\0' || limit
< MIN_LIMIT
|| limit
> MAX_LIMIT
)
3189 fprintf(stderr
, "Error: Got unexpected value \"%s\" for option -%c, expected between %.0f to %.0f.\n", optarg
, opt
, MIN_LIMIT
, MAX_LIMIT
);
3196 truncSize
= strtoul(optarg
, &end
, 10);
3197 if(end
[0] != '\0' || truncSize
< MIN_TRUNCSIZE
|| truncSize
> MAX_TRUNCSIZE
|| (truncSize
%MOD_TRUNCSIZE
))
3199 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
);
3205 if(strcmp(optarg
, "dataset") == 0)
3207 else if(strcmp(optarg
, "sphere") == 0)
3211 fprintf(stderr
, "Error: Got unexpected value \"%s\" for option -%c, expected dataset or sphere.\n", optarg
, opt
);
3217 radius
= strtod(optarg
, &end
);
3218 if(end
[0] != '\0' || radius
< MIN_CUSTOM_RADIUS
|| radius
> MAX_CUSTOM_RADIUS
)
3220 fprintf(stderr
, "Error: Got unexpected value \"%s\" for option -%c, expected between %.2f to %.2f.\n", optarg
, opt
, MIN_CUSTOM_RADIUS
, MAX_CUSTOM_RADIUS
);
3234 PrintHelp(argv
[0], stdout
);
3238 PrintHelp(argv
[0], stderr
);
3243 if(!ProcessDefinition(inName
, outRate
, fftSize
, equalize
, surface
, limit
,
3244 truncSize
, model
, radius
, outName
))
3246 fprintf(stdout
, "Operation completed.\n");
3248 return EXIT_SUCCESS
;