23 changes: 21 additions & 2 deletions src/hrtf/easy.c
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
#include <stdio.h>
#include <stdbool.h>
#include <assert.h>
#include <string.h>

#include "mysofa_export.h"
#include "mysofa.h"

Expand Down Expand Up @@ -138,9 +140,9 @@ MYSOFA_EXPORT void mysofa_getfilter_short(struct MYSOFA_EASY* easy, float x,
}
}

MYSOFA_EXPORT void mysofa_getfilter_float(struct MYSOFA_EASY* easy, float x,
MYSOFA_EXPORT void mysofa_getfilter_float_advanced(struct MYSOFA_EASY* easy, float x,
float y, float z, float *IRleft, float *IRright, float *delayLeft,
float *delayRight) {
float *delayRight, bool interpolate) {
float c[3];
float delays[2];
float *fl;
Expand All @@ -156,6 +158,11 @@ MYSOFA_EXPORT void mysofa_getfilter_float(struct MYSOFA_EASY* easy, float x,
assert(nearest >= 0);
neighbors = mysofa_neighborhood(easy->neighborhood, nearest);

// bypass interpolate by forcing current cooordinates to nearest's
if (!interpolate) {
memcpy(c, easy->hrtf->SourcePosition.values + nearest * easy->hrtf->C, sizeof(float) * easy->hrtf->C);
}

float *res = mysofa_interpolate(easy->hrtf, c, nearest, neighbors,
easy->fir, delays);

Expand All @@ -170,6 +177,18 @@ MYSOFA_EXPORT void mysofa_getfilter_float(struct MYSOFA_EASY* easy, float x,
}
}

MYSOFA_EXPORT void mysofa_getfilter_float(struct MYSOFA_EASY* easy, float x,
float y, float z, float *IRleft, float *IRright, float *delayLeft,
float *delayRight) {
mysofa_getfilter_float_advanced(easy, x, y, z, IRleft, IRright, delayLeft, delayRight, true);
}

MYSOFA_EXPORT void mysofa_getfilter_float_nointerp(struct MYSOFA_EASY* easy, float x,
float y, float z, float *IRleft, float *IRright, float *delayLeft,
float *delayRight) {
mysofa_getfilter_float_advanced(easy, x, y, z, IRleft, IRright, delayLeft, delayRight, false);
}

MYSOFA_EXPORT void mysofa_close(struct MYSOFA_EASY* easy) {
if (easy) {
if (easy->fir)
Expand Down
2 changes: 2 additions & 0 deletions src/hrtf/mysofa.h
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,8 @@ void mysofa_getfilter_short(struct MYSOFA_EASY* easy, float x, float y, float z,
short *IRleft, short *IRright, int *delayLeft, int *delayRight);
void mysofa_getfilter_float(struct MYSOFA_EASY* easy, float x, float y, float z,
float *IRleft, float *IRright, float *delayLeft, float *delayRight);
void mysofa_getfilter_float_nointerp(struct MYSOFA_EASY* easy, float x, float y, float z,
float *IRleft, float *IRright, float *delayLeft, float *delayRight);
void mysofa_close(struct MYSOFA_EASY* easy);
void mysofa_close_cached(struct MYSOFA_EASY* easy);

Expand Down
23 changes: 13 additions & 10 deletions src/hrtf/neighbors.c
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,11 @@ MYSOFA_EXPORT struct MYSOFA_NEIGHBORHOOD *mysofa_neighborhood_init_withstepdefin
int i, index;
float *origin, *test;
float radius, radius2;
float theta, theta2;
float phi, phi2;
float theta;
float phi;

// distance (degree) beyond which neighbor search is abandonned
float maxNeighborSearchAngle = 45;

struct MYSOFA_NEIGHBORHOOD *neighbor = malloc(
sizeof(struct MYSOFA_NEIGHBORHOOD));
Expand All @@ -53,7 +56,7 @@ MYSOFA_EXPORT struct MYSOFA_NEIGHBORHOOD *mysofa_neighborhood_init_withstepdefin
if ((lookup->phi_max - lookup->phi_min) > FLT_MIN) {
phi = angleStep;
do {
phi2 = test[0] = origin[0] + phi;
test[0] = origin[0] + phi;
test[1] = origin[1];
test[2] = origin[2];
convertSphericalToCartesian(test, 3);
Expand All @@ -63,11 +66,11 @@ MYSOFA_EXPORT struct MYSOFA_NEIGHBORHOOD *mysofa_neighborhood_init_withstepdefin
break;
}
phi += angleStep;
} while (phi2 <= lookup->phi_max + angleStep);
} while (phi <= maxNeighborSearchAngle);

phi = -angleStep;
do {
phi2 = test[0] = origin[0] + phi;
test[0] = origin[0] + phi;
test[1] = origin[1];
test[2] = origin[2];
convertSphericalToCartesian(test, 3);
Expand All @@ -77,14 +80,14 @@ MYSOFA_EXPORT struct MYSOFA_NEIGHBORHOOD *mysofa_neighborhood_init_withstepdefin
break;
}
phi -= angleStep;
} while (phi2 >= lookup->phi_min - angleStep);
} while (phi >= -maxNeighborSearchAngle);
}

