2 * HRTF utility for producing and demonstrating the process of creating an
3 * OpenAL Soft compatible HRIR data set.
5 * Copyright (C) 2011-2019 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. This has since been modified to use coverage volume for
37 * multi-field HRIR data sets.
38 * 3. Take the diffuse-field average and limit its magnitude range.
39 * 4. Equalize the responses by using the inverse of the diffuse-field
41 * 5. Reconstruct the minimum-phase responses.
42 * 5. Zero the DC component.
43 * 6. IFFT the result and truncate to the desired-length minimum-phase FIR.
45 * The spherical head algorithm for calculating propagation delay was adapted
48 * Modeling Interaural Time Difference Assuming a Spherical Head
50 * Music 150, Musical Acoustics, Stanford University
53 * The formulae for calculating the Kaiser window metrics are from the
56 * Discrete-Time Signal Processing
57 * Alan V. Oppenheim and Ronald W. Schafer
58 * Prentice-Hall Signal Processing Series
88 #include "../getopt.h"
91 #include "alcomplex.h"
92 #include "alfstream.h"
98 #include "win_main_utf8.h"
103 using namespace std::placeholders;
108 #define M_PI (3.14159265358979323846)
112 HrirDataT::~HrirDataT() = default;
114 // Head model used for calculating the impulse delays.
117 HM_DATASET, // Measure the onset from the dataset.
118 HM_SPHERE // Calculate the onset using a spherical head model.
122 // The epsilon used to maintain signal stability.
123 #define EPSILON (1e-9)
125 // The limits to the FFT window size override on the command line.
126 #define MIN_FFTSIZE (65536)
127 #define MAX_FFTSIZE (131072)
129 // The limits to the equalization range limit on the command line.
130 #define MIN_LIMIT (2.0)
131 #define MAX_LIMIT (120.0)
133 // The limits to the truncation window size on the command line.
134 #define MIN_TRUNCSIZE (16)
135 #define MAX_TRUNCSIZE (128)
137 // The limits to the custom head radius on the command line.
138 #define MIN_CUSTOM_RADIUS (0.05)
139 #define MAX_CUSTOM_RADIUS (0.15)
141 // The defaults for the command line options.
142 #define DEFAULT_FFTSIZE (65536)
143 #define DEFAULT_EQUALIZE (1)
144 #define DEFAULT_SURFACE (1)
145 #define DEFAULT_LIMIT (24.0)
146 #define DEFAULT_TRUNCSIZE (64)
147 #define DEFAULT_HEAD_MODEL (HM_DATASET)
148 #define DEFAULT_CUSTOM_RADIUS (0.0)
150 // The maximum propagation delay value supported by OpenAL Soft.
151 #define MAX_HRTD (63.0)
153 // The OpenAL Soft HRTF format marker. It stands for minimum-phase head
154 // response protocol 03.
155 #define MHR_FORMAT ("MinPHR03")
157 /* Channel index enums. Mono uses LeftChannel only. */
158 enum ChannelIndex : uint {
164 /* Performs a string substitution. Any case-insensitive occurrences of the
165 * pattern string are replaced with the replacement string. The result is
166 * truncated if necessary.
168 static std::string StrSubst(al::span<const char> in, const al::span<const char> pat,
169 const al::span<const char> rep)
172 ret.reserve(in.size() + pat.size());
174 while(in.size() >= pat.size())
176 if(al::strncasecmp(in.data(), pat.data(), pat.size()) == 0)
178 in = in.subspan(pat.size());
179 ret.append(rep.data(), rep.size());
184 while(endpos < in.size() && in[endpos] != pat.front())
186 ret.append(in.data(), endpos);
187 in = in.subspan(endpos);
190 ret.append(in.data(), in.size());
196 /*********************
197 *** Math routines ***
198 *********************/
200 // Simple clamp routine.
201 static double Clamp(const double val, const double lower, const double upper)
203 return std::min(std::max(val, lower), upper);
206 static inline uint dither_rng(uint *seed)
208 *seed = *seed * 96314165 + 907633515;
212 // Performs a triangular probability density function dither. The input samples
213 // should be normalized (-1 to +1).
214 static void TpdfDither(double *RESTRICT out, const double *RESTRICT in, const double scale,
215 const uint count, const uint step, uint *seed)
217 static constexpr double PRNG_SCALE = 1.0 / std::numeric_limits<uint>::max();
219 for(uint i{0};i < count;i++)
221 uint prn0{dither_rng(seed)};
222 uint prn1{dither_rng(seed)};
223 *out = std::round(*(in++)*scale + (prn0*PRNG_SCALE - prn1*PRNG_SCALE));
229 /* Calculate the complex helical sequence (or discrete-time analytical signal)
230 * of the given input using the Hilbert transform. Given the natural logarithm
231 * of a signal's magnitude response, the imaginary components can be used as
232 * the angles for minimum-phase reconstruction.
234 inline static void Hilbert(const uint n, complex_d *inout)
235 { complex_hilbert({inout, n}); }
237 /* Calculate the magnitude response of the given input. This is used in
238 * place of phase decomposition, since the phase residuals are discarded for
239 * minimum phase reconstruction. The mirrored half of the response is also
242 void MagnitudeResponse(const uint n, const complex_d *in, double *out)
244 const uint m = 1 + (n / 2);
247 out[i] = std::max(std::abs(in[i]), EPSILON);
250 /* Apply a range limit (in dB) to the given magnitude response. This is used
251 * to adjust the effects of the diffuse-field average on the equalization
254 static void LimitMagnitudeResponse(const uint n, const uint m, const double limit, const double *in, double *out)
257 uint i, lower, upper;
260 halfLim = limit / 2.0;
261 // Convert the response to dB.
263 out[i] = 20.0 * std::log10(in[i]);
264 // Use six octaves to calculate the average magnitude of the signal.
265 lower = (static_cast<uint>(std::ceil(n / std::pow(2.0, 8.0)))) - 1;
266 upper = (static_cast<uint>(std::floor(n / std::pow(2.0, 2.0)))) - 1;
268 for(i = lower;i <= upper;i++)
270 ave /= upper - lower + 1;
271 // Keep the response within range of the average magnitude.
273 out[i] = Clamp(out[i], ave - halfLim, ave + halfLim);
274 // Convert the response back to linear magnitude.
276 out[i] = std::pow(10.0, out[i] / 20.0);
279 /* Reconstructs the minimum-phase component for the given magnitude response
280 * of a signal. This is equivalent to phase recomposition, sans the missing
281 * residuals (which were discarded). The mirrored half of the response is
284 static void MinimumPhase(const uint n, double *mags, complex_d *out)
286 const uint m{(n/2) + 1};
290 out[i] = std::log(mags[i]);
293 mags[i] = mags[n - i];
297 // Remove any DC offset the filter has.
301 auto a = std::exp(complex_d{0.0, out[i].imag()});
302 out[i] = a * mags[i];
307 /***************************
308 *** File storage output ***
309 ***************************/
311 // Write an ASCII string to a file.
312 static int WriteAscii(const char *out, FILE *fp, const char *filename)
317 if(fwrite(out, 1, len, fp) != len)
320 fprintf(stderr, "\nError: Bad write to file '%s'.\n", filename);
326 // Write a binary value of the given byte order and byte size to a file,
327 // loading it from a 32-bit unsigned integer.
328 static int WriteBin4(const uint bytes, const uint32_t in, FILE *fp, const char *filename)
333 for(i = 0;i < bytes;i++)
334 out[i] = (in>>(i*8)) & 0x000000FF;
336 if(fwrite(out, 1, bytes, fp) != bytes)
338 fprintf(stderr, "\nError: Bad write to file '%s'.\n", filename);
344 // Store the OpenAL Soft HRTF data set.
345 static int StoreMhr(const HrirDataT *hData, const char *filename)
347 const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
348 const uint n{hData->mIrPoints};
349 uint dither_seed{22222};
353 if((fp=fopen(filename, "wb")) == nullptr)
355 fprintf(stderr, "\nError: Could not open MHR file '%s'.\n", filename);
358 if(!WriteAscii(MHR_FORMAT, fp, filename))
360 if(!WriteBin4(4, hData->mIrRate, fp, filename))
362 if(!WriteBin4(1, static_cast<uint32_t>(hData->mChannelType), fp, filename))
364 if(!WriteBin4(1, hData->mIrPoints, fp, filename))
366 if(!WriteBin4(1, static_cast<uint>(hData->mFds.size()), fp, filename))
368 for(fi = static_cast<uint>(hData->mFds.size()-1);fi < hData->mFds.size();fi--)
370 auto fdist = static_cast<uint32_t>(std::round(1000.0 * hData->mFds[fi].mDistance));
371 if(!WriteBin4(2, fdist, fp, filename))
373 if(!WriteBin4(1, static_cast<uint32_t>(hData->mFds[fi].mEvs.size()), fp, filename))
375 for(ei = 0;ei < hData->mFds[fi].mEvs.size();ei++)
377 const auto &elev = hData->mFds[fi].mEvs[ei];
378 if(!WriteBin4(1, static_cast<uint32_t>(elev.mAzs.size()), fp, filename))
383 for(fi = static_cast<uint>(hData->mFds.size()-1);fi < hData->mFds.size();fi--)
385 constexpr double scale{8388607.0};
386 constexpr uint bps{3u};
388 for(ei = 0;ei < hData->mFds[fi].mEvs.size();ei++)
390 for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzs.size();ai++)
392 HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
393 double out[2 * MAX_TRUNCSIZE];
395 TpdfDither(out, azd->mIrs[0], scale, n, channels, &dither_seed);
396 if(hData->mChannelType == CT_STEREO)
397 TpdfDither(out+1, azd->mIrs[1], scale, n, channels, &dither_seed);
398 for(i = 0;i < (channels * n);i++)
400 const auto v = static_cast<int>(Clamp(out[i], -scale-1.0, scale));
401 if(!WriteBin4(bps, static_cast<uint32_t>(v), fp, filename))
407 for(fi = static_cast<uint>(hData->mFds.size()-1);fi < hData->mFds.size();fi--)
409 /* Delay storage has 2 bits of extra precision. */
410 constexpr double DelayPrecScale{4.0};
411 for(ei = 0;ei < hData->mFds[fi].mEvs.size();ei++)
413 for(const auto &azd : hData->mFds[fi].mEvs[ei].mAzs)
415 auto v = static_cast<uint>(std::round(azd.mDelays[0]*DelayPrecScale));
416 if(!WriteBin4(1, v, fp, filename)) return 0;
417 if(hData->mChannelType == CT_STEREO)
419 v = static_cast<uint>(std::round(azd.mDelays[1]*DelayPrecScale));
420 if(!WriteBin4(1, v, fp, filename)) return 0;
430 /***********************
431 *** HRTF processing ***
432 ***********************/
434 /* Balances the maximum HRIR magnitudes of multi-field data sets by
435 * independently normalizing each field in relation to the overall maximum.
436 * This is done to ignore distance attenuation.
438 static void BalanceFieldMagnitudes(const HrirDataT *hData, const uint channels, const uint m)
440 double maxMags[MAX_FD_COUNT];
444 for(fi = 0;fi < hData->mFds.size();fi++)
448 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
450 for(const auto &azd : hData->mFds[fi].mEvs[ei].mAzs)
452 for(ti = 0;ti < channels;ti++)
455 maxMags[fi] = std::max(azd.mIrs[ti][i], maxMags[fi]);
460 maxMag = std::max(maxMags[fi], maxMag);
463 for(fi = 0;fi < hData->mFds.size();fi++)
465 const double magFactor{maxMag / maxMags[fi]};
467 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
469 for(const auto &azd : hData->mFds[fi].mEvs[ei].mAzs)
471 for(ti = 0;ti < channels;ti++)
474 azd.mIrs[ti][i] *= magFactor;
481 /* Calculate the contribution of each HRIR to the diffuse-field average based
482 * on its coverage volume. All volumes are centered at the spherical HRIR
483 * coordinates and measured by extruded solid angle.
485 static void CalculateDfWeights(const HrirDataT *hData, double *weights)
487 double sum, innerRa, outerRa, evs, ev, upperEv, lowerEv;
488 double solidAngle, solidVolume;
492 // The head radius acts as the limit for the inner radius.
493 innerRa = hData->mRadius;
494 for(fi = 0;fi < hData->mFds.size();fi++)
496 // Each volume ends half way between progressive field measurements.
497 if((fi + 1) < hData->mFds.size())
498 outerRa = 0.5f * (hData->mFds[fi].mDistance + hData->mFds[fi + 1].mDistance);
499 // The final volume has its limit extended to some practical value.
500 // This is done to emphasize the far-field responses in the average.
504 const double raPowDiff{std::pow(outerRa, 3.0) - std::pow(innerRa, 3.0)};
505 evs = M_PI / 2.0 / static_cast<double>(hData->mFds[fi].mEvs.size() - 1);
506 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
508 const auto &elev = hData->mFds[fi].mEvs[ei];
509 // For each elevation, calculate the upper and lower limits of
511 ev = elev.mElevation;
512 lowerEv = std::max(-M_PI / 2.0, ev - evs);
513 upperEv = std::min(M_PI / 2.0, ev + evs);
514 // Calculate the surface area of the patch band.
515 solidAngle = 2.0 * M_PI * (std::sin(upperEv) - std::sin(lowerEv));
516 // Then the volume of the extruded patch band.
517 solidVolume = solidAngle * raPowDiff / 3.0;
518 // Each weight is the volume of one extruded patch.
519 weights[(fi*MAX_EV_COUNT) + ei] = solidVolume / static_cast<double>(elev.mAzs.size());
520 // Sum the total coverage volume of the HRIRs for all fields.
527 for(fi = 0;fi < hData->mFds.size();fi++)
529 // Normalize the weights given the total surface coverage for all
531 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
532 weights[(fi * MAX_EV_COUNT) + ei] /= sum;
536 /* Calculate the diffuse-field average from the given magnitude responses of
537 * the HRIR set. Weighting can be applied to compensate for the varying
538 * coverage of each HRIR. The final average can then be limited by the
539 * specified magnitude range (in positive dB; 0.0 to skip).
541 static void CalculateDiffuseFieldAverage(const HrirDataT *hData, const uint channels, const uint m,
542 const int weighted, const double limit, double *dfa)
544 std::vector<double> weights(hData->mFds.size() * MAX_EV_COUNT);
545 uint count, ti, fi, ei, i, ai;
549 // Use coverage weighting to calculate the average.
550 CalculateDfWeights(hData, weights.data());
556 // If coverage weighting is not used, the weights still need to be
557 // averaged by the number of existing HRIRs.
558 count = hData->mIrCount;
559 for(fi = 0;fi < hData->mFds.size();fi++)
561 for(ei = 0;ei < hData->mFds[fi].mEvStart;ei++)
562 count -= static_cast<uint>(hData->mFds[fi].mEvs[ei].mAzs.size());
564 weight = 1.0 / count;
566 for(fi = 0;fi < hData->mFds.size();fi++)
568 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
569 weights[(fi * MAX_EV_COUNT) + ei] = weight;
572 for(ti = 0;ti < channels;ti++)
575 dfa[(ti * m) + i] = 0.0;
576 for(fi = 0;fi < hData->mFds.size();fi++)
578 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
580 for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzs.size();ai++)
582 HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
583 // Get the weight for this HRIR's contribution.
584 double weight = weights[(fi * MAX_EV_COUNT) + ei];
586 // Add this HRIR's weighted power average to the total.
588 dfa[(ti * m) + i] += weight * azd->mIrs[ti][i] * azd->mIrs[ti][i];
592 // Finish the average calculation and keep it from being too small.
594 dfa[(ti * m) + i] = std::max(sqrt(dfa[(ti * m) + i]), EPSILON);
595 // Apply a limit to the magnitude range of the diffuse-field average
598 LimitMagnitudeResponse(hData->mFftSize, m, limit, &dfa[ti * m], &dfa[ti * m]);
602 // Perform diffuse-field equalization on the magnitude responses of the HRIR
603 // set using the given average response.
604 static void DiffuseFieldEqualize(const uint channels, const uint m, const double *dfa, const HrirDataT *hData)
608 for(fi = 0;fi < hData->mFds.size();fi++)
610 for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvs.size();ei++)
612 for(auto &azd : hData->mFds[fi].mEvs[ei].mAzs)
614 for(ti = 0;ti < channels;ti++)
617 azd.mIrs[ti][i] /= dfa[(ti * m) + i];
624 /* Given field and elevation indices and an azimuth, calculate the indices of
625 * the two HRIRs that bound the coordinate along with a factor for
626 * calculating the continuous HRIR using interpolation.
628 static void CalcAzIndices(const HrirFdT &field, const uint ei, const double az, uint *a0, uint *a1, double *af)
630 double f{(2.0*M_PI + az) * static_cast<double>(field.mEvs[ei].mAzs.size()) / (2.0*M_PI)};
631 const uint i{static_cast<uint>(f) % static_cast<uint>(field.mEvs[ei].mAzs.size())};
635 *a1 = (i + 1) % static_cast<uint>(field.mEvs[ei].mAzs.size());
639 /* Synthesize any missing onset timings at the bottom elevations of each field.
640 * This just mirrors some top elevations for the bottom, and blends the
641 * remaining elevations (not an accurate model).
643 static void SynthesizeOnsets(HrirDataT *hData)
645 const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
647 auto proc_field = [channels](HrirFdT &field) -> void
649 /* Get the starting elevation from the measurements, and use it as the
650 * upper elevation limit for what needs to be calculated.
652 const uint upperElevReal{field.mEvStart};
653 if(upperElevReal <= 0) return;
655 /* Get the lowest half of the missing elevations' delays by mirroring
656 * the top elevation delays. The responses are on a spherical grid
657 * centered between the ears, so these should align.
662 /* Take the polar opposite position of the desired measurement and
665 field.mEvs[0].mAzs[0].mDelays[0] = field.mEvs[field.mEvs.size()-1].mAzs[0].mDelays[1];
666 field.mEvs[0].mAzs[0].mDelays[1] = field.mEvs[field.mEvs.size()-1].mAzs[0].mDelays[0];
667 for(ei = 1u;ei < (upperElevReal+1)/2;++ei)
669 const uint topElev{static_cast<uint>(field.mEvs.size()-ei-1)};
671 for(uint ai{0u};ai < field.mEvs[ei].mAzs.size();ai++)
676 /* Rotate this current azimuth by a half-circle, and lookup
677 * the mirrored elevation to find the indices for the polar
678 * opposite position (may need blending).
680 const double az{field.mEvs[ei].mAzs[ai].mAzimuth + M_PI};
681 CalcAzIndices(field, topElev, az, &a0, &a1, &af);
683 /* Blend the delays, and again, swap the ears. */
684 field.mEvs[ei].mAzs[ai].mDelays[0] = Lerp(
685 field.mEvs[topElev].mAzs[a0].mDelays[1],
686 field.mEvs[topElev].mAzs[a1].mDelays[1], af);
687 field.mEvs[ei].mAzs[ai].mDelays[1] = Lerp(
688 field.mEvs[topElev].mAzs[a0].mDelays[0],
689 field.mEvs[topElev].mAzs[a1].mDelays[0], af);
695 field.mEvs[0].mAzs[0].mDelays[0] = field.mEvs[field.mEvs.size()-1].mAzs[0].mDelays[0];
696 for(ei = 1u;ei < (upperElevReal+1)/2;++ei)
698 const uint topElev{static_cast<uint>(field.mEvs.size()-ei-1)};
700 for(uint ai{0u};ai < field.mEvs[ei].mAzs.size();ai++)
705 /* For mono data sets, mirror the azimuth front<->back
706 * since the other ear is a mirror of what we have (e.g.
707 * the left ear's back-left is simulated with the right
708 * ear's front-right, which uses the left ear's front-left
711 double az{field.mEvs[ei].mAzs[ai].mAzimuth};
712 if(az <= M_PI) az = M_PI - az;
713 else az = (M_PI*2.0)-az + M_PI;
714 CalcAzIndices(field, topElev, az, &a0, &a1, &af);
716 field.mEvs[ei].mAzs[ai].mDelays[0] = Lerp(
717 field.mEvs[topElev].mAzs[a0].mDelays[0],
718 field.mEvs[topElev].mAzs[a1].mDelays[0], af);
722 /* Record the lowest elevation filled in with the mirrored top. */
723 const uint lowerElevFake{ei-1u};
725 /* Fill in the remaining delays using bilinear interpolation. This
726 * helps smooth the transition back to the real delays.
728 for(;ei < upperElevReal;++ei)
730 const double ef{(field.mEvs[upperElevReal].mElevation - field.mEvs[ei].mElevation) /
731 (field.mEvs[upperElevReal].mElevation - field.mEvs[lowerElevFake].mElevation)};
733 for(uint ai{0u};ai < field.mEvs[ei].mAzs.size();ai++)
738 double az{field.mEvs[ei].mAzs[ai].mAzimuth};
739 CalcAzIndices(field, upperElevReal, az, &a0, &a1, &af0);
740 CalcAzIndices(field, lowerElevFake, az, &a2, &a3, &af1);
742 (1.0-ef) * (1.0-af0),
748 for(uint ti{0u};ti < channels;ti++)
750 field.mEvs[ei].mAzs[ai].mDelays[ti] =
751 field.mEvs[upperElevReal].mAzs[a0].mDelays[ti]*blend[0] +
752 field.mEvs[upperElevReal].mAzs[a1].mDelays[ti]*blend[1] +
753 field.mEvs[lowerElevFake].mAzs[a2].mDelays[ti]*blend[2] +
754 field.mEvs[lowerElevFake].mAzs[a3].mDelays[ti]*blend[3];
759 std::for_each(hData->mFds.begin(), hData->mFds.end(), proc_field);
762 /* Attempt to synthesize any missing HRIRs at the bottom elevations of each
763 * field. Right now this just blends the lowest elevation HRIRs together and
764 * applies a low-pass filter to simulate body occlusion. It is a simple, if
767 static void SynthesizeHrirs(HrirDataT *hData)
769 const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
770 auto htemp = std::vector<complex_d>(hData->mFftSize);
771 const uint m{hData->mFftSize/2u + 1u};
772 auto filter = std::vector<double>(m);
773 const double beta{3.5e-6 * hData->mIrRate};
775 auto proc_field = [channels,m,beta,&htemp,&filter](HrirFdT &field) -> void
777 const uint oi{field.mEvStart};
780 for(uint ti{0u};ti < channels;ti++)
785 /* Use the lowest immediate-left response for the left ear and
786 * lowest immediate-right response for the right ear. Given no comb
787 * effects as a result of the left response reaching the right ear
788 * and vice-versa, this produces a decent phantom-center response
789 * underneath the head.
791 CalcAzIndices(field, oi, ((ti==0) ? -M_PI : M_PI) / 2.0, &a0, &a1, &af);
792 for(uint i{0u};i < m;i++)
794 field.mEvs[0].mAzs[0].mIrs[ti][i] = Lerp(field.mEvs[oi].mAzs[a0].mIrs[ti][i],
795 field.mEvs[oi].mAzs[a1].mIrs[ti][i], af);
799 for(uint ei{1u};ei < field.mEvStart;ei++)
801 const double of{static_cast<double>(ei) / field.mEvStart};
802 const double b{(1.0 - of) * beta};
805 /* Calculate a low-pass filter to simulate body occlusion. */
806 lp[0] = Lerp(1.0, lp[0], b);
807 lp[1] = Lerp(lp[0], lp[1], b);
808 lp[2] = Lerp(lp[1], lp[2], b);
809 lp[3] = Lerp(lp[2], lp[3], b);
811 for(size_t i{1u};i < htemp.size();i++)
813 lp[0] = Lerp(0.0, lp[0], b);
814 lp[1] = Lerp(lp[0], lp[1], b);
815 lp[2] = Lerp(lp[1], lp[2], b);
816 lp[3] = Lerp(lp[2], lp[3], b);
819 /* Get the filter's frequency-domain response and extract the
820 * frequency magnitudes (phase will be reconstructed later)).
822 FftForward(static_cast<uint>(htemp.size()), htemp.data());
823 std::transform(htemp.cbegin(), htemp.cbegin()+m, filter.begin(),
824 [](const complex_d &c) -> double { return std::abs(c); });
826 for(uint ai{0u};ai < field.mEvs[ei].mAzs.size();ai++)
831 CalcAzIndices(field, oi, field.mEvs[ei].mAzs[ai].mAzimuth, &a0, &a1, &af);
832 for(uint ti{0u};ti < channels;ti++)
834 for(uint i{0u};i < m;i++)
836 /* Blend the two defined HRIRs closest to this azimuth,
837 * then blend that with the synthesized -90 elevation.
839 const double s1{Lerp(field.mEvs[oi].mAzs[a0].mIrs[ti][i],
840 field.mEvs[oi].mAzs[a1].mIrs[ti][i], af)};
841 const double s{Lerp(field.mEvs[0].mAzs[0].mIrs[ti][i], s1, of)};
842 field.mEvs[ei].mAzs[ai].mIrs[ti][i] = s * filter[i];
847 const double b{beta};
849 lp[0] = Lerp(1.0, lp[0], b);
850 lp[1] = Lerp(lp[0], lp[1], b);
851 lp[2] = Lerp(lp[1], lp[2], b);
852 lp[3] = Lerp(lp[2], lp[3], b);
854 for(size_t i{1u};i < htemp.size();i++)
856 lp[0] = Lerp(0.0, lp[0], b);
857 lp[1] = Lerp(lp[0], lp[1], b);
858 lp[2] = Lerp(lp[1], lp[2], b);
859 lp[3] = Lerp(lp[2], lp[3], b);
862 FftForward(static_cast<uint>(htemp.size()), htemp.data());
863 std::transform(htemp.cbegin(), htemp.cbegin()+m, filter.begin(),
864 [](const complex_d &c) -> double { return std::abs(c); });
866 for(uint ti{0u};ti < channels;ti++)
868 for(uint i{0u};i < m;i++)
869 field.mEvs[0].mAzs[0].mIrs[ti][i] *= filter[i];
872 std::for_each(hData->mFds.begin(), hData->mFds.end(), proc_field);
875 // The following routines assume a full set of HRIRs for all elevations.
877 /* Perform minimum-phase reconstruction using the magnitude responses of the
878 * HRIR set. Work is delegated to this struct, which runs asynchronously on one
879 * or more threads (sharing the same reconstructor object).
881 struct HrirReconstructor {
882 std::vector<double*> mIrs;
883 std::atomic<size_t> mCurrent;
884 std::atomic<size_t> mDone;
890 auto h = std::vector<complex_d>(mFftSize);
891 auto mags = std::vector<double>(mFftSize);
892 size_t m{(mFftSize/2) + 1};
896 /* Load the current index to process. */
897 size_t idx{mCurrent.load()};
899 /* If the index is at the end, we're done. */
900 if(idx >= mIrs.size())
902 /* Otherwise, increment the current index atomically so other
903 * threads know to go to the next one. If this call fails, the
904 * current index was just changed by another thread and the new
905 * value is loaded into idx, which we'll recheck.
907 } while(!mCurrent.compare_exchange_weak(idx, idx+1, std::memory_order_relaxed));
909 /* Now do the reconstruction, and apply the inverse FFT to get the
910 * time-domain response.
912 for(size_t i{0};i < m;++i)
913 mags[i] = std::max(mIrs[idx][i], EPSILON);
914 MinimumPhase(mFftSize, mags.data(), h.data());
915 FftInverse(mFftSize, h.data());
916 for(uint i{0u};i < mIrPoints;++i)
917 mIrs[idx][i] = h[i].real();
919 /* Increment the number of IRs done. */
925 static void ReconstructHrirs(const HrirDataT *hData, const uint numThreads)
927 const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
929 /* Set up the reconstructor with the needed size info and pointers to the
932 HrirReconstructor reconstructor;
933 reconstructor.mCurrent.store(0, std::memory_order_relaxed);
934 reconstructor.mDone.store(0, std::memory_order_relaxed);
935 reconstructor.mFftSize = hData->mFftSize;
936 reconstructor.mIrPoints = hData->mIrPoints;
937 for(const auto &field : hData->mFds)
939 for(auto &elev : field.mEvs)
941 for(const auto &azd : elev.mAzs)
943 for(uint ti{0u};ti < channels;ti++)
944 reconstructor.mIrs.push_back(azd.mIrs[ti]);
949 /* Launch threads to work on reconstruction. */
950 std::vector<std::thread> thrds;
951 thrds.reserve(numThreads);
952 for(size_t i{0};i < numThreads;++i)
953 thrds.emplace_back(std::mem_fn(&HrirReconstructor::Worker), &reconstructor);
955 /* Keep track of the number of IRs done, periodically reporting it. */
958 std::this_thread::sleep_for(std::chrono::milliseconds{50});
960 count = reconstructor.mDone.load();
961 size_t pcdone{count * 100 / reconstructor.mIrs.size()};
963 printf("\r%3zu%% done (%zu of %zu)", pcdone, count, reconstructor.mIrs.size());
965 } while(count < reconstructor.mIrs.size());
968 for(auto &thrd : thrds)
975 // Normalize the HRIR set and slightly attenuate the result.
976 static void NormalizeHrirs(HrirDataT *hData)
978 const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
979 const uint irSize{hData->mIrPoints};
981 /* Find the maximum amplitude and RMS out of all the IRs. */
982 struct LevelPair { double amp, rms; };
983 auto mesasure_channel = [irSize](const LevelPair levels, const double *ir)
985 /* Calculate the peak amplitude and RMS of this IR. */
986 auto current = std::accumulate(ir, ir+irSize, LevelPair{0.0, 0.0},
987 [](const LevelPair cur, const double impulse)
989 return LevelPair{std::max(std::abs(impulse), cur.amp), cur.rms + impulse*impulse};
991 current.rms = std::sqrt(current.rms / irSize);
993 /* Accumulate levels by taking the maximum amplitude and RMS. */
994 return LevelPair{std::max(current.amp, levels.amp), std::max(current.rms, levels.rms)};
996 auto measure_azi = [channels,mesasure_channel](const LevelPair levels, const HrirAzT &azi)
997 { return std::accumulate(azi.mIrs, azi.mIrs+channels, levels, mesasure_channel); };
998 auto measure_elev = [measure_azi](const LevelPair levels, const HrirEvT &elev)
999 { return std::accumulate(elev.mAzs.cbegin(), elev.mAzs.cend(), levels, measure_azi); };
1000 auto measure_field = [measure_elev](const LevelPair levels, const HrirFdT &field)
1001 { return std::accumulate(field.mEvs.cbegin(), field.mEvs.cend(), levels, measure_elev); };
1003 const auto maxlev = std::accumulate(hData->mFds.begin(), hData->mFds.end(),
1004 LevelPair{0.0, 0.0}, measure_field);
1006 /* Normalize using the maximum RMS of the HRIRs. The RMS measure for the
1007 * non-filtered signal is of an impulse with equal length (to the filter):
1009 * rms_impulse = sqrt(sum([ 1^2, 0^2, 0^2, ... ]) / n)
1012 * This helps keep a more consistent volume between the non-filtered signal
1013 * and various data sets.
1015 double factor{std::sqrt(1.0 / irSize) / maxlev.rms};
1017 /* Also ensure the samples themselves won't clip. */
1018 factor = std::min(factor, 0.99/maxlev.amp);
1020 /* Now scale all IRs by the given factor. */
1021 auto proc_channel = [irSize,factor](double *ir)
1022 { std::transform(ir, ir+irSize, ir, [factor](double s){ return s * factor; }); };
1023 auto proc_azi = [channels,proc_channel](HrirAzT &azi)
1024 { std::for_each(azi.mIrs, azi.mIrs+channels, proc_channel); };
1025 auto proc_elev = [proc_azi](HrirEvT &elev)
1026 { std::for_each(elev.mAzs.begin(), elev.mAzs.end(), proc_azi); };
1027 auto proc1_field = [proc_elev](HrirFdT &field)
1028 { std::for_each(field.mEvs.begin(), field.mEvs.end(), proc_elev); };
1030 std::for_each(hData->mFds.begin(), hData->mFds.end(), proc1_field);
1033 // Calculate the left-ear time delay using a spherical head model.
1034 static double CalcLTD(const double ev, const double az, const double rad, const double dist)
1036 double azp, dlp, l, al;
1038 azp = std::asin(std::cos(ev) * std::sin(az));
1039 dlp = std::sqrt((dist*dist) + (rad*rad) + (2.0*dist*rad*sin(azp)));
1040 l = std::sqrt((dist*dist) - (rad*rad));
1041 al = (0.5 * M_PI) + azp;
1043 dlp = l + (rad * (al - std::acos(rad / dist)));
1047 // Calculate the effective head-related time delays for each minimum-phase
1048 // HRIR. This is done per-field since distance delay is ignored.
1049 static void CalculateHrtds(const HeadModelT model, const double radius, HrirDataT *hData)
1051 uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
1052 double customRatio{radius / hData->mRadius};
1055 if(model == HM_SPHERE)
1057 for(auto &field : hData->mFds)
1059 for(auto &elev : field.mEvs)
1061 for(auto &azd : elev.mAzs)
1063 for(ti = 0;ti < channels;ti++)
1064 azd.mDelays[ti] = CalcLTD(elev.mElevation, azd.mAzimuth, radius, field.mDistance);
1069 else if(customRatio != 1.0)
1071 for(auto &field : hData->mFds)
1073 for(auto &elev : field.mEvs)
1075 for(auto &azd : elev.mAzs)
1077 for(ti = 0;ti < channels;ti++)
1078 azd.mDelays[ti] *= customRatio;
1084 double maxHrtd{0.0};
1085 for(auto &field : hData->mFds)
1087 double minHrtd{std::numeric_limits<double>::infinity()};
1088 for(auto &elev : field.mEvs)
1090 for(auto &azd : elev.mAzs)
1092 for(ti = 0;ti < channels;ti++)
1093 minHrtd = std::min(azd.mDelays[ti], minHrtd);
1097 for(auto &elev : field.mEvs)
1099 for(auto &azd : elev.mAzs)
1101 for(ti = 0;ti < channels;ti++)
1103 azd.mDelays[ti] = (azd.mDelays[ti]-minHrtd) * hData->mIrRate;
1104 maxHrtd = std::max(maxHrtd, azd.mDelays[ti]);
1109 if(maxHrtd > MAX_HRTD)
1111 fprintf(stdout, " Scaling for max delay of %f samples to %f\n...\n", maxHrtd, MAX_HRTD);
1112 const double scale{MAX_HRTD / maxHrtd};
1113 for(auto &field : hData->mFds)
1115 for(auto &elev : field.mEvs)
1117 for(auto &azd : elev.mAzs)
1119 for(ti = 0;ti < channels;ti++)
1120 azd.mDelays[ti] *= scale;
1127 // Allocate and configure dynamic HRIR structures.
1128 bool PrepareHrirData(const al::span<const double> distances,
1129 const al::span<const uint,MAX_FD_COUNT> evCounts,
1130 const al::span<const std::array<uint,MAX_EV_COUNT>,MAX_FD_COUNT> azCounts, HrirDataT *hData)
1132 uint evTotal{0}, azTotal{0};
1134 for(size_t fi{0};fi < distances.size();++fi)
1136 evTotal += evCounts[fi];
1137 for(size_t ei{0};ei < evCounts[fi];++ei)
1138 azTotal += azCounts[fi][ei];
1140 if(!evTotal || !azTotal)
1143 hData->mEvsBase.resize(evTotal);
1144 hData->mAzsBase.resize(azTotal);
1145 hData->mFds.resize(distances.size());
1146 hData->mIrCount = azTotal;
1149 for(size_t fi{0};fi < distances.size();++fi)
1151 hData->mFds[fi].mDistance = distances[fi];
1152 hData->mFds[fi].mEvStart = 0;
1153 hData->mFds[fi].mEvs = {&hData->mEvsBase[evTotal], evCounts[fi]};
1154 evTotal += evCounts[fi];
1155 for(uint ei{0};ei < evCounts[fi];++ei)
1157 uint azCount = azCounts[fi][ei];
1159 hData->mFds[fi].mEvs[ei].mElevation = -M_PI / 2.0 + M_PI * ei / (evCounts[fi] - 1);
1160 hData->mFds[fi].mEvs[ei].mAzs = {&hData->mAzsBase[azTotal], azCount};
1161 for(uint ai{0};ai < azCount;ai++)
1163 hData->mFds[fi].mEvs[ei].mAzs[ai].mAzimuth = 2.0 * M_PI * ai / azCount;
1164 hData->mFds[fi].mEvs[ei].mAzs[ai].mIndex = azTotal + ai;
1165 hData->mFds[fi].mEvs[ei].mAzs[ai].mDelays[0] = 0.0;
1166 hData->mFds[fi].mEvs[ei].mAzs[ai].mDelays[1] = 0.0;
1167 hData->mFds[fi].mEvs[ei].mAzs[ai].mIrs[0] = nullptr;
1168 hData->mFds[fi].mEvs[ei].mAzs[ai].mIrs[1] = nullptr;
1177 /* Parse the data set definition and process the source data, storing the
1178 * resulting data set as desired. If the input name is NULL it will read
1179 * from standard input.
1181 static int ProcessDefinition(const char *inName, const uint outRate, const ChannelModeT chanMode,
1182 const bool farfield, const uint numThreads, const uint fftSize, const int equalize,
1183 const int surface, const double limit, const uint truncSize, const HeadModelT model,
1184 const double radius, const char *outName)
1188 fprintf(stdout, "Using %u thread%s.\n", numThreads, (numThreads==1)?"":"s");
1192 fprintf(stdout, "Reading HRIR definition from %s...\n", inName);
1193 if(!LoadDefInput(std::cin, nullptr, 0, inName, fftSize, truncSize, outRate, chanMode, &hData))
1198 std::unique_ptr<al::ifstream> input{new al::ifstream{inName}};
1199 if(!input->is_open())
1201 fprintf(stderr, "Error: Could not open input file '%s'\n", inName);
1205 char startbytes[4]{};
1206 input->read(startbytes, sizeof(startbytes));
1207 std::streamsize startbytecount{input->gcount()};
1208 if(startbytecount != sizeof(startbytes) || !input->good())
1210 fprintf(stderr, "Error: Could not read input file '%s'\n", inName);
1214 if(startbytes[0] == '\x89' && startbytes[1] == 'H' && startbytes[2] == 'D'
1215 && startbytes[3] == 'F')
1218 fprintf(stdout, "Reading HRTF data from %s...\n", inName);
1219 if(!LoadSofaFile(inName, numThreads, fftSize, truncSize, outRate, chanMode, &hData))
1224 fprintf(stdout, "Reading HRIR definition from %s...\n", inName);
1225 if(!LoadDefInput(*input, startbytes, startbytecount, inName, fftSize, truncSize, outRate, chanMode, &hData))
1232 uint c{(hData.mChannelType == CT_STEREO) ? 2u : 1u};
1233 uint m{hData.mFftSize/2u + 1u};
1234 auto dfa = std::vector<double>(c * m);
1236 if(hData.mFds.size() > 1)
1238 fprintf(stdout, "Balancing field magnitudes...\n");
1239 BalanceFieldMagnitudes(&hData, c, m);
1241 fprintf(stdout, "Calculating diffuse-field average...\n");
1242 CalculateDiffuseFieldAverage(&hData, c, m, surface, limit, dfa.data());
1243 fprintf(stdout, "Performing diffuse-field equalization...\n");
1244 DiffuseFieldEqualize(c, m, dfa.data(), &hData);
1246 if(hData.mFds.size() > 1)
1248 fprintf(stdout, "Sorting %zu fields...\n", hData.mFds.size());
1249 std::sort(hData.mFds.begin(), hData.mFds.end(),
1250 [](const HrirFdT &lhs, const HrirFdT &rhs) noexcept
1251 { return lhs.mDistance < rhs.mDistance; });
1254 fprintf(stdout, "Clearing %zu near field%s...\n", hData.mFds.size()-1,
1255 (hData.mFds.size()-1 != 1) ? "s" : "");
1256 hData.mFds.erase(hData.mFds.cbegin(), hData.mFds.cend()-1);
1259 fprintf(stdout, "Synthesizing missing elevations...\n");
1260 if(model == HM_DATASET)
1261 SynthesizeOnsets(&hData);
1262 SynthesizeHrirs(&hData);
1263 fprintf(stdout, "Performing minimum phase reconstruction...\n");
1264 ReconstructHrirs(&hData, numThreads);
1265 fprintf(stdout, "Truncating minimum-phase HRIRs...\n");
1266 hData.mIrPoints = truncSize;
1267 fprintf(stdout, "Normalizing final HRIRs...\n");
1268 NormalizeHrirs(&hData);
1269 fprintf(stdout, "Calculating impulse delays...\n");
1270 CalculateHrtds(model, (radius > DEFAULT_CUSTOM_RADIUS) ? radius : hData.mRadius, &hData);
1272 const auto rateStr = std::to_string(hData.mIrRate);
1273 const auto expName = StrSubst({outName, strlen(outName)}, {"%r", 2},
1274 {rateStr.data(), rateStr.size()});
1275 fprintf(stdout, "Creating MHR data set %s...\n", expName.c_str());
1276 return StoreMhr(&hData, expName.c_str());
1279 static void PrintHelp(const char *argv0, FILE *ofile)
1281 fprintf(ofile, "Usage: %s [<option>...]\n\n", argv0);
1282 fprintf(ofile, "Options:\n");
1283 fprintf(ofile, " -r <rate> Change the data set sample rate to the specified value and\n");
1284 fprintf(ofile, " resample the HRIRs accordingly.\n");
1285 fprintf(ofile, " -m Change the data set to mono, mirroring the left ear for the\n");
1286 fprintf(ofile, " right ear.\n");
1287 fprintf(ofile, " -a Change the data set to single field, using the farthest field.\n");
1288 fprintf(ofile, " -j <threads> Number of threads used to process HRIRs (default: 2).\n");
1289 fprintf(ofile, " -f <points> Override the FFT window size (default: %u).\n", DEFAULT_FFTSIZE);
1290 fprintf(ofile, " -e {on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE ? "on" : "off"));
1291 fprintf(ofile, " -s {on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE ? "on" : "off"));
1292 fprintf(ofile, " -l {<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
1293 fprintf(ofile, " average (default: %.2f).\n", DEFAULT_LIMIT);
1294 fprintf(ofile, " -w <points> Specify the size of the truncation window that's applied\n");
1295 fprintf(ofile, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE);
1296 fprintf(ofile, " -d {dataset| Specify the model used for calculating the head-delay timing\n");
1297 fprintf(ofile, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL == HM_DATASET) ? "dataset" : "sphere"));
1298 fprintf(ofile, " -c <radius> Use a customized head radius measured to-ear in meters.\n");
1299 fprintf(ofile, " -i <filename> Specify an HRIR definition file to use (defaults to stdin).\n");
1300 fprintf(ofile, " -o <filename> Specify an output file. Use of '%%r' will be substituted with\n");
1301 fprintf(ofile, " the data set sample rate.\n");
1304 // Standard command line dispatch.
1305 int main(int argc, char *argv[])
1307 const char *inName = nullptr, *outName = nullptr;
1308 uint outRate, fftSize;
1309 int equalize, surface;
1310 char *end = nullptr;
1311 ChannelModeT chanMode;
1322 fprintf(stdout, "HRTF Processing and Composition Utility\n\n");
1323 PrintHelp(argv[0], stdout);
1327 outName = "./oalsoft_hrtf_%r.mhr";
1329 chanMode = CM_AllowStereo;
1330 fftSize = DEFAULT_FFTSIZE;
1331 equalize = DEFAULT_EQUALIZE;
1332 surface = DEFAULT_SURFACE;
1333 limit = DEFAULT_LIMIT;
1335 truncSize = DEFAULT_TRUNCSIZE;
1336 model = DEFAULT_HEAD_MODEL;
1337 radius = DEFAULT_CUSTOM_RADIUS;
1340 while((opt=getopt(argc, argv, "r:maj:f:e:s:l:w:d:c:e:i:o:h")) != -1)
1345 outRate = static_cast<uint>(strtoul(optarg, &end, 10));
1346 if(end[0] != '\0' || outRate < MIN_RATE || outRate > MAX_RATE)
1348 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected between %u to %u.\n", optarg, opt, MIN_RATE, MAX_RATE);
1354 chanMode = CM_ForceMono;
1362 numThreads = static_cast<uint>(strtoul(optarg, &end, 10));
1363 if(end[0] != '\0' || numThreads > 64)
1365 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected between %u to %u.\n", optarg, opt, 0, 64);
1369 numThreads = std::thread::hardware_concurrency();
1373 fftSize = static_cast<uint>(strtoul(optarg, &end, 10));
1374 if(end[0] != '\0' || (fftSize&(fftSize-1)) || fftSize < MIN_FFTSIZE || fftSize > MAX_FFTSIZE)
1376 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected a power-of-two between %u to %u.\n", optarg, opt, MIN_FFTSIZE, MAX_FFTSIZE);
1382 if(strcmp(optarg, "on") == 0)
1384 else if(strcmp(optarg, "off") == 0)
1388 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg, opt);
1394 if(strcmp(optarg, "on") == 0)
1396 else if(strcmp(optarg, "off") == 0)
1400 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg, opt);
1406 if(strcmp(optarg, "none") == 0)
1410 limit = strtod(optarg, &end);
1411 if(end[0] != '\0' || limit < MIN_LIMIT || limit > MAX_LIMIT)
1413 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected between %.0f to %.0f.\n", optarg, opt, MIN_LIMIT, MAX_LIMIT);
1420 truncSize = static_cast<uint>(strtoul(optarg, &end, 10));
1421 if(end[0] != '\0' || truncSize < MIN_TRUNCSIZE || truncSize > MAX_TRUNCSIZE)
1423 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected between %u to %u.\n", optarg, opt, MIN_TRUNCSIZE, MAX_TRUNCSIZE);
1429 if(strcmp(optarg, "dataset") == 0)
1431 else if(strcmp(optarg, "sphere") == 0)
1435 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected dataset or sphere.\n", optarg, opt);
1441 radius = strtod(optarg, &end);
1442 if(end[0] != '\0' || radius < MIN_CUSTOM_RADIUS || radius > MAX_CUSTOM_RADIUS)
1444 fprintf(stderr, "\nError: Got unexpected value \"%s\" for option -%c, expected between %.2f to %.2f.\n", optarg, opt, MIN_CUSTOM_RADIUS, MAX_CUSTOM_RADIUS);
1458 PrintHelp(argv[0], stdout);
1462 PrintHelp(argv[0], stderr);
1467 int ret = ProcessDefinition(inName, outRate, chanMode, farfield, numThreads, fftSize, equalize,
1468 surface, limit, truncSize, model, radius, outName);
1470 fprintf(stdout, "Operation completed.\n");
1472 return EXIT_SUCCESS;