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.

common/polyphase_resampler.h

@@ -4,6 +4,8 @@
 #include <vector>
 
 
+using uint = unsigned int;
+
 /* This is a polyphase sinc-filtered resampler. It is built for very high
  * quality results, rather than real-time performance.
  *
@@ -32,8 +34,6 @@
  */
 
 struct PPhaseResampler {
-    using uint = unsigned int;
-
     void init(const uint srcRate, const uint dstRate);
     void process(const uint inN, const double *in, const uint outN, double *out);
 

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);
+    rs.process(n, hrir, static_cast<uint>(upsampled.size()), upsampled.data());
-    {
-        PPhaseResampler rs;
-        rs.init(rate, 10 * rate);
-        rs.process(n, hrir, 10 * n, 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]);
+                azd->mDelays[0] = AverageHrirOnset(onsetResampler, onsetSamples, hData->mIrRate,
-                AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir.data(), 1.0, azd->mIrs[0]);
+                    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]);
+                    azd->mDelays[1] = AverageHrirOnset(onsetResampler, onsetSamples,
-                    AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir.data(), 1.0, azd->mIrs[1]);
+                        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]);
+            azd->mDelays[ti] = AverageHrirOnset(onsetResampler, onsetSamples, hData->mIrRate,
-            AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir.data(), 1.0 / factor[ti], azd->mIrs[ti]);
+                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,
+    const char *filename, const uint fftSize, const uint truncSize, const uint outRate,
-    HrirDataT *hData)
+    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;

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,
+    const char *filename, const uint fftSize, const uint truncSize, const uint outRate,
-    HrirDataT *hData);
+    const ChannelModeT chanMode, HrirDataT *hData);
 
 #endif /* LOADDEF_H */

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,
+static void CalcHrirMagnitude(const uint points, const uint n, al::span<complex_d> h, double *hrir)
-    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],
+                if(!resampler)
-                    hData->mIrPoints, azd->mIrs[ti]);
+                    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;
 

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 */

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;
         }
     }