Resample before frequency analysis

We want to resample before minimum phase reconstruction since that changes the
phase relationship of the sampled signal, introducing a slight bit of noise
from truncated sampling. It's not clear that the frequency domain resampling
method is accurate, so resampling prior to frequency analysis is an alternative
to ensure the resulting frequencies are given the proper phase for sampling.

This also cleans up some micro allocations in loops.

1 files changed, 53 insertions, 23 deletions

diff --git a/utils/makemhr/loaddef.cpp b/utils/makemhr/loaddef.cpp
index d325eda1..ab505f47 100644
--- a/utils/makemhr/loaddef.cpp
+++ b/utils/makemhr/loaddef.cpp
@@ -37,8 +37,11 @@
 #include <vector>
 
 #include "alfstream.h"
+#include "aloptional.h"
+#include "alspan.h"
 #include "alstring.h"
 #include "makemhr.h"
+#include "polyphase_resampler.h"
 
 #include "mysofa.h"
 
@@ -1707,14 +1710,11 @@ static int MatchTargetEar(const char *ident)
 
 // Calculate the onset time of an HRIR and average it with any existing
 // timing for its field, elevation, azimuth, and ear.
-static double AverageHrirOnset(const uint rate, const uint n, const double *hrir, const double f, const double onset)
+static constexpr int OnsetRateMultiple{10};
+static double AverageHrirOnset(PPhaseResampler &rs, al::span<double> upsampled, const uint rate,
+    const uint n, const double *hrir, const double f, const double onset)
 {
-    std::vector<double> upsampled(10 * n);
-    {
-        PPhaseResampler rs;
-        rs.init(rate, 10 * rate);
-        rs.process(n, hrir, 10 * n, upsampled.data());
-    }
+    rs.process(n, hrir, static_cast<uint>(upsampled.size()), upsampled.data());
 
     auto abs_lt = [](const double &lhs, const double &rhs) -> bool
     { return std::abs(lhs) < std::abs(rhs); };
@@ -1731,9 +1731,9 @@ static void AverageHrirMagnitude(const uint points, const uint n, const double *
     std::vector<double> r(n);
 
     for(i = 0;i < points;i++)
-        h[i] = complex_d{hrir[i], 0.0};
+        h[i] = hrir[i];
     for(;i < n;i++)
-        h[i] = complex_d{0.0, 0.0};
+        h[i] = 0.0;
     FftForward(n, h.data());
     MagnitudeResponse(n, h.data(), r.data());
     for(i = 0;i < m;i++)
@@ -1741,15 +1741,27 @@ static void AverageHrirMagnitude(const uint points, const uint n, const double *
 }
 
 // Process the list of sources in the data set definition.
-static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
+static int ProcessSources(TokenReaderT *tr, HrirDataT *hData, const uint outRate)
 {
     uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
     hData->mHrirsBase.resize(channels * hData->mIrCount * hData->mIrSize);
     double *hrirs = hData->mHrirsBase.data();
-    std::vector<double> hrir(hData->mIrPoints);
+    auto hrir = std::make_unique<double[]>(hData->mIrSize);
     uint line, col, fi, ei, ai;
     int count;
 
+    std::vector<double> onsetSamples(OnsetRateMultiple * hData->mIrPoints);
+    PPhaseResampler onsetResampler;
+    onsetResampler.init(hData->mIrRate, OnsetRateMultiple*hData->mIrRate);
+
+    al::optional<PPhaseResampler> resampler;
+    if(outRate && outRate != hData->mIrRate)
+        resampler.emplace().init(hData->mIrRate, outRate);
+    const double rateScale{outRate ? static_cast<double>(outRate) / hData->mIrRate : 1.0};
+    const uint irPoints{outRate
+        ? std::min(static_cast<uint>(std::ceil(hData->mIrPoints*rateScale)), hData->mIrPoints)
+        : hData->mIrPoints};
+
     printf("Loading sources...");
     fflush(stdout);
     count = 0;
@@ -1862,17 +1874,24 @@ static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
                     return 0;
                 }
 
-                ExtractSofaHrir(sofa, si, 0, src.mOffset, hData->mIrPoints, hrir.data());
+                ExtractSofaHrir(sofa, si, 0, src.mOffset, hData->mIrPoints, hrir.get());
                 azd->mIrs[0] = &hrirs[hData->mIrSize * azd->mIndex];
-                azd->mDelays[0] = AverageHrirOnset(hData->mIrRate, hData->mIrPoints, hrir.data(), 1.0, azd->mDelays[0]);
-                AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir.data(), 1.0, azd->mIrs[0]);
+                azd->mDelays[0] = AverageHrirOnset(onsetResampler, onsetSamples, hData->mIrRate,
+                    hData->mIrPoints, hrir.get(), 1.0, azd->mDelays[0]);
+                if(resampler)
+                    resampler->process(hData->mIrPoints, hrir.get(), hData->mIrSize, hrir.get());
+                AverageHrirMagnitude(irPoints, hData->mFftSize, hrir.get(), 1.0, azd->mIrs[0]);
 
                 if(src.mChannel == 1)
                 {
-                    ExtractSofaHrir(sofa, si, 1, src.mOffset, hData->mIrPoints, hrir.data());
+                    ExtractSofaHrir(sofa, si, 1, src.mOffset, hData->mIrPoints, hrir.get());
                     azd->mIrs[1] = &hrirs[hData->mIrSize * (hData->mIrCount + azd->mIndex)];
-                    azd->mDelays[1] = AverageHrirOnset(hData->mIrRate, hData->mIrPoints, hrir.data(), 1.0, azd->mDelays[1]);
-                    AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir.data(), 1.0, azd->mIrs[1]);
+                    azd->mDelays[1] = AverageHrirOnset(onsetResampler, onsetSamples,
+                        hData->mIrRate, hData->mIrPoints, hrir.get(), 1.0, azd->mDelays[1]);
+                    if(resampler)
+                        resampler->process(hData->mIrPoints, hrir.get(), hData->mIrSize,
+                            hrir.get());
+                    AverageHrirMagnitude(irPoints, hData->mFftSize, hrir.get(), 1.0, azd->mIrs[1]);
                 }
 
                 // TODO: Since some SOFA files contain minimum phase HRIRs,