if ((lookup->theta_max - lookup->theta_min) > FLT_MIN) {
theta = angleStep;
do {
test[0] = origin[0];
theta2 = test[1] = origin[1] + theta;
test[1] = origin[1] + theta;
test[2] = origin[2];
convertSphericalToCartesian(test, 3);
index = mysofa_lookup(lookup, test);
Expand All @@ -93,12 +96,12 @@ MYSOFA_EXPORT struct MYSOFA_NEIGHBORHOOD *mysofa_neighborhood_init_withstepdefin
break;
}
theta += angleStep;
} while (theta2 <= lookup->theta_max + angleStep);
} while (theta <= maxNeighborSearchAngle);

theta = -angleStep;
do {
test[0] = origin[0];
theta2 = test[1] = origin[1] + theta;
test[1] = origin[1] + theta;
test[2] = origin[2];
convertSphericalToCartesian(test, 3);
index = mysofa_lookup(lookup, test);
Expand All @@ -107,7 +110,7 @@ MYSOFA_EXPORT struct MYSOFA_NEIGHBORHOOD *mysofa_neighborhood_init_withstepdefin
break;
}
theta -= angleStep;
} while (theta2 >= lookup->theta_min - angleStep);
} while (theta >= -maxNeighborSearchAngle);
}

if ((lookup->radius_max - lookup->radius_min) > FLT_MIN) {
Expand Down
7 changes: 5 additions & 2 deletions src/hrtf/reader.c
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ static int getDimension(unsigned *dim, struct DATAOBJECT *dataobject) {
log(" %s=%s\n",attr->name,attr->value);

if (!strcmp(attr->name, "NAME")
&& !strncmp(attr->value,
&& attr->value && !strncmp(attr->value,
"This is a netCDF dimension but not a netCDF variable.",
53)) {
char *p = attr->value + strlen(attr->value) - 1;
Expand Down Expand Up @@ -189,7 +189,10 @@ static struct MYSOFA_HRTF *getHrtf(struct READER *reader, int *err) {
dir = reader->superblock.dataobject.directory;
while (dir) {

if (!strcmp(dir->dataobject.name, "ListenerPosition")) {
if(!dir->dataobject.name) {
log("SOFA VARIABLE IS NULL.\n");
}
else if (!strcmp(dir->dataobject.name, "ListenerPosition")) {
*err = getArray(&hrtf->ListenerPosition, &dir->dataobject);
} else if (!strcmp(dir->dataobject.name, "ReceiverPosition")) {
*err = getArray(&hrtf->ReceiverPosition, &dir->dataobject);
Expand Down
4 changes: 2 additions & 2 deletions src/hrtf/resample.c
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,8 @@ MYSOFA_EXPORT int mysofa_resample(struct MYSOFA_HRTF *hrtf, float samplerate) {
/*
* update delay values
*/
for (i = 0; i < hrtf->DataIR.elements; i++)
hrtf->DataIR.values[i] /= factor;
for (i = 0; i < hrtf->DataDelay.elements; i++)
hrtf->DataDelay.values[i] *= factor;

/*
* update sample rate
Expand Down
2 changes: 1 addition & 1 deletion src/tests/json.c
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ static void printDimensions(FILE *out, struct MYSOFA_HRTF *hrtf,

fprintf(out, " \"DimensionNames\":[");
s = found->value;
while (s[0] && dims < 4) {
while (s && s[0] && dims < 4) {
switch (s[0]) {
case 'I':
dimensions[dims++] = hrtf->I;
Expand Down
60 changes: 46 additions & 14 deletions src/tests/resample.c
Original file line number Diff line number Diff line change
Expand Up @@ -10,40 +10,72 @@
void test_resample() {
struct MYSOFA_HRTF *hrtf;
int err = 0, i;
float *backup;
float *backupIr;
float *backupDelays;
int resampFactor = 2; // integer, >= 1
int num_ir_to_test = 10; // number of filters to test in the IR, to avoid going through all potential positions / channels
float irPeakValue = 0.0;

hrtf = mysofa_load("tests/Pulse.sofa", &err);
// hrtf = mysofa_load("tests/Pulse.sofa", &err);
hrtf = mysofa_load("tests/CIPIC_subject_003_hrir_final_itdInDelayField.sofa", &err);
if (!hrtf) {
CU_FAIL_FATAL("Error reading file.");
return;
}

backup = malloc(sizeof(float) * hrtf->N * 3);
if (!backup) {
CU_FAIL_FATAL("No memory, N is too large.");
// backup (non resampled) IR
backupIr = malloc(sizeof(float) * hrtf->N * num_ir_to_test);
if (!backupIr) {
CU_FAIL_FATAL("No memory, N (IR) is too large.");
mysofa_free(hrtf);
return;
}

for (i = 0; i < hrtf->N * 3; i++) {
backup[i] = hrtf->DataIR.values[i];
for (i = 0; i < hrtf->N * num_ir_to_test; i++) {
backupIr[i] = hrtf->DataIR.values[i];
irPeakValue = fmax(irPeakValue, backupIr[i]);
}

err = mysofa_resample(hrtf, 96000.);
// backup (non resampled) Delays
backupDelays = malloc(sizeof(float) * hrtf->DataDelay.elements);
if (!backupDelays) {
CU_FAIL_FATAL("No memory, N (Delay) is too large.");
mysofa_free(hrtf);
return;
}

for (i = 0; i < hrtf->DataDelay.elements; i++) {
backupDelays[i] = hrtf->DataDelay.values[i];
}

float fsOld = hrtf->DataSamplingRate.values[0];
float fsNew = resampFactor * fsOld;
err = mysofa_resample(hrtf, fsNew);
CU_ASSERT_FATAL(err == MYSOFA_OK);

for (i = 0; i < hrtf->N * 3; i++) {
// loop over delays (in samples), check resampling factor correctly applied
float delayThresh = 0.001; // in samples, maximum acceptable error threshold
for (i = 0; i < hrtf->DataDelay.elements; i+= 2*resampFactor) {
#ifdef VDEBUG
printf("%f %f ", backupDelays[i]/fsOld, hrtf->DataDelay.values[i]/fsNew);
#endif
CU_ASSERT( fabs( backupDelays[i]/fsOld - hrtf->DataDelay.values[i]/fsNew) <= delayThresh);
}

// loop over IR to check if paired values still match after resampling
float irThresh = 0.001 * irPeakValue; // relative acceptable error threshold
for (i = 0; i < hrtf->N * num_ir_to_test; i+=resampFactor) {

#ifdef VDEBUG
printf("%6.3f~%6.3f ", hrtf->DataIR.values[i], backup[i/2]);
printf("%f %f ", backupIr[i/resampFactor], hrtf->DataIR.values[i]);
if ((i % hrtf->N) == (hrtf->N - 1))
printf("\n");
#endif
CU_ASSERT(
!((hrtf->DataIR.values[i] > 0.4 && backup[i / 2] == 0.)
|| (hrtf->DataIR.values[i] <= 0.3 && backup[i / 2] == 1.)));
CU_ASSERT( ( fabs(hrtf->DataIR.values[i] - backupIr[i / resampFactor]) <= irThresh ) ) ;
}

free(backup);
free(backupIr);
free(backupDelays);
mysofa_free(hrtf);
}

Binary file not shown.