2 * HRTF utility for producing and demonstrating the process of creating an
3 * OpenAL Soft compatible HRIR data set.
5 * Copyright (C) 2018-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
42 #include "alnumeric.h"
44 #include "polyphase_resampler.h"
45 #include "sofa-support.h"
50 using uint
= unsigned int;
52 /* Attempts to produce a compatible layout. Most data sets tend to be
53 * uniform and have the same major axis as used by OpenAL Soft's HRTF model.
54 * This will remove outliers and produce a maximally dense layout when
55 * possible. Those sets that contain purely random measurements or use
56 * different major axes will fail.
58 static bool PrepareLayout(const uint m
, const float *xyzs
, HrirDataT
*hData
)
60 fprintf(stdout
, "Detecting compatible layout...\n");
62 auto fds
= GetCompatibleLayout(m
, xyzs
);
63 if(fds
.size() > MAX_FD_COUNT
)
65 fprintf(stdout
, "Incompatible layout (inumerable radii).\n");
69 std::array
<double,MAX_FD_COUNT
> distances
{};
70 std::array
<uint
,MAX_FD_COUNT
> evCounts
{};
71 auto azCounts
= std::vector
<std::array
<uint
,MAX_EV_COUNT
>>(MAX_FD_COUNT
);
72 for(auto &azs
: azCounts
) azs
.fill(0u);
74 uint fi
{0u}, ir_total
{0u};
75 for(const auto &field
: fds
)
77 distances
[fi
] = field
.mDistance
;
78 evCounts
[fi
] = field
.mEvCount
;
80 for(uint ei
{0u};ei
< field
.mEvStart
;ei
++)
81 azCounts
[fi
][ei
] = field
.mAzCounts
[field
.mEvCount
-ei
-1];
82 for(uint ei
{field
.mEvStart
};ei
< field
.mEvCount
;ei
++)
84 azCounts
[fi
][ei
] = field
.mAzCounts
[ei
];
85 ir_total
+= field
.mAzCounts
[ei
];
90 fprintf(stdout
, "Using %u of %u IRs.\n", ir_total
, m
);
91 const auto azs
= al::span
{azCounts
}.first
<MAX_FD_COUNT
>();
92 return PrepareHrirData(al::span
{distances
}.first(fi
), evCounts
, azs
, hData
);
96 float GetSampleRate(MYSOFA_HRTF
*sofaHrtf
)
98 const char *srate_dim
{nullptr};
99 const char *srate_units
{nullptr};
100 MYSOFA_ARRAY
*srate_array
{&sofaHrtf
->DataSamplingRate
};
101 MYSOFA_ATTRIBUTE
*srate_attrs
{srate_array
->attributes
};
104 if(std::string
{"DIMENSION_LIST"} == srate_attrs
->name
)
108 fprintf(stderr
, "Duplicate SampleRate.DIMENSION_LIST\n");
111 srate_dim
= srate_attrs
->value
;
113 else if(std::string
{"Units"} == srate_attrs
->name
)
117 fprintf(stderr
, "Duplicate SampleRate.Units\n");
120 srate_units
= srate_attrs
->value
;
123 fprintf(stderr
, "Unexpected sample rate attribute: %s = %s\n", srate_attrs
->name
,
125 srate_attrs
= srate_attrs
->next
;
129 fprintf(stderr
, "Missing sample rate dimensions\n");
132 if(srate_dim
!= std::string
{"I"})
134 fprintf(stderr
, "Unsupported sample rate dimensions: %s\n", srate_dim
);
139 fprintf(stderr
, "Missing sample rate unit type\n");
142 if(srate_units
!= std::string
{"hertz"})
144 fprintf(stderr
, "Unsupported sample rate unit type: %s\n", srate_units
);
147 /* I dimensions guarantees 1 element, so just extract it. */
148 if(srate_array
->values
[0] < float{MIN_RATE
} || srate_array
->values
[0] > float{MAX_RATE
})
150 fprintf(stderr
, "Sample rate out of range: %f (expected %u to %u)", srate_array
->values
[0],
154 return srate_array
->values
[0];
157 enum class DelayType
: uint8_t {
159 I_R
, /* [1][Channels] */
160 M_R
, /* [HRIRs][Channels] */
163 DelayType
PrepareDelay(MYSOFA_HRTF
*sofaHrtf
)
165 const char *delay_dim
{nullptr};
166 MYSOFA_ARRAY
*delay_array
{&sofaHrtf
->DataDelay
};
167 MYSOFA_ATTRIBUTE
*delay_attrs
{delay_array
->attributes
};
170 if(std::string
{"DIMENSION_LIST"} == delay_attrs
->name
)
174 fprintf(stderr
, "Duplicate Delay.DIMENSION_LIST\n");
175 return DelayType::Invalid
;
177 delay_dim
= delay_attrs
->value
;
180 fprintf(stderr
, "Unexpected delay attribute: %s = %s\n", delay_attrs
->name
,
181 delay_attrs
->value
? delay_attrs
->value
: "<null>");
182 delay_attrs
= delay_attrs
->next
;
186 fprintf(stderr
, "Missing delay dimensions\n");
187 return DelayType::None
;
189 if(delay_dim
== std::string
{"I,R"})
190 return DelayType::I_R
;
191 if(delay_dim
== std::string
{"M,R"})
192 return DelayType::M_R
;
194 fprintf(stderr
, "Unsupported delay dimensions: %s\n", delay_dim
);
195 return DelayType::Invalid
;
198 bool CheckIrData(MYSOFA_HRTF
*sofaHrtf
)
200 const char *ir_dim
{nullptr};
201 MYSOFA_ARRAY
*ir_array
{&sofaHrtf
->DataIR
};
202 MYSOFA_ATTRIBUTE
*ir_attrs
{ir_array
->attributes
};
205 if(std::string
{"DIMENSION_LIST"} == ir_attrs
->name
)
209 fprintf(stderr
, "Duplicate IR.DIMENSION_LIST\n");
212 ir_dim
= ir_attrs
->value
;
215 fprintf(stderr
, "Unexpected IR attribute: %s = %s\n", ir_attrs
->name
,
216 ir_attrs
->value
? ir_attrs
->value
: "<null>");
217 ir_attrs
= ir_attrs
->next
;
221 fprintf(stderr
, "Missing IR dimensions\n");
224 if(ir_dim
!= std::string
{"M,R,N"})
226 fprintf(stderr
, "Unsupported IR dimensions: %s\n", ir_dim
);
233 /* Calculate the onset time of a HRIR. */
234 static constexpr int OnsetRateMultiple
{10};
235 static double CalcHrirOnset(PPhaseResampler
&rs
, const uint rate
, const uint n
,
236 al::span
<double> upsampled
, const double *hrir
)
238 rs
.process(n
, hrir
, static_cast<uint
>(upsampled
.size()), upsampled
.data());
240 auto abs_lt
= [](const double &lhs
, const double &rhs
) -> bool
241 { return std::abs(lhs
) < std::abs(rhs
); };
242 auto iter
= std::max_element(upsampled
.cbegin(), upsampled
.cend(), abs_lt
);
243 return static_cast<double>(std::distance(upsampled
.cbegin(), iter
)) /
244 (double{OnsetRateMultiple
}*rate
);
247 /* Calculate the magnitude response of a HRIR. */
248 static void CalcHrirMagnitude(const uint points
, const uint n
, al::span
<complex_d
> h
, double *hrir
)
250 auto iter
= std::copy_n(hrir
, points
, h
.begin());
251 std::fill(iter
, h
.end(), complex_d
{0.0, 0.0});
253 FftForward(n
, h
.data());
254 MagnitudeResponse(n
, h
.data(), hrir
);
257 static bool LoadResponses(MYSOFA_HRTF
*sofaHrtf
, HrirDataT
*hData
, const DelayType delayType
,
260 std::atomic
<uint
> loaded_count
{0u};
262 auto load_proc
= [sofaHrtf
,hData
,delayType
,outRate
,&loaded_count
]() -> bool
264 const uint channels
{(hData
->mChannelType
== CT_STEREO
) ? 2u : 1u};
265 hData
->mHrirsBase
.resize(channels
* size_t{hData
->mIrCount
} * hData
->mIrSize
, 0.0);
266 double *hrirs
= hData
->mHrirsBase
.data();
268 std::vector
<double> restmp
;
269 std::optional
<PPhaseResampler
> resampler
;
270 if(outRate
&& outRate
!= hData
->mIrRate
)
272 resampler
.emplace().init(hData
->mIrRate
, outRate
);
273 restmp
.resize(sofaHrtf
->N
);
276 for(uint si
{0u};si
< sofaHrtf
->M
;++si
)
278 loaded_count
.fetch_add(1u);
281 sofaHrtf
->SourcePosition
.values
[3_uz
*si
],
282 sofaHrtf
->SourcePosition
.values
[3_uz
*si
+ 1],
283 sofaHrtf
->SourcePosition
.values
[3_uz
*si
+ 2]
285 mysofa_c2s(aer
.data());
287 if(std::abs(aer
[1]) >= 89.999f
)
290 aer
[0] = std::fmod(360.0f
- aer
[0], 360.0f
);
292 auto field
= std::find_if(hData
->mFds
.cbegin(), hData
->mFds
.cend(),
293 [&aer
](const HrirFdT
&fld
) -> bool
294 { return (std::abs(aer
[2] - fld
.mDistance
) < 0.001); });
295 if(field
== hData
->mFds
.cend())
298 const double evscale
{180.0 / static_cast<double>(field
->mEvs
.size()-1)};
299 double ef
{(90.0 + aer
[1]) / evscale
};
300 auto ei
= static_cast<uint
>(std::round(ef
));
301 ef
= (ef
- ei
) * evscale
;
302 if(std::abs(ef
) >= 0.1) continue;
304 const double azscale
{360.0 / static_cast<double>(field
->mEvs
[ei
].mAzs
.size())};
305 double af
{aer
[0] / azscale
};
306 auto ai
= static_cast<uint
>(std::round(af
));
307 af
= (af
-ai
) * azscale
;
308 ai
%= static_cast<uint
>(field
->mEvs
[ei
].mAzs
.size());
309 if(std::abs(af
) >= 0.1) continue;
311 HrirAzT
*azd
= &field
->mEvs
[ei
].mAzs
[ai
];
312 if(azd
->mIrs
[0] != nullptr)
314 fprintf(stderr
, "\nMultiple measurements near [ a=%f, e=%f, r=%f ].\n",
315 aer
[0], aer
[1], aer
[2]);
319 for(uint ti
{0u};ti
< channels
;++ti
)
321 azd
->mIrs
[ti
] = &hrirs
[(size_t{hData
->mIrCount
}*ti
+ azd
->mIndex
)*hData
->mIrSize
];
323 std::copy_n(&sofaHrtf
->DataIR
.values
[(size_t{si
}*sofaHrtf
->R
+ ti
)*sofaHrtf
->N
],
324 sofaHrtf
->N
, azd
->mIrs
[ti
]);
327 std::copy_n(&sofaHrtf
->DataIR
.values
[(size_t{si
}*sofaHrtf
->R
+ ti
)*sofaHrtf
->N
],
328 sofaHrtf
->N
, restmp
.data());
329 resampler
->process(sofaHrtf
->N
, restmp
.data(), hData
->mIrSize
, azd
->mIrs
[ti
]);
333 /* Include any per-channel or per-HRIR delays. */
334 if(delayType
== DelayType::I_R
)
336 const float *delayValues
{sofaHrtf
->DataDelay
.values
};
337 for(uint ti
{0u};ti
< channels
;++ti
)
338 azd
->mDelays
[ti
] = delayValues
[ti
] / static_cast<float>(hData
->mIrRate
);
340 else if(delayType
== DelayType::M_R
)
342 const float *delayValues
{sofaHrtf
->DataDelay
.values
};
343 for(uint ti
{0u};ti
< channels
;++ti
)
344 azd
->mDelays
[ti
] = delayValues
[si
*sofaHrtf
->R
+ ti
] /
345 static_cast<float>(hData
->mIrRate
);
349 if(outRate
&& outRate
!= hData
->mIrRate
)
351 const double scale
{static_cast<double>(outRate
) / hData
->mIrRate
};
352 hData
->mIrRate
= outRate
;
353 hData
->mIrPoints
= std::min(static_cast<uint
>(std::ceil(hData
->mIrPoints
*scale
)),
359 std::future_status load_status
{};
360 auto load_future
= std::async(std::launch::async
, load_proc
);
362 load_status
= load_future
.wait_for(std::chrono::milliseconds
{50});
363 printf("\rLoading HRIRs... %u of %u", loaded_count
.load(), sofaHrtf
->M
);
365 } while(load_status
!= std::future_status::ready
);
367 return load_future
.get();
371 /* Calculates the frequency magnitudes of the HRIR set. Work is delegated to
372 * this struct, which runs asynchronously on one or more threads (sharing the
373 * same calculator object).
375 struct MagCalculator
{
376 const uint mFftSize
{};
377 const uint mIrPoints
{};
378 std::vector
<double*> mIrs
{};
379 std::atomic
<size_t> mCurrent
{};
380 std::atomic
<size_t> mDone
{};
384 auto htemp
= std::vector
<complex_d
>(mFftSize
);
388 /* Load the current index to process. */
389 size_t idx
{mCurrent
.load()};
391 /* If the index is at the end, we're done. */
392 if(idx
>= mIrs
.size())
394 /* Otherwise, increment the current index atomically so other
395 * threads know to go to the next one. If this call fails, the
396 * current index was just changed by another thread and the new
397 * value is loaded into idx, which we'll recheck.
399 } while(!mCurrent
.compare_exchange_weak(idx
, idx
+1, std::memory_order_relaxed
));
401 CalcHrirMagnitude(mIrPoints
, mFftSize
, htemp
, mIrs
[idx
]);
403 /* Increment the number of IRs done. */
409 bool LoadSofaFile(const char *filename
, const uint numThreads
, const uint fftSize
,
410 const uint truncSize
, const uint outRate
, const ChannelModeT chanMode
, HrirDataT
*hData
)
413 MySofaHrtfPtr sofaHrtf
{mysofa_load(filename
, &err
)};
416 fprintf(stdout
, "Error: Could not load %s: %s\n", filename
, SofaErrorStr(err
));
420 /* NOTE: Some valid SOFA files are failing this check. */
421 err
= mysofa_check(sofaHrtf
.get());
423 fprintf(stderr
, "Warning: Supposedly malformed source file '%s' (%s).\n", filename
,
426 mysofa_tocartesian(sofaHrtf
.get());
428 /* Make sure emitter and receiver counts are sane. */
431 fprintf(stderr
, "%u emitters not supported\n", sofaHrtf
->E
);
434 if(sofaHrtf
->R
> 2 || sofaHrtf
->R
< 1)
436 fprintf(stderr
, "%u receivers not supported\n", sofaHrtf
->R
);
439 /* Assume R=2 is a stereo measurement, and R=1 is mono left-ear-only. */
440 if(sofaHrtf
->R
== 2 && chanMode
== CM_AllowStereo
)
441 hData
->mChannelType
= CT_STEREO
;
443 hData
->mChannelType
= CT_MONO
;
445 /* Check and set the FFT and IR size. */
446 if(sofaHrtf
->N
> fftSize
)
448 fprintf(stderr
, "Sample points exceeds the FFT size.\n");
451 if(sofaHrtf
->N
< truncSize
)
453 fprintf(stderr
, "Sample points is below the truncation size.\n");
456 hData
->mIrPoints
= sofaHrtf
->N
;
457 hData
->mFftSize
= fftSize
;
458 hData
->mIrSize
= std::max(1u + (fftSize
/2u), sofaHrtf
->N
);
460 /* Assume a default head radius of 9cm. */
461 hData
->mRadius
= 0.09;
463 hData
->mIrRate
= static_cast<uint
>(std::lround(GetSampleRate(sofaHrtf
.get())));
467 DelayType delayType
= PrepareDelay(sofaHrtf
.get());
468 if(delayType
== DelayType::Invalid
)
471 if(!CheckIrData(sofaHrtf
.get()))
473 if(!PrepareLayout(sofaHrtf
->M
, sofaHrtf
->SourcePosition
.values
, hData
))
475 if(!LoadResponses(sofaHrtf
.get(), hData
, delayType
, outRate
))
479 for(uint fi
{0u};fi
< hData
->mFds
.size();fi
++)
482 for(;ei
< hData
->mFds
[fi
].mEvs
.size();ei
++)
485 for(;ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzs
.size();ai
++)
487 HrirAzT
&azd
= hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
488 if(azd
.mIrs
[0] != nullptr) break;
490 if(ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzs
.size())
493 if(ei
>= hData
->mFds
[fi
].mEvs
.size())
495 fprintf(stderr
, "Missing source references [ %d, *, * ].\n", fi
);
498 hData
->mFds
[fi
].mEvStart
= ei
;
499 for(;ei
< hData
->mFds
[fi
].mEvs
.size();ei
++)
501 for(uint ai
{0u};ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzs
.size();ai
++)
503 HrirAzT
&azd
= hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
504 if(azd
.mIrs
[0] == nullptr)
506 fprintf(stderr
, "Missing source reference [ %d, %d, %d ].\n", fi
, ei
, ai
);
514 size_t hrir_total
{0};
515 const uint channels
{(hData
->mChannelType
== CT_STEREO
) ? 2u : 1u};
516 double *hrirs
= hData
->mHrirsBase
.data();
517 for(uint fi
{0u};fi
< hData
->mFds
.size();fi
++)
519 for(uint ei
{0u};ei
< hData
->mFds
[fi
].mEvStart
;ei
++)
521 for(uint ai
{0u};ai
< hData
->mFds
[fi
].mEvs
[ei
].mAzs
.size();ai
++)
523 HrirAzT
&azd
= hData
->mFds
[fi
].mEvs
[ei
].mAzs
[ai
];
524 for(size_t ti
{0u};ti
< channels
;ti
++)
525 azd
.mIrs
[ti
] = &hrirs
[hData
->mIrSize
* (hData
->mIrCount
*ti
+ azd
.mIndex
)];
529 for(uint ei
{hData
->mFds
[fi
].mEvStart
};ei
< hData
->mFds
[fi
].mEvs
.size();ei
++)
530 hrir_total
+= hData
->mFds
[fi
].mEvs
[ei
].mAzs
.size() * channels
;
533 std::atomic
<size_t> hrir_done
{0};
534 auto onset_proc
= [hData
,channels
,&hrir_done
]() -> bool
536 /* Temporary buffer used to calculate the IR's onset. */
537 auto upsampled
= std::vector
<double>(size_t{OnsetRateMultiple
} * hData
->mIrPoints
);
538 /* This resampler is used to help detect the response onset. */
540 rs
.init(hData
->mIrRate
, OnsetRateMultiple
*hData
->mIrRate
);
542 for(auto &field
: hData
->mFds
)
544 for(auto &elev
: field
.mEvs
.subspan(field
.mEvStart
))
546 for(auto &azd
: elev
.mAzs
)
548 for(uint ti
{0};ti
< channels
;ti
++)
550 hrir_done
.fetch_add(1u, std::memory_order_acq_rel
);
551 azd
.mDelays
[ti
] += CalcHrirOnset(rs
, hData
->mIrRate
, hData
->mIrPoints
,
552 upsampled
, azd
.mIrs
[ti
]);
560 std::future_status load_status
{};
561 auto load_future
= std::async(std::launch::async
, onset_proc
);
563 load_status
= load_future
.wait_for(std::chrono::milliseconds
{50});
564 printf("\rCalculating HRIR onsets... %zu of %zu", hrir_done
.load(), hrir_total
);
566 } while(load_status
!= std::future_status::ready
);
568 if(!load_future
.get())
571 MagCalculator calculator
{hData
->mFftSize
, hData
->mIrPoints
};
572 for(auto &field
: hData
->mFds
)
574 for(auto &elev
: field
.mEvs
.subspan(field
.mEvStart
))
576 for(auto &azd
: elev
.mAzs
)
578 for(uint ti
{0};ti
< channels
;ti
++)
579 calculator
.mIrs
.push_back(azd
.mIrs
[ti
]);
584 std::vector
<std::thread
> thrds
;
585 thrds
.reserve(numThreads
);
586 for(size_t i
{0};i
< numThreads
;++i
)
587 thrds
.emplace_back(std::mem_fn(&MagCalculator::Worker
), &calculator
);
590 std::this_thread::sleep_for(std::chrono::milliseconds
{50});
591 count
= calculator
.mDone
.load();
593 printf("\rCalculating HRIR magnitudes... %zu of %zu", count
, calculator
.mIrs
.size());
595 } while(count
!= calculator
.mIrs
.size());
598 for(auto &thrd
: thrds
)