Dima353/openal-soft
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Resample HRIRs prior to minimum phase reconstruction

Browse Source
This commit is contained in:
Chris Robinson
2017-06-01 04:06:08 -07:00
parent f54946f9cb
commit 1c04cbad04
1 changed files with 41 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files
utils/makehrtf.c
+41 -36
View File
@@ -1267,9 +1267,9 @@ static void ResamplerSetup(ResamplerT *rs, const uint srcRate, const uint dstRat
rs->mP = dstRate / gcd;
rs->mQ = srcRate / gcd;
/* The cutoff is adjusted by half the transition width, so the transition
* ends before the nyquist (0.5). Both are scaled by the downsampling
* factor.
*/
* ends before the nyquist (0.5). Both are scaled by the downsampling
* factor.
*/
if(rs->mP > rs->mQ)
{
cutoff = 0.45 / rs->mP;
@@ -1923,13 +1923,12 @@ static int StoreMhr(const HrirDataT *hData, const char *filename)
// Calculate the onset time of an HRIR and average it with any existing
// timing for its elevation and azimuth.
static void AverageHrirOnset(const double *hrir, const double f, const uint ei, const uint ai, const HrirDataT *hData)
static void AverageHrirOnset(const double *hrir, const uint n, const double f, const uint ei, const uint ai, const HrirDataT *hData)
{
double mag;
uint n, i, j;
uint i, j;
mag = 0.0;
n = hData->mIrPoints;
for(i = 0;i < n;i++)
mag = fmax(fabs(hrir[i]), mag);
mag *= 0.15;
@@ -1944,7 +1943,7 @@ static void AverageHrirOnset(const double *hrir, const double f, const uint ei,
// Calculate the magnitude response of an HRIR and average it with any
// existing responses for its elevation and azimuth.
static void AverageHrirMagnitude(const double *hrir, const double f, const uint ei, const uint ai, const HrirDataT *hData)
static void AverageHrirMagnitude(const double *hrir, const uint npoints, const double f, const uint ei, const uint ai, const HrirDataT *hData)
{
double *re, *im;
uint n, m, i, j;
@@ -1952,7 +1951,7 @@ static void AverageHrirMagnitude(const double *hrir, const double f, const uint
n = hData->mFftSize;
re = CreateArray(n);
im = CreateArray(n);
for(i = 0;i < hData->mIrPoints;i++)
for(i = 0;i < npoints;i++)
{
re[i] = hrir[i];
im[i] = 0.0;
@@ -2108,23 +2107,6 @@ static void ReconstructHrirs(const HrirDataT *hData)
DestroyArray (re);
}
// Resamples the HRIRs for use at the given sampling rate.
static void ResampleHrirs(const uint rate, HrirDataT *hData)
{
uint n, step, start, end, j;
ResamplerT rs;
ResamplerSetup(&rs, hData->mIrRate, rate);
n = hData->mIrPoints;
step = hData->mIrSize;
start = hData->mEvOffset[hData->mEvStart] * step;
end = hData->mIrCount * step;
for(j = start;j < end;j += step)
ResamplerRun(&rs, n, &hData->mHrirs[j], n, &hData->mHrirs[j]);
ResamplerClear(&rs);
hData->mIrRate = rate;
}
/* Given an elevation index and an azimuth, calculate the indices of the two
* HRIRs that bound the coordinate along with a factor for calculating the
* continous HRIR using interpolation.
@@ -2614,17 +2596,37 @@ static int ReadSourceRef(TokenReaderT *tr, SourceRefT *src)
}
// Process the list of sources in the data set definition.
static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *hData)
static int ProcessSources(const HeadModelT model, const uint dstRate, TokenReaderT *tr, HrirDataT *hData)
{
uint *setCount, *setFlag;
uint line, col, ei, ai;
uint res_points;
SourceRefT src;
double factor;
double *hrir;
ResamplerT rs;
ResamplerSetup(&rs, hData->mIrRate, dstRate);
/* Scale the number of IR points for resampling. This could be improved by
* also including space for the resampler build-up and fall-off (rs.mL*2),
* instead of clipping them off, but that could affect the HRTDs. It's not
* a big deal to exclude them for sources that aren't already minimum-
* phase).
*/
res_points = (uint)(((uint64)hData->mIrPoints*dstRate + hData->mIrRate-1) /
hData->mIrRate);
/* Clamp to the IR size to prevent overflow, and don't go less than the
* original point count.
*/
if(res_points > hData->mIrSize)
res_points = hData->mIrSize;
else if(res_points < hData->mIrPoints)
res_points = hData->mIrPoints;
setCount = (uint*)calloc(hData->mEvCount, sizeof(uint));
setFlag = (uint*)calloc(hData->mIrCount, sizeof(uint));
hrir = CreateArray(hData->mIrPoints);
hrir = CreateArray(res_points);
while(TrIsOperator(tr, "["))
{
TrIndication(tr, & line, & col);
@@ -2649,9 +2651,13 @@ static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *h
if(!LoadSource(&src, hData->mIrRate, hData->mIrPoints, hrir))
goto error;
if(hData->mIrRate != dstRate)
ResamplerRun(&rs, hData->mIrPoints, hrir,
res_points, hrir);
if(model == HM_DATASET)
AverageHrirOnset(hrir, 1.0 / factor, ei, ai, hData);
AverageHrirMagnitude(hrir, 1.0 / factor, ei, ai, hData);
AverageHrirOnset(hrir, res_points, 1.0 / factor, ei, ai, hData);
AverageHrirMagnitude(hrir, res_points, 1.0 / factor, ei, ai, hData);
factor += 1.0;
if(!TrIsOperator(tr, "+"))
break;
@@ -2660,6 +2666,8 @@ static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *h
setFlag[hData->mEvOffset[ei] + ai] = 1;
setCount[ei]++;
}
hData->mIrPoints = res_points;
hData->mIrRate = dstRate;
ei = 0;
while(ei < hData->mEvCount && setCount[ei] < 1)
@@ -2673,6 +2681,7 @@ static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *h
{
if(!TrLoad(tr))
{
ResamplerClear(&rs);
DestroyArray(hrir);
free(setFlag);
free(setCount);
@@ -2687,6 +2696,7 @@ static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *h
TrError(tr, "Missing source references.\n");
error:
ResamplerClear(&rs);
DestroyArray(hrir);
free(setFlag);
free(setCount);
@@ -2735,9 +2745,9 @@ static int ProcessDefinition(const char *inName, const uint outRate, const uint
fclose(fp);
return 0;
}
hData.mHrirs = CreateArray(hData.mIrCount * hData . mIrSize);
hData.mHrirs = CreateArray(hData.mIrCount * hData.mIrSize);
hData.mHrtds = CreateArray(hData.mIrCount);
if(!ProcessSources(model, &tr, &hData))
if(!ProcessSources(model, outRate ? outRate : hData.mIrRate, &tr, &hData))
{
DestroyArray(hData.mHrtds);
DestroyArray(hData.mHrirs);
@@ -2758,11 +2768,6 @@ static int ProcessDefinition(const char *inName, const uint outRate, const uint
}
fprintf(stdout, "Performing minimum phase reconstruction...\n");
ReconstructHrirs(&hData);
if(outRate != 0 && outRate != hData.mIrRate)
{
fprintf(stdout, "Resampling HRIRs...\n");
ResampleHrirs(outRate, &hData);
}
fprintf(stdout, "Truncating minimum-phase HRIRs...\n");
hData.mIrPoints = truncSize;
fprintf(stdout, "Synthesizing missing elevations...\n");