@@ -1911,7 +1930,7 @@ static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
             printf("\rLoading sources... %d file%s", count, (count==1)?"":"s");
             fflush(stdout);
 
-            if(!LoadSource(&src, hData->mIrRate, hData->mIrPoints, hrir.data()))
+            if(!LoadSource(&src, hData->mIrRate, hData->mIrPoints, hrir.get()))
                 return 0;
 
             uint ti{0};
@@ -1929,8 +1948,12 @@ static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
                 }
             }
             azd->mIrs[ti] = &hrirs[hData->mIrSize * (ti * hData->mIrCount + azd->mIndex)];
-            azd->mDelays[ti] = AverageHrirOnset(hData->mIrRate, hData->mIrPoints, hrir.data(), 1.0 / factor[ti], azd->mDelays[ti]);
-            AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir.data(), 1.0 / factor[ti], azd->mIrs[ti]);
+            azd->mDelays[ti] = AverageHrirOnset(onsetResampler, onsetSamples, hData->mIrRate,
+                hData->mIrPoints, hrir.get(), 1.0 / factor[ti], azd->mDelays[ti]);
+            if(resampler)
+                resampler->process(hData->mIrPoints, hrir.get(), hData->mIrSize, hrir.get());
+            AverageHrirMagnitude(irPoints, hData->mFftSize, hrir.get(), 1.0 / factor[ti],
+                azd->mIrs[ti]);
             factor[ti] += 1.0;
             if(!TrIsOperator(tr, "+"))
                 break;
@@ -1951,6 +1974,13 @@ static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
         }
     }
     printf("\n");
+    hrir = nullptr;
+    if(resampler)
+    {
+        hData->mIrRate = outRate;
+        hData->mIrPoints = irPoints;
+        resampler.reset();
+    }
     for(fi = 0;fi < hData->mFdCount;fi++)
     {
         for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
@@ -2012,14 +2042,14 @@ static int ProcessSources(TokenReaderT *tr, HrirDataT *hData)
 
 
 bool LoadDefInput(std::istream &istream, const char *startbytes, std::streamsize startbytecount,
-    const char *filename, const uint fftSize, const uint truncSize, const ChannelModeT chanMode,
-    HrirDataT *hData)
+    const char *filename, const uint fftSize, const uint truncSize, const uint outRate,
+    const ChannelModeT chanMode, HrirDataT *hData)
 {
     TokenReaderT tr{istream};
 
     TrSetup(startbytes, startbytecount, filename, &tr);
     if(!ProcessMetrics(&tr, fftSize, truncSize, chanMode, hData)
-        || !ProcessSources(&tr, hData))
+        || !ProcessSources(&tr, hData, outRate))
         return false;
 
     return true;