diff --git a/common/polyphase_resampler.h b/common/polyphase_resampler.h index e9833d12..0b343d5a 100644 --- a/common/polyphase_resampler.h +++ b/common/polyphase_resampler.h @@ -4,6 +4,8 @@ #include +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); 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 #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 upsampled, const uint rate, + const uint n, const double *hrir, const double f, const double onset) { - std::vector upsampled(10 * n); - { - PPhaseResampler rs; - rs.init(rate, 10 * rate); - rs.process(n, hrir, 10 * n, upsampled.data()); - } + rs.process(n, hrir, static_cast(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 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 hrir(hData->mIrPoints); + auto hrir = std::make_unique(hData->mIrSize); uint line, col, fi, ei, ai; int count; + std::vector onsetSamples(OnsetRateMultiple * hData->mIrPoints); + PPhaseResampler onsetResampler; + onsetResampler.init(hData->mIrRate, OnsetRateMultiple*hData->mIrRate); + + al::optional resampler; + if(outRate && outRate != hData->mIrRate) + resampler.emplace().init(hData->mIrRate, outRate); + const double rateScale{outRate ? static_cast(outRate) / hData->mIrRate : 1.0}; + const uint irPoints{outRate + ? std::min(static_cast(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 #include +#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 &upsampled, const double *hrir) + al::span upsampled, const double *hrir) { rs.process(n, hrir, static_cast(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 &h, - double *hrir) +static void CalcHrirMagnitude(const uint points, const uint n, al::span 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 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 restmp; + al::optional resampler; + if(outRate && outRate != hData->mIrRate) + { + resampler.emplace().init(hData->mIrRate, outRate); + restmp = std::make_unique(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(outRate) / hData->mIrRate}; + hData->mIrRate = outRate; + hData->mIrPoints = std::min(static_cast(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; } }