Convert filtering method to frequency domain

- Add evalIIRTransferFunction and associated unit tests
- Processing function to only update coeffs on changed distatance/direction.
- BInauraliser: remove unused var src_dirs.
- add `src_dirs_cur` pointer to switch between non/rotated directions.
- `doaToIpsiInteraural` returns lateral angle `alpha` to properly parameterize the DVF with changing elevation.
- Add binauraliser_nf to README.md

Author: mtmccrea

README.md

@@ -156,6 +156,7 @@ Several **examples** have also been included in the repository, which may serve
 * **array2sh** - converts microphone array signals into spherical harmonic signals (aka Ambisonic signals), based on theoretical descriptions of the array configuration and construction.
 * **beamformer** - a beamformer/virtual microphone generator for Ambisonic signals, with several different beam pattern options.
 * **binauraliser** - convolves input audio with interpolated HRTFs, which can be optionally loaded from a SOFA file.
+* **binauraliser_nf** - binauraliser, with the addition of proximity filtering for near field sound sources.
 * **decorrelator** - a basic multi-channel signal decorrelator.
 * **dirass** - a sound-field visualiser based on re-assigning the energy of beamformers. This re-assignment is based on the DoA estimates extracted from spatially-localised active-intensity vectors, which are biased towards each beamformer direction.
 * **matrixconv** - a basic matrix convolver with an optional partitioned convolution mode. 

examples/include/binauraliser_nf.h

@@ -132,6 +132,17 @@ void binauraliserNF_process(void* const hBin,
                           int nOutputs,
                           int nSamples);
 
+/**
+ * Alternate version of binauraliserNF_process() that performs frequency-domain
+ * DVF filtering. Not used but kept for posterity.
+ */
+void binauraliserNF_processFD(void* const hBin,
+                          const float *const * inputs,
+                          float** const outputs,
+                          int nInputs,
+                          int nOutputs,
+                          int nSamples);
+
 
 /* ========================================================================== */
 /*                                Set Functions                               */
