Resample HRIRs after equalization

1 files changed, 93 insertions, 93 deletions

diff --git a/utils/makemhr/makemhr.cpp b/utils/makemhr/makemhr.cpp
index 3ef2a9ef..07ac5f88 100644
--- a/utils/makemhr/makemhr.cpp
+++ b/utils/makemhr/makemhr.cpp
@@ -523,94 +523,6 @@ static int StoreMhr(const HrirDataT *hData, const char *filename)
  *** HRTF processing ***
  ***********************/
 
-// Resamples the HRIRs for use at the given sampling rate.
-static void ResampleHrirs(const uint rate, HrirDataT *hData)
-{
-    struct Resampler {
-        const double scale;
-        const size_t m;
-        double *const resampled;
-
-        /* Resampling from a lower rate to a higher rate. This likely only
-         * works properly when 1 <= scale <= 2.
-         */
-        void upsample(double *ir) const
-        {
-            std::fill_n(resampled, m, 0.0);
-            resampled[0] = ir[0];
-            for(size_t in{1};in < m;++in)
-            {
-                const auto offset = static_cast<double>(in) / scale;
-                const auto out = static_cast<size_t>(offset);
-
-                const double a{offset - static_cast<double>(out)};
-                if(out == m-1)
-                    resampled[out] += ir[in]*(1.0-a);
-                else
-                {
-                    resampled[out  ] += ir[in]*(1.0-a);
-                    resampled[out+1] += ir[in]*a;
-                }
-            }
-            /* This should probably be rescaled according to the scale, however
-             * it'll all be normalized in the end so a constant scalar is fine
-             * to leave.
-             */
-            std::copy_n(resampled, m, ir);
-        }
-
-        /* Resampling from a higher rate to a lower rate. This likely only
-         * works properly when 0.5 <= scale <= 1.0.
-         */
-        void downsample(double *ir) const
-        {
-            std::fill_n(resampled, m, 0.0);
-            resampled[0] = ir[0];
-            for(size_t out{1};out < m;++out)
-            {
-                const auto offset = static_cast<double>(out) * scale;
-                const auto in = static_cast<size_t>(offset);
-
-                const double a{offset - static_cast<double>(in)};
-                if(in == m-1)
-                    resampled[out] = ir[in]*(1.0-a);
-                else
-                    resampled[out] = ir[in]*(1.0-a) + ir[in+1]*a;
-            }
-            std::copy_n(resampled, m, ir);
-        }
-    };
-
-    while(rate > hData->mIrRate*2)
-        ResampleHrirs(hData->mIrRate*2, hData);
-    while(rate < (hData->mIrRate+1)/2)
-        ResampleHrirs((hData->mIrRate+1)/2, hData);
-
-    const auto scale = static_cast<double>(rate) / hData->mIrRate;
-    const size_t m{hData->mFftSize/2u + 1u};
-    auto resampled = std::vector<double>(m);
-
-    const Resampler resampler{scale, m, resampled.data()};
-    auto do_resample = std::bind(
-        std::mem_fn((scale > 1.0) ? &Resampler::upsample : &Resampler::downsample), &resampler,
-        _1);
-
-    const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
-    for(uint fi{0};fi < hData->mFdCount;++fi)
-    {
-        for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvCount;++ei)
-        {
-            for(uint ai{0};ai < hData->mFds[fi].mEvs[ei].mAzCount;++ai)
-            {
-                HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
-                for(uint ti{0};ti < channels;++ti)
-                    do_resample(azd->mIrs[ti]);
-            }
-        }
-    }
-    hData->mIrRate = rate;
-}
-
 /* Balances the maximum HRIR magnitudes of multi-field data sets by
  * independently normalizing each field in relation to the overall maximum.
  * This is done to ignore distance attenuation.
@@ -802,6 +714,94 @@ static void DiffuseFieldEqualize(const uint channels, const uint m, const double
     }
 }
 
+// Resamples the HRIRs for use at the given sampling rate.
+static void ResampleHrirs(const uint rate, HrirDataT *hData)
+{
+    struct Resampler {
+        const double scale;
+        const size_t m;
+
+        /* Resampling from a lower rate to a higher rate. This likely only
+         * works properly when 1 <= scale <= 2.
+         */
+        void upsample(double *resampled, const double *ir) const
+        {
+            std::fill_n(resampled, m, 0.0);
+            resampled[0] = ir[0];
+            for(size_t in{1};in < m;++in)
+            {
+                const auto offset = static_cast<double>(in) / scale;
+                const auto out = static_cast<size_t>(offset);
+
+                const double a{offset - static_cast<double>(out)};
+                if(out == m-1)
+                    resampled[out] += ir[in]*(1.0-a);
+                else
+                {
+                    resampled[out  ] += ir[in]*(1.0-a);
+                    resampled[out+1] += ir[in]*a;
+                }
+            }
+        }
+
+        /* Resampling from a higher rate to a lower rate. This likely only
+         * works properly when 0.5 <= scale <= 1.0.
+         */
+        void downsample(double *resampled, const double *ir) const
+        {
+            resampled[0] = ir[0];
+            for(size_t out{1};out < m;++out)
+            {
+                const auto offset = static_cast<double>(out) * scale;
+                const auto in = static_cast<size_t>(offset);
+
+                const double a{offset - static_cast<double>(in)};
+                if(in == m-1)
+                    resampled[out] = ir[in]*(1.0-a);
+                else
+                    resampled[out] = ir[in]*(1.0-a) + ir[in+1]*a;
+            }
+        }
+    };
+
+    while(rate > hData->mIrRate*2)
+        ResampleHrirs(hData->mIrRate*2, hData);
+    while(rate < (hData->mIrRate+1)/2)
+        ResampleHrirs((hData->mIrRate+1)/2, hData);
+
+    const auto scale = static_cast<double>(rate) / hData->mIrRate;
+    const size_t m{hData->mFftSize/2u + 1u};
+    auto resampled = std::vector<double>(m);
+
+    const Resampler resampler{scale, m};
+    auto do_resample = std::bind(
+        std::mem_fn((scale > 1.0) ? &Resampler::upsample : &Resampler::downsample), &resampler,
+        _1, _2);
+
+    const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
+    for(uint fi{0};fi < hData->mFdCount;++fi)
+    {
+        for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvCount;++ei)
+        {
+            for(uint ai{0};ai < hData->mFds[fi].mEvs[ei].mAzCount;++ai)
+            {
+                HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
+                for(uint ti{0};ti < channels;++ti)
+                {
+                    do_resample(resampled.data(), azd->mIrs[ti]);
+                    /* This should probably be rescaled according to the scale,
+                     * however it'll all be normalized in the end so a constant
+                     * scalar is fine to leave.
+                     */
+                    std::transform(resampled.cbegin(), resampled.cend(), azd->mIrs[ti],
+                        [](const double d) { return std::max(d, EPSILON); });
+                }
+            }
+        }
+    }
+    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).
@@ -1424,11 +1424,6 @@ static int ProcessDefinition(const char *inName, const uint outRate, const Chann
         }
     }
 
-    if(outRate != 0 && outRate != hData.mIrRate)
-    {
-        fprintf(stdout, "Resampling HRIRs...\n");
-        ResampleHrirs(outRate, &hData);
-    }
     if(equalize)
     {
         uint c{(hData.mChannelType == CT_STEREO) ? 2u : 1u};
@@ -1445,6 +1440,11 @@ static int ProcessDefinition(const char *inName, const uint outRate, const Chann
         fprintf(stdout, "Performing diffuse-field equalization...\n");
         DiffuseFieldEqualize(c, m, dfa.data(), &hData);
     }
+    if(outRate != 0 && outRate != hData.mIrRate)
+    {
+        fprintf(stdout, "Resampling HRIRs...\n");
+        ResampleHrirs(outRate, &hData);
+    }
     fprintf(stdout, "Performing minimum phase reconstruction...\n");
     ReconstructHrirs(&hData);
     fprintf(stdout, "Truncating minimum-phase HRIRs...\n");