Resample HRIRs after equalization

This commit is contained in:
Chris Robinson 2020-06-17 13:03:26 -07:00
parent 4cfd63a1c7
commit beba71a6f9

View File

@ -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");