From 24393ab192a15cbbcd12dab7d9a6268d59bc9fae Mon Sep 17 00:00:00 2001
From: Chris Robinson <chris.kcat@gmail.com>
Date: Sun, 21 Jun 2020 00:29:57 -0700
Subject: Synthesize missing elevations in the frequency domain
This should help avoid destructive phase interference. The occlusion low-pass
filter is still applied in the time domain due to no clear topology (cutoff
frequency, slope, bandwidth, etc).
---
utils/makemhr/makemhr.cpp | 288 +++++++++++++++++++++++-----------------------
1 file changed, 147 insertions(+), 141 deletions(-)
(limited to 'utils/makemhr')
diff --git a/utils/makemhr/makemhr.cpp b/utils/makemhr/makemhr.cpp
index e9e950fb..c209c5c1 100644
--- a/utils/makemhr/makemhr.cpp
+++ b/utils/makemhr/makemhr.cpp
@@ -802,115 +802,6 @@ static void ResampleHrirs(const uint rate, HrirDataT *hData)
hData->mIrRate = rate;
}
-/* Perform minimum-phase reconstruction using the magnitude responses of the
- * HRIR set. Work is delegated to this struct, which runs asynchronously on one
- * or more threads (sharing the same reconstructor object).
- */
-struct HrirReconstructor {
- std::vector<double*> mIrs;
- std::atomic<size_t> mCurrent;
- std::atomic<size_t> mDone;
- uint mFftSize;
- uint mIrPoints;
-
- void Worker()
- {
- auto h = std::vector<complex_d>(mFftSize);
-
- while(1)
- {
- /* Load the current index to process. */
- size_t idx{mCurrent.load()};
- do {
- /* If the index is at the end, we're done. */
- if(idx >= mIrs.size())
- return;
- /* Otherwise, increment the current index atomically so other
- * threads know to go to the next one. If this call fails, the
- * current index was just changed by another thread and the new
- * value is loaded into idx, which we'll recheck.
- */
- } while(!mCurrent.compare_exchange_weak(idx, idx+1, std::memory_order_relaxed));
-
- /* Now do the reconstruction, and apply the inverse FFT to get the
- * time-domain response.
- */
- MinimumPhase(mFftSize, mIrs[idx], h.data());
- FftInverse(mFftSize, h.data());
- for(uint i{0u};i < mIrPoints;++i)
- mIrs[idx][i] = h[i].real();
-
- /* Increment the number of IRs done. */
- mDone.fetch_add(1);
- }
- }
-};
-
-static void ReconstructHrirs(const HrirDataT *hData, const uint numThreads)
-{
- const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
-
- /* Count the number of IRs to process (excluding elevations that will be
- * synthesized later).
- */
- size_t total{0};
- for(uint fi{0u};fi < hData->mFdCount;fi++)
- {
- for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvCount;ei++)
- total += hData->mFds[fi].mEvs[ei].mAzCount;
- }
- total *= channels;
-
- /* Set up the reconstructor with the needed size info and pointers to the
- * IRs to process.
- */
- HrirReconstructor reconstructor;
- reconstructor.mIrs.reserve(total);
- reconstructor.mCurrent.store(0, std::memory_order_relaxed);
- reconstructor.mDone.store(0, std::memory_order_relaxed);
- reconstructor.mFftSize = hData->mFftSize;
- reconstructor.mIrPoints = hData->mIrPoints;
- for(uint fi{0u};fi < hData->mFdCount;fi++)
- {
- const HrirFdT &field = hData->mFds[fi];
- for(uint ei{field.mEvStart};ei < field.mEvCount;ei++)
- {
- const HrirEvT &elev = field.mEvs[ei];
- for(uint ai{0u};ai < elev.mAzCount;ai++)
- {
- const HrirAzT &azd = elev.mAzs[ai];
- for(uint ti{0u};ti < channels;ti++)
- reconstructor.mIrs.push_back(azd.mIrs[ti]);
- }
- }
- }
-
- /* Launch threads to work on reconstruction. */
- std::vector<std::thread> thrds;
- thrds.reserve(numThreads);
- for(size_t i{0};i < numThreads;++i)
- thrds.emplace_back(std::mem_fn(&HrirReconstructor::Worker), &reconstructor);
-
- /* Keep track of the number of IRs done, periodically reporting it. */
- size_t count;
- do {
- std::this_thread::sleep_for(std::chrono::milliseconds{50});
-
- count = reconstructor.mDone.load();
- size_t pcdone{count * 100 / total};
-
- printf("\r%3zu%% done (%zu of %zu)", pcdone, count, total);
- fflush(stdout);
- } while(count != total);
- fputc('\n', stdout);
-
- for(auto &thrd : thrds)
- {
- if(thrd.joinable())
- thrd.join();
- }
-}
-
/* Given field and elevation indices and an azimuth, calculate the indices of
* the two HRIRs that bound the coordinate along with a factor for
* calculating the continuous HRIR using interpolation.
@@ -1050,17 +941,15 @@ static void SynthesizeOnsets(HrirDataT *hData)
}
/* Attempt to synthesize any missing HRIRs at the bottom elevations of each
- * field. Right now this just blends the lowest elevation HRIRs together and
- * applies some attenuation and high frequency damping. It is a simple, if
- * inaccurate model.
+ * field. Right now this just blends the lowest elevation HRIRs together. It
+ * is a simple, if inaccurate model.
*/
static void SynthesizeHrirs(HrirDataT *hData)
{
const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
- const uint irSize{hData->mIrPoints};
- const double beta{3.5e-6 * hData->mIrRate};
+ const uint m{hData->mFftSize/2u + 1u};
- auto proc_field = [channels,irSize,beta](HrirFdT &field) -> void
+ auto proc_field = [channels,m](HrirFdT &field) -> void
{
const uint oi{field.mEvStart};
if(oi <= 0) return;
@@ -1074,11 +963,10 @@ static void SynthesizeHrirs(HrirDataT *hData)
* lowest immediate-right response for the right ear. Given no comb
* effects as a result of the left response reaching the right ear
* and vice-versa, this produces a decent phantom-center response
- * underneath the head (with a low-pass filter applied to simulate
- * body occlusion).
+ * underneath the head.
*/
CalcAzIndices(field, oi, ((ti==0) ? -M_PI : M_PI) / 2.0, &a0, &a1, &af);
- for(uint i{0u};i < irSize;i++)
+ for(uint i{0u};i < m;i++)
{
field.mEvs[0].mAzs[0].mIrs[ti][i] = Lerp(field.mEvs[oi].mAzs[a0].mIrs[ti][i],
field.mEvs[oi].mAzs[a1].mIrs[ti][i], af);
@@ -1088,7 +976,6 @@ static void SynthesizeHrirs(HrirDataT *hData)
for(uint ei{1u};ei < field.mEvStart;ei++)
{
const double of{static_cast<double>(ei) / field.mEvStart};
- const double b{(1.0 - of) * beta};
for(uint ai{0u};ai < field.mEvs[ei].mAzCount;ai++)
{
uint a0, a1;
@@ -1097,15 +984,147 @@ static void SynthesizeHrirs(HrirDataT *hData)
CalcAzIndices(field, oi, field.mEvs[ei].mAzs[ai].mAzimuth, &a0, &a1, &af);
for(uint ti{0u};ti < channels;ti++)
{
- double lp[4]{};
- for(uint i{0u};i < irSize;i++)
+ for(uint i{0u};i < m;i++)
{
/* Blend the two defined HRIRs closest to this azimuth,
* then blend that with the synthesized -90 elevation.
*/
const double s1{Lerp(field.mEvs[oi].mAzs[a0].mIrs[ti][i],
field.mEvs[oi].mAzs[a1].mIrs[ti][i], af)};
- const double s0{Lerp(field.mEvs[0].mAzs[0].mIrs[ti][i], s1, of)};
+ field.mEvs[ei].mAzs[ai].mIrs[ti][i] = Lerp(
+ field.mEvs[0].mAzs[0].mIrs[ti][i], s1, of);
+ }
+ }
+ }
+ }
+ };
+ std::for_each(hData->mFds.begin(), hData->mFds.begin()+hData->mFdCount, proc_field);
+}
+
+/* Perform minimum-phase reconstruction using the magnitude responses of the
+ * HRIR set. Work is delegated to this struct, which runs asynchronously on one
+ * or more threads (sharing the same reconstructor object).
+ */
+struct HrirReconstructor {
+ std::vector<double*> mIrs;
+ std::atomic<size_t> mCurrent;
+ std::atomic<size_t> mDone;
+ uint mFftSize;
+ uint mIrPoints;
+
+ void Worker()
+ {
+ auto h = std::vector<complex_d>(mFftSize);
+
+ while(1)
+ {
+ /* Load the current index to process. */
+ size_t idx{mCurrent.load()};
+ do {
+ /* If the index is at the end, we're done. */
+ if(idx >= mIrs.size())
+ return;
+ /* Otherwise, increment the current index atomically so other
+ * threads know to go to the next one. If this call fails, the
+ * current index was just changed by another thread and the new
+ * value is loaded into idx, which we'll recheck.
+ */
+ } while(!mCurrent.compare_exchange_weak(idx, idx+1, std::memory_order_relaxed));
+
+ /* Now do the reconstruction, and apply the inverse FFT to get the
+ * time-domain response.
+ */
+ MinimumPhase(mFftSize, mIrs[idx], h.data());
+ FftInverse(mFftSize, h.data());
+ for(uint i{0u};i < mIrPoints;++i)
+ mIrs[idx][i] = h[i].real();
+
+ /* Increment the number of IRs done. */
+ mDone.fetch_add(1);
+ }
+ }
+};
+
+static void ReconstructHrirs(const HrirDataT *hData, const uint numThreads)
+{
+ const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
+
+ /* Set up the reconstructor with the needed size info and pointers to the
+ * IRs to process.
+ */
+ HrirReconstructor reconstructor;
+ reconstructor.mCurrent.store(0, std::memory_order_relaxed);
+ reconstructor.mDone.store(0, std::memory_order_relaxed);
+ reconstructor.mFftSize = hData->mFftSize;
+ reconstructor.mIrPoints = hData->mIrPoints;
+ for(uint fi{0u};fi < hData->mFdCount;fi++)
+ {
+ const HrirFdT &field = hData->mFds[fi];
+ for(uint ei{0};ei < field.mEvCount;ei++)
+ {
+ const HrirEvT &elev = field.mEvs[ei];
+ for(uint ai{0u};ai < elev.mAzCount;ai++)
+ {
+ const HrirAzT &azd = elev.mAzs[ai];
+ for(uint ti{0u};ti < channels;ti++)
+ reconstructor.mIrs.push_back(azd.mIrs[ti]);
+ }
+ }
+ }
+
+ /* Launch threads to work on reconstruction. */
+ std::vector<std::thread> thrds;
+ thrds.reserve(numThreads);
+ for(size_t i{0};i < numThreads;++i)
+ thrds.emplace_back(std::mem_fn(&HrirReconstructor::Worker), &reconstructor);
+
+ /* Keep track of the number of IRs done, periodically reporting it. */
+ size_t count;
+ do {
+ std::this_thread::sleep_for(std::chrono::milliseconds{50});
+
+ count = reconstructor.mDone.load();
+ size_t pcdone{count * 100 / reconstructor.mIrs.size()};
+
+ printf("\r%3zu%% done (%zu of %zu)", pcdone, count, reconstructor.mIrs.size());
+ fflush(stdout);
+ } while(count < reconstructor.mIrs.size());
+ fputc('\n', stdout);
+
+ for(auto &thrd : thrds)
+ {
+ if(thrd.joinable())
+ thrd.join();
+ }
+}
+
+/* Attempt to occlude the synthesized HRIRs. Right now this just applies some
+ * attenuation and high frequency damping. It is a simple, if inaccurate
+ * model.
+ */
+static void OccludeHrirs(HrirDataT *hData)
+{
+ const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
+ const uint irSize{hData->mIrPoints};
+ const double beta{3.5e-6 * hData->mIrRate};
+
+ auto proc_field = [channels,irSize,beta](HrirFdT &field) -> void
+ {
+ const uint oi{field.mEvStart};
+ if(oi <= 0) return;
+
+ for(uint ei{0u};ei < field.mEvStart;ei++)
+ {
+ const double of{static_cast<double>(ei) / field.mEvStart};
+ const double b{(1.0 - of) * beta};
+ for(uint ai{0u};ai < field.mEvs[ei].mAzCount;ai++)
+ {
+ for(uint ti{0u};ti < channels;ti++)
+ {
+ double lp[4]{};
+ for(uint i{0u};i < irSize;i++)
+ {
+ const double s0{field.mEvs[ei].mAzs[ai].mIrs[ti][i]};
/* Apply a low-pass to simulate body occlusion. */
lp[0] = Lerp(s0, lp[0], b);
lp[1] = Lerp(lp[0], lp[1], b);
@@ -1116,21 +1135,6 @@ static void SynthesizeHrirs(HrirDataT *hData)
}
}
}
- for(uint ti{0u};ti < channels;ti++)
- {
- const double b{beta};
- double lp[4]{};
- for(uint i{0u};i < irSize;i++)
- {
- const double s0{field.mEvs[0].mAzs[0].mIrs[ti][i]};
- lp[0] = Lerp(s0, lp[0], b);
- lp[1] = Lerp(lp[0], lp[1], b);
- lp[2] = Lerp(lp[1], lp[2], b);
- lp[3] = Lerp(lp[2], lp[3], b);
- field.mEvs[0].mAzs[0].mIrs[ti][i] = lp[3];
- }
- }
- field.mEvStart = 0;
};
std::for_each(hData->mFds.begin(), hData->mFds.begin()+hData->mFdCount, proc_field);
}
@@ -1460,14 +1464,16 @@ static int ProcessDefinition(const char *inName, const uint outRate, const Chann
fprintf(stdout, "Resampling HRIRs...\n");
ResampleHrirs(outRate, &hData);
}
- fprintf(stdout, "Performing minimum phase reconstruction...\n");
- ReconstructHrirs(&hData, numThreads);
- fprintf(stdout, "Truncating minimum-phase HRIRs...\n");
- hData.mIrPoints = truncSize;
fprintf(stdout, "Synthesizing missing elevations...\n");
if(model == HM_DATASET)
SynthesizeOnsets(&hData);
SynthesizeHrirs(&hData);
+ fprintf(stdout, "Performing minimum phase reconstruction...\n");
+ ReconstructHrirs(&hData, numThreads);
+ fprintf(stdout, "Truncating minimum-phase HRIRs...\n");
+ hData.mIrPoints = truncSize;
+ fprintf(stdout, "Occluding synthesized HRIRs...\n");
+ OccludeHrirs(&hData);
fprintf(stdout, "Normalizing final HRIRs...\n");
NormalizeHrirs(&hData);
fprintf(stdout, "Calculating impulse delays...\n");
--
cgit v1.2.3