@@ -170,12 +181,14 @@ void binauraliserNF_setInputConfigPreset(void* const hBin,
 float binauraliserNF_getSourceDist_m(void* const hBin, int index);
 
 /**
- * Returns the distance considered to be the far field (beyond which no near field filtering is applied), in METERS
+ * Returns the distance considered to be the far field (beyond which no near
+ * field filtering is applied), in METERS
  */
 float binauraliserNF_getFarfieldThresh_m(void* const hBin);
 
 /**
- * Returns the scaling factor to give the far field threshold headroom (useful for UI range limits)
+ * Returns the scaling factor to give the far field threshold headroom (useful
+ * for UI range limits)
  */
 float binauraliserNF_getFarfieldHeadroom(void* const hBin);
 

examples/src/binauraliser/binauraliser.c

@@ -200,14 +200,13 @@ void binauraliser_process
 {
     binauraliser_data *pData = (binauraliser_data*)(hBin);
     int ch, ear, i, band, nSources;
-    float src_dirs[MAX_NUM_INPUTS][2], Rxyz[3][3], hypotxy;
+    float Rxyz[3][3], hypotxy;
     int enableRotation;
 
     /* copy user parameters to local variables */
     nSources = pData->nSources;
     enableRotation = pData->enableRotation;
-    memcpy(src_dirs, pData->src_dirs_deg, MAX_NUM_INPUTS*2*sizeof(float));
-
+    
     /* apply binaural panner */
     if ((nSamples == BINAURALISER_FRAME_SIZE) && (pData->hrtf_fb!=NULL) && (pData->codecStatus==CODEC_STATUS_INITIALISED) ){
         pData->procStatus = PROC_STATUS_ONGOING;

examples/src/binauraliser/binauraliser_internal.h

@@ -67,7 +67,8 @@ extern "C" {
 
 /**
  * Main structure for binauraliser. Contains variables for audio buffers,
- * afSTFT, HRTFs, internal variables, flags, user parameters
+ * afSTFT, HRTFs, internal variables, flags, user parameters.
+ * Note: if this is modified, identically modify _binauraliserNF struct.
  */
 typedef struct _binauraliser
 {

examples/src/binauraliser_nf/binauraliser_nf.c

@@ -64,12 +64,13 @@ extern "C" {
  * Main structure for binauraliserNF. Contains all variables from binauraliser
  * (audio buffers, afSTFT, HRTFs, internal variables, flags, user parameters) plus
  * those specific to the near field variant.
+ * FREQUENCY DOMAIN implementation.
  */
 typedef struct _binauraliserNF {
-    //struct _binauraliser;               /**< "inherit" member vars of binauraliser struct */
+    //struct _binauraliser; /**< "inherit" member vars of binauraliser struct - requires -fms-extensions compile flag */
 
     /* The following variables MUST match those of the _binauraliser struct */
-    
+
     /* audio buffers */
     float** inputFrameTD;            /**< time-domain input frame; #MAX_NUM_INPUTS x #BINAURALISER_FRAME_SIZE */
     float** outframeTD;              /**< time-domain output frame; #NUM_EARS x #BINAURALISER_FRAME_SIZE */
@@ -78,7 +79,7 @@ typedef struct _binauraliserNF {
     int fs;                          /**< Host sampling rate, in Hz */
     float freqVector[HYBRID_BANDS];  /**< Frequency vector (filterbank centre frequencies) */
     void* hSTFT;                     /**< afSTFT handle */
-    
+
     /* sofa file info */
     char* sofa_filepath;             /**< absolute/relevative file path for a sofa file */
     float* hrirs;                    /**< time domain HRIRs; FLAT: N_hrir_dirs x #NUM_EARS x hrir_len */
@@ -89,19 +90,19 @@ typedef struct _binauraliserNF {
     int hrir_loaded_fs;              /**< sampling rate of the loaded HRIRs  */
     int hrir_runtime_fs;             /**< sampling rate of the HRIRs being used for processing (after any resampling) */
     float* weights;                  /**< Integration weights for the HRIR measurement grid */
-    
+
     /* vbap gain table */
     int hrtf_vbapTableRes[2];        /**< [0] azimuth, and [1] elevation grid resolution, in degrees */
     int N_hrtf_vbap_gtable;          /**< Number of interpolation weights/directions */
     int* hrtf_vbap_gtableIdx;        /**< N_hrtf_vbap_gtable x 3 */
     float* hrtf_vbap_gtableComp;     /**< N_hrtf_vbap_gtable x 3 */
-    
+
     /* hrir filterbank coefficients */
     float* itds_s;                   /**< interaural-time differences for each HRIR (in seconds); nBands x 1 */
     float_complex* hrtf_fb;          /**< hrtf filterbank coefficients; nBands x nCH x N_hrirs */
     float* hrtf_fb_mag;              /**< magnitudes of the hrtf filterbank coefficients; nBands x nCH x N_hrirs */
     float_complex hrtf_interp[MAX_NUM_INPUTS][HYBRID_BANDS][NUM_EARS]; /**< Interpolated HRTFs */
-    
+
     /* flags/status */
     CODEC_STATUS codecStatus;        /**< see #CODEC_STATUS */
     float progressBar0_1;            /**< Current (re)initialisation progress, between [0..1] */
@@ -110,7 +111,7 @@ typedef struct _binauraliserNF {
     int recalc_hrtf_interpFLAG[MAX_NUM_INPUTS]; /**< 1: re-calculate/interpolate the HRTF, 0: do not */
     int reInitHRTFsAndGainTables;    /**< 1: reinitialise the HRTFs and interpolation tables, 0: do not */
     int recalc_M_rotFLAG;            /**< 1: re-calculate the rotation matrix, 0: do not */
-    
+
     /* misc. */
     float src_dirs_rot_deg[MAX_NUM_INPUTS][2]; /**< Intermediate rotated source directions, in degrees */
     float src_dirs_rot_xyz[MAX_NUM_INPUTS][3]; /**< Intermediate rotated source directions, as unit-length Cartesian coordinates */
@@ -133,23 +134,29 @@ typedef struct _binauraliserNF {
     int bFlipRoll;                           /**< flag to flip the sign of the roll rotation angle */
     int useRollPitchYawFlag;                 /**< rotation order flag, 1: r-p-y, 0: y-p-r */
     float src_gains[MAX_NUM_INPUTS];         /**< Gains applied per source */
-    
+
     /* End copied _binauraliser struct members. The following are unique to the _binauraliserNF struct */
-    
-    /* audio buffers */
-    float**          binsrcsTD;         /**< near field DVF-filtered sources frame; (#MAX_NUM_INPUTS * #NUM_EARS) x #BINAURALISER_FRAME_SIZE */
-    float_complex*** binauralTF;        /**< time-frequency domain output frame; #HYBRID_BANDS x (#MAX_NUM_INPUTS * #NUM_EARS) x #TIME_SLOTS
-                                         *   Note: (#MAX_NUM_INPUTS * #NUM_EARS) dimensions are combined because afSTFT requires type float_complex*** */
+
+    float b_dvf[MAX_NUM_INPUTS][NUM_EARS][2];                   /**< shelf IIR numerator coefficients for each input, left and right. */
+    float a_dvf[MAX_NUM_INPUTS][NUM_EARS][2];                   /**< shelf IIR denominator coefficients for each input, left and right. */
+    float dvfmags[MAX_NUM_INPUTS][NUM_EARS][HYBRID_BANDS];      /**< DVF filter frequency band magnitudes. */
+    float dvfphases[MAX_NUM_INPUTS][NUM_EARS][HYBRID_BANDS];    /**< DVF filter frequency band phases. */
+
     /* misc. */
-    float src_dists_m[MAX_NUM_INPUTS];  /**< source distance,  meters */
-    float farfield_thresh_m;            /**< distance considered to be far field (no near field filtering),  meters */
-    float farfield_headroom;            /**< scale factor applied to farfield_thresh_m when resetting to the far field, and for UI range, meters */
-    float nearfield_limit_m;            /**< minimum distance allowed for near-field filtering, from head _center_,  meters, def. 0.15 */
-    float head_radius;                  /**< head radius, used calculate normalized source distance meters, def. 0.09096 */
-    float head_radius_recip;            /**< reciprocal of head radius */
+    float src_dists_m[MAX_NUM_INPUTS];  /**< Source distance,  meters. */
+    float farfield_thresh_m;            /**< Distance considered to be far field (no near field filtering),  meters. */
+    float farfield_headroom;            /**< Scale factor applied to farfield_thresh_m when resetting to the far field, and for UI range, meters. */
+    float nearfield_limit_m;            /**< Minimum distance allowed for near-field filtering, from head _center_, meters, def. 0.15. */
+    float head_radius;                  /**< Head radius, used calculate normalized source distance meters, def. 0.09096. */
+    float head_radius_recip;            /**< Reciprocal of head radius. */
+    float (*src_dirs_cur)[2];           /**< Pointer to assign to the current HRTF directions being operated on (non/rotated directions switch). */
+
+    /* flags/status */
+    int recalc_dvfCoeffFLAG[MAX_NUM_INPUTS]; /**< 1: re-calculate the DVF coefficients on change in distance, 0: do not. */
 
 } binauraliserNF_data;
 
+
 /* ========================================================================== */
 /*                             Internal Functions                             */
 /* ========================================================================== */

examples/src/binauraliser_nf/binauraliser_nf_internal.h

@@ -48,12 +48,12 @@ static const double p23[19] = { -0.67, 0.142, 3404., -0.91, -0.67, -1.21, -1.76,
 static const double p33[19] = { 0.174, -0.11, -1699., 0.437, 0.658, 2.02, 6.815, 0.614, 589.3, 16818., -9.39, -44.1, -23.6, -92.3, -1.81, 10.54, 0.532, 0.285, -2.08 };
 static const double q13[19] = { -1.75, -0.01, 7354., -2.18, -1.2, -1.59, -1.23, -0.89, 29.23, 1945., -0.06, 5.635, 3.308, 13.88, -0.88, -2.23, -0.96, -0.9, -0.57 };
 static const double q23[19] = { 0.699, -0.35, -5350., 1.188, 0.256, 0.816, 1.166, 0.76, 59.51, 1707., -1.12, -6.18, -3.39, -12.7, -0.19, 1.295, -0.02, -0.08, -0.4 };
-
 static const int numAz_table = sizeof(q23);
-//static const float a_0 = 0.0875f;       /**< Reference head size, 8.75 centimeters, used in the generation of the coeff lookup table. */
-//static const float a_head = 0.09096f;   /**< This head size, See note for head_radius in binauraliser_nf. */
-static const float headDim = SAF_PI * (0.0875f/*a_0*/ / 0.09096f /*a_head*/);
-static const float sosDiv2PiA = 343.0f / (2.0f * SAF_PI * 0.09096f /*a_head*/);
+
+/* a_0 = 0.0875f; Reference head size, 8.75 centimeters, used in the generation of the coeff lookup table. */
+/* a_head = 0.09096f; This head size, see note for head_radius in binauraliser_nf. */
+static const float headDim = SAF_PI * (0.0875f / 0.09096f);             /* pi * (a_0 / a_head) */
+static const float sosDiv2PiA = 343.0f / (2.0f * SAF_PI * 0.09096f);    /*  c / (2pi * a_head) */
 
 /** Linear interpolation between two values */
 static float interpolate_lin
@@ -78,18 +78,18 @@ static float db2mag
 /*
  * Calculate high-shelf parameters, g0, gInf, fc, from the lookup table coefficients (10 degree steps).
  * Called twice per update as the returned values are subsequently interpolated to exact azimuth. */
-void calcHighShelfParams
+void calcDVFShelfParams
 (
-    int i,          /* index into the coefficient table, dictated by azimuth */
-    float rhoIn,    /* normalized source distance */
-    float* g0,      /* high shelf gain at DC */
-    float* gInf,    /* high shelf gain at inf */
-    float* fc       /* high shelf cutoff frequency */
+    int i,          /* Index into the coefficient table, dictated by azimuth. */
+    float rhoIn,    /* Normalized source distance. */
+    float* g0,      /* High shelf gain at DC. */
+    float* gInf,    /* High shelf gain at inf. */
+    float* fc       /* High shelf cutoff frequency. */
 )
 {
     float fc_tmp;
     double rho = (double)rhoIn;
-	double rhoSq = rho*rho;
+	double rhoSq = rho * rho;
 
     /*  Eq (8), (13) and (14) */
     *g0    = (float)((p11[i] * rho   + p21[i]) / (rhoSq + q11[i] * rho + q21[i]));
@@ -101,69 +101,65 @@ void calcHighShelfParams
 }
 
 /*
- * Interpolate (linear) the high shelf parameters generated by calcHighShelfParams()
- * which is called twice to generate the high shelf parameters for the nearest thetas
- * in the lookup table. */
-void interpHighShelfParams
+ * Interpolate (linearly) the high shelf parameters generated by
+ * calcDVFShelfParams() which is called twice to generate the high shelf
+ * parameters for the nearest thetas in the lookup table. */
+void interpDVFShelfParams
 (
-    float theta,   /* ipsilateral azimuth, on the inter-aural axis [0, 180] (deg) */
-    float rho,     /* distance, normalized to head radius, >= 1 */
+    float theta,   /* Ipsilateral azimuth, on the inter-aural axis [0, 180] (deg). */
+    float rho,     /* Distance, normalized to head radius, >= 1. */
     /* output */
     float* iG0,
     float* iGInf,
     float* iFc
 )
 {
     int theta_idx_lower, theta_idx_upper;
-    float ifac;
-    float thetaDiv10;
+    float ifac, thetaDiv10;
     float g0_1, g0_2;       /* high shelf gain at DC */
     float gInf_1, gInf_2;   /* high shelf gain at inf */
     float fc_1, fc_2;       /* high shelf cutoff frequency */
 
-    // TBD: add range checking? - clip theta and rho to valid range
+    // TBD: could add range checking, clipping theta and rho to valid range
 
-    /* linearly interpolate DC gain, HF gain, center freq at theta */
+    /* Linearly interpolate DC gain, HF gain, center freq at theta.
+     * Table is in 10 degree steps, floor(x/10) gets lower index. */
     thetaDiv10 = theta / 10.f;
-    theta_idx_lower = (int)thetaDiv10;      /* Table is in 10 degree steps, floor(x/10) gets lower index */
+    theta_idx_lower = (int)thetaDiv10;
     theta_idx_upper = theta_idx_lower + 1;
-    if(theta_idx_upper == numAz_table) {    // alternatively, instead check theta_idx_upper => numAz_table, could clip the value > 180 here
+    if(theta_idx_upper == numAz_table) {    // could alternatively check theta_idx_upper => numAz_table, clip the value > 180 here
         theta_idx_upper = theta_idx_lower;
         theta_idx_lower = theta_idx_lower - 1;
     }
 
-    calcHighShelfParams(theta_idx_lower, rho, &g0_1, &gInf_1, &fc_1);
-    calcHighShelfParams(theta_idx_upper, rho, &g0_2, &gInf_2, &fc_2);
+    calcDVFShelfParams(theta_idx_lower, rho, &g0_1, &gInf_1, &fc_1);
+    calcDVFShelfParams(theta_idx_upper, rho, &g0_2, &gInf_2, &fc_2);
 
-    ifac = thetaDiv10 - theta_idx_lower;  /* interpolation factor between table steps */
+    /* interpolation factor between table steps */
+    ifac   = thetaDiv10 - theta_idx_lower;
     *iG0   = interpolate_lin(g0_1,   g0_2,   ifac);
     *iGInf = interpolate_lin(gInf_1, gInf_2, ifac);
     *iFc   = interpolate_lin(fc_1,   fc_2,   ifac);
 }
 
 /*
- * Generate IIR coefficients from the shelf parameters generated by calcHighShelfParams() and
- * interpHighShelfParams(). */
-void calcIIRCoeffs
+ * Calculate the DVF filter coefficients from shelving filter parameters.
+ */
+void dvfShelfCoeffs
 (
     /* Input */
-    float g0,      /* high shelf dc gain */
-    float gInf,    /* high shelf high gain */
-    float fc,      /* high shelf center freq */
-    float fs,      /* sample rate */
+    float g0,      /* High shelf dc gain */
+    float gInf,    /* High shelf high gain */
+    float fc,      /* High shelf center freq */
+    float fs,      /* Sample rate */
     /* Output */
     float* b0,     /* IIR coeffs */
     float* b1,
     float* a1
 )
 {
-    float v0;
-    float g0_mag;
-    float tanF;
-    float v0tanF;
-    float a_c;
-    float v;
-    float va_c;
+    float v0, g0_mag, tanF, v0tanF;
+    float a_c, v, va_c;
 
     v0     = db2mag(gInf);              /* Eq. (12), (10), and (11) */
     g0_mag = db2mag(g0);                // optim TBD: revisit; does g0, gInf need to be in dB?
@@ -178,40 +174,56 @@ void calcIIRCoeffs
     *a1  = a_c;
 }
 
-void applyDVF
+void calcDVFCoeffs
 (
-    float theta,        /* ipsilateral azimuth, on the inter-aural axis [0, 180] (deg) */
-    float rho,          /* distance, normalized to head radius, >= 1 */
-    float* in_signal,
-    int nSamples,       /* Number of samples to process */
-    float fs,           /* Sample rate */
-    float* wz,          /* (&) Filter coefficients to be passed to the next block */
-    float* out_signal
+    float alpha,    /* Lateral angle, on the inter-aural axis [0, 180] (deg) */
+    float rho,      /* Distance, normalized to head radius, >= 1 */
+    float fs,       /* Sample rate */
+    float* b,
+    float* a
 )
 {
-    float b[2] = {0.f, 0.f};
-    float a[2] = {1.f, 0.f};
     float iG0, iGInf, iFc;
-    
-    interpHighShelfParams(theta, rho, &iG0, &iGInf, &iFc);
-    calcIIRCoeffs(iG0, iGInf, iFc, fs, &b[0], &b[1], &a[1]);
-    applyIIR(in_signal, nSamples, 2, &b[0], &a[0], wz, out_signal);
+    interpDVFShelfParams(alpha, rho, &iG0, &iGInf, &iFc);
+    dvfShelfCoeffs(iG0, iGInf, iFc, fs, &b[0], &b[1], &a[1]);
 }
 
-void convertFrontalDoAToIpsilateral
+void doaToIpsiInteraural
 (
-    float thetaFront, /* DOA wrt 0˚ forward-facing  [deg, (-180, 180)] */
-    float* ipsiDoaLR  /* pointer to 2-element array of left and right ear DoAs wrt interaural axis */
+    float azimuth,      /* azimuth wrt 0˚ forward-facing  (deg, [-360, 360]) */
+    float elevation,    /* elevation from the horizontal plane (deg, [-180, 180]) */
+    float* alphaLR,     /* (&) 2-element array of lateral angle alpha for left and right ear (deg, [0,180]) */
+    float* betaLR       /* (&) 2-element array of vertical angle beta for left and right ear (deg, [0,90]) */
 )
 {
-    float thetaL;
+    float az_rad, el_rad, alpha_ipsi, beta_ipsi, alpha_ipsi_deg, beta_ipsi_deg;
+    float sinaz, sinel, cosaz, cosel;
 
-    thetaFront = SAF_CLAMP(thetaFront, -180.f, 180.f);
-    thetaL = fabsf(90.f - thetaFront);
-    if (thetaL > 180.f) {
-        thetaL = 360.f - thetaL;
+    az_rad = DEG2RAD(azimuth);
+    el_rad = DEG2RAD(elevation);
+    sinaz = sinf(az_rad);
+    sinel = sinf(el_rad);
+    cosaz = cosf(az_rad);
+    cosel = cosf(el_rad);
+    alpha_ipsi = SAF_PI/2.f - acosf(sinaz * cosel);
+    beta_ipsi  = asinf(sinel / sqrtf(powf(sinel,2.f) + (powf(cosaz,2.f) * powf(cosel,2.f))));
+    
+    /* Convert interaural-polar coordinates back to spherical _ranges_. */
+    if(beta_ipsi > SAF_PI/2.f) {
+        alpha_ipsi = SAF_PI - alpha_ipsi;   /* [0,pi] */
+        beta_ipsi = SAF_PI - beta_ipsi;     /* [0,pi/2] */
     }
+
+    /* Convert to ipsilateral offset from the left ear: [-360, 360]->[0, 180] */
+    alpha_ipsi = fabsf(SAF_PI/2.f - alpha_ipsi);
+    if (alpha_ipsi > SAF_PI)
+        alpha_ipsi = 2*SAF_PI - alpha_ipsi;
 
-    ipsiDoaLR[0] = thetaL;
-    ipsiDoaLR[1] = 180.f - thetaL; // thetaR
+    alphaLR[0] = alpha_ipsi_deg = RAD2DEG(alpha_ipsi);
+    alphaLR[1] = 180.f - alpha_ipsi_deg;
+    
+    if (betaLR != NULL) {
+        betaLR[0] = beta_ipsi_deg = RAD2DEG(beta_ipsi);
+        betaLR[1] = 180.f - beta_ipsi_deg;
+    }
 }