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.

5 files changed, 92 insertions, 38 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;
diff --git a/utils/makemhr/loaddef.h b/utils/makemhr/loaddef.h
index 34fbb832..63600dcd 100644
--- a/utils/makemhr/loaddef.h
+++ b/utils/makemhr/loaddef.h
@@ -7,7 +7,7 @@
 
 
 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);
 
 #endif /* LOADDEF_H */
diff --git a/utils/makemhr/loadsofa.cpp b/utils/makemhr/loadsofa.cpp
index 7d091be8..18e84637 100644
--- a/utils/makemhr/loadsofa.cpp
+++ b/utils/makemhr/loadsofa.cpp
@@ -37,6 +37,8 @@
 #include <thread>
 #include <vector>
 
+#include "aloptional.h"
+#include "alspan.h"
 #include "makemhr.h"
 #include "polyphase_resampler.h"
 #include "sofa-support.h"
@@ -234,7 +236,7 @@ bool CheckIrData(MYSOFA_HRTF *sofaHrtf)
 /* Calculate the onset time of a HRIR. */
 static constexpr int OnsetRateMultiple{10};
 static double CalcHrirOnset(PPhaseResampler &rs, const uint rate, const uint n,
-    std::vector<double> &upsampled, const double *hrir)
+    al::span<double> upsampled, const double *hrir)
 {
     rs.process(n, hrir, static_cast<uint>(upsampled.size()), upsampled.data());
 
@@ -246,8 +248,7 @@ static double CalcHrirOnset(PPhaseResampler &rs, const uint rate, const uint n,
 }
 
 /* Calculate the magnitude response of a HRIR. */
-static void CalcHrirMagnitude(const uint points, const uint n, std::vector<complex_d> &h,
-    double *hrir)
+static void CalcHrirMagnitude(const uint points, const uint n, al::span<complex_d> h, double *hrir)
 {
     auto iter = std::copy_n(hrir, points, h.begin());
     std::fill(iter, h.end(), complex_d{0.0, 0.0});
@@ -256,16 +257,24 @@ static void CalcHrirMagnitude(const uint points, const uint n, std::vector<compl
     MagnitudeResponse(n, h.data(), hrir);
 }
 
-static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
+static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData, const uint outRate)
 {
     std::atomic<uint> loaded_count{0u};
 
-    auto load_proc = [sofaHrtf,hData,&loaded_count]() -> bool
+    auto load_proc = [sofaHrtf,hData,outRate,&loaded_count]() -> bool
     {
         const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
         hData->mHrirsBase.resize(channels * hData->mIrCount * hData->mIrSize, 0.0);
         double *hrirs = hData->mHrirsBase.data();
 
+        std::unique_ptr<double[]> restmp;
+        al::optional<PPhaseResampler> resampler;
+        if(outRate && outRate != hData->mIrRate)
+        {
+            resampler.emplace().init(hData->mIrRate, outRate);
+            restmp = std::make_unique<double[]>(sofaHrtf->N);
+        }
+
         for(uint si{0u};si < sofaHrtf->M;++si)
         {
             loaded_count.fetch_add(1u);
@@ -313,8 +322,15 @@ static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
             for(uint ti{0u};ti < channels;++ti)
             {
                 azd->mIrs[ti] = &hrirs[hData->mIrSize * (hData->mIrCount*ti + azd->mIndex)];
-                std::copy_n(&sofaHrtf->DataIR.values[(si*sofaHrtf->R + ti)*sofaHrtf->N],
-                    hData->mIrPoints, azd->mIrs[ti]);
+                if(!resampler)
+                    std::copy_n(&sofaHrtf->DataIR.values[(si*sofaHrtf->R + ti)*sofaHrtf->N],
+                        sofaHrtf->N, azd->mIrs[ti]);
+                else
+                {
+                    std::copy_n(&sofaHrtf->DataIR.values[(si*sofaHrtf->R + ti)*sofaHrtf->N],
+                        sofaHrtf->N, restmp.get());
+                    resampler->process(sofaHrtf->N, restmp.get(), hData->mIrSize, azd->mIrs[ti]);
+                }
             }
 
             /* TODO: Since some SOFA files contain minimum phase HRIRs,
@@ -322,6 +338,14 @@ static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
              * (when available) to reconstruct the HRTDs.
              */
         }
+
+        if(outRate && outRate != hData->mIrRate)
+        {
+            const double scale{static_cast<double>(outRate) / hData->mIrRate};
+            hData->mIrRate = outRate;
+            hData->mIrPoints = std::min(static_cast<uint>(std::ceil(hData->mIrPoints*scale)),
+                hData->mIrSize);
+        }
         return true;
     };
 
@@ -376,7 +400,7 @@ struct MagCalculator {
 };
 
 bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSize,
-    const uint truncSize, const ChannelModeT chanMode, HrirDataT *hData)
+    const uint truncSize, const uint outRate, const ChannelModeT chanMode, HrirDataT *hData)
 {
     int err;
     MySofaHrtfPtr sofaHrtf{mysofa_load(filename, &err)};
@@ -435,7 +459,7 @@ bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSiz
     if(!PrepareLayout(sofaHrtf->M, sofaHrtf->SourcePosition.values, hData))
         return false;
 
-    if(!LoadResponses(sofaHrtf.get(), hData))
+    if(!LoadResponses(sofaHrtf.get(), hData, outRate))
         return false;
     sofaHrtf = nullptr;
 
diff --git a/utils/makemhr/loadsofa.h b/utils/makemhr/loadsofa.h
index 803bcf88..82dce85a 100644
--- a/utils/makemhr/loadsofa.h
+++ b/utils/makemhr/loadsofa.h
@@ -5,6 +5,6 @@
 
 
 bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSize,
-    const uint truncSize, const ChannelModeT chanMode, HrirDataT *hData);
+    const uint truncSize, const uint outRate, const ChannelModeT chanMode, HrirDataT *hData);
 
 #endif /* LOADSOFA_H */
diff --git a/utils/makemhr/makemhr.cpp b/utils/makemhr/makemhr.cpp
index 554bdfe0..8f7ae630 100644
--- a/utils/makemhr/makemhr.cpp
+++ b/utils/makemhr/makemhr.cpp
@@ -1392,7 +1392,7 @@ static int ProcessDefinition(const char *inName, const uint outRate, const Chann
     {
         inName = "stdin";
         fprintf(stdout, "Reading HRIR definition from %s...\n", inName);
-        if(!LoadDefInput(std::cin, nullptr, 0, inName, fftSize, truncSize, chanMode, &hData))
+        if(!LoadDefInput(std::cin, nullptr, 0, inName, fftSize, truncSize, outRate, chanMode, &hData))
             return 0;
     }
     else
@@ -1418,13 +1418,13 @@ static int ProcessDefinition(const char *inName, const uint outRate, const Chann
         {
             input = nullptr;
             fprintf(stdout, "Reading HRTF data from %s...\n", inName);
-            if(!LoadSofaFile(inName, numThreads, fftSize, truncSize, chanMode, &hData))
+            if(!LoadSofaFile(inName, numThreads, fftSize, truncSize, outRate, chanMode, &hData))
                 return 0;
         }
         else
         {
             fprintf(stdout, "Reading HRIR definition from %s...\n", inName);
-            if(!LoadDefInput(*input, startbytes, startbytecount, inName, fftSize, truncSize, chanMode, &hData))
+            if(!LoadDefInput(*input, startbytes, startbytecount, inName, fftSize, truncSize, outRate, chanMode, &hData))
                 return 0;
         }
     }