diff --git a/README.md b/README.md index adf394dba..99cf09296 100644 --- a/README.md +++ b/README.md @@ -156,7 +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** - same as the binaural example, except including a distance control with near-field/proximity effect support. +* **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. diff --git a/examples/include/binauraliser_nf.h b/examples/include/binauraliser_nf.h index 1f7716d25..a426db07e 100644 --- a/examples/include/binauraliser_nf.h +++ b/examples/include/binauraliser_nf.h @@ -128,6 +128,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 */ @@ -166,12 +177,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); diff --git a/examples/src/binauraliser/binauraliser.c b/examples/src/binauraliser/binauraliser.c index 977e16290..e8025b134 100644 --- a/examples/src/binauraliser/binauraliser.c +++ b/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; diff --git a/examples/src/binauraliser/binauraliser_internal.h b/examples/src/binauraliser/binauraliser_internal.h index 4ba85e2d2..d311a784d 100644 --- a/examples/src/binauraliser/binauraliser_internal.h +++ b/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 { diff --git a/examples/src/binauraliser_nf/binauraliser_nf.c b/examples/src/binauraliser_nf/binauraliser_nf.c index 1a046605c..773ce1c08 100644 --- a/examples/src/binauraliser_nf/binauraliser_nf.c +++ b/examples/src/binauraliser_nf/binauraliser_nf.c @@ -38,14 +38,14 @@ #include "binauraliser_nf_internal.h" -void binauraliserNF_create +void binauraliserNF_create /* FREQUENCY DOMAIN version */ ( void ** const phBin ) { binauraliserNF_data* pData = (binauraliserNF_data*)malloc1d(sizeof(binauraliserNF_data)); *phBin = (void*)pData; - int ch, nDim; // nDim not actually used + int nDim; // nDim not actually used /* user parameters */ pData->useDefaultHRIRsFLAG = 1; /* pars->sofa_filepath must be valid to set this to 0 */ @@ -90,16 +90,13 @@ void binauraliserNF_create pData->weights = NULL; pData->N_hrir_dirs = pData->hrir_loaded_len = pData->hrir_runtime_len = 0; pData->hrir_loaded_fs = pData->hrir_runtime_fs = -1; /* unknown */ + /* time domain buffers */ - pData->inputFrameTD = (float**)calloc2d(MAX_NUM_INPUTS, BINAURALISER_FRAME_SIZE, sizeof(float)); - pData->binsrcsTD = (float**)calloc2d(MAX_NUM_INPUTS * NUM_EARS, BINAURALISER_FRAME_SIZE, sizeof(float)); - pData->outframeTD = NULL; /* outframeTD inherited from _binauraliser but not used by binauraliserNF functions */ - /* frequency domain buffers */ - pData->inputframeTF = (float_complex***)calloc3d(HYBRID_BANDS, MAX_NUM_INPUTS, TIME_SLOTS, sizeof(float_complex)); - pData->binauralTF = (float_complex***)calloc3d(HYBRID_BANDS, MAX_NUM_INPUTS * NUM_EARS, TIME_SLOTS, sizeof(float_complex)); /* MAX_NUM_INPUTS*NUM_EARS dimensions are combined because afSTFT requires float_complex*** */ - pData->outputframeTF = NULL; /* outputframeTF inherited from _binauraliser but not used by binauraliserNF functions */ + pData->inputFrameTD = (float**)malloc2d(MAX_NUM_INPUTS, BINAURALISER_FRAME_SIZE, sizeof(float)); + pData->outframeTD = (float**)malloc2d(NUM_EARS, BINAURALISER_FRAME_SIZE, sizeof(float)); + pData->inputframeTF = (float_complex***)malloc3d(HYBRID_BANDS, MAX_NUM_INPUTS, TIME_SLOTS, sizeof(float_complex)); + pData->outputframeTF = (float_complex***)malloc3d(HYBRID_BANDS, NUM_EARS, TIME_SLOTS, sizeof(float_complex)); - /* vbap (amplitude normalised) */ pData->hrtf_vbap_gtableIdx = NULL; pData->hrtf_vbap_gtableComp = NULL; @@ -109,6 +106,17 @@ void binauraliserNF_create pData->itds_s = NULL; pData->hrtf_fb = NULL; pData->hrtf_fb_mag = NULL; + + /* Initialize DVF filter parameters */ + memset(FLATTEN3D(pData->dvfmags), 1.f, MAX_NUM_INPUTS * NUM_EARS * HYBRID_BANDS * sizeof(float)); + memset(FLATTEN3D(pData->dvfphases), 0.f, MAX_NUM_INPUTS * NUM_EARS * HYBRID_BANDS * sizeof(float)); + memset(FLATTEN3D(pData->b_dvf), 0.f, MAX_NUM_INPUTS * NUM_EARS * 2 * sizeof(float)); + memset(FLATTEN3D(pData->a_dvf), 0.f, MAX_NUM_INPUTS * NUM_EARS * 2 * sizeof(float)); + for(int ch = 0; ch < MAX_NUM_INPUTS; ch++) { + for(int ear = 0; ear < NUM_EARS; ear++) { + pData->a_dvf[ch][ear][0] = 1.f; /* a_0 = 1.0, always */ + } + } /* flags/status */ pData->progressBar0_1 = 0.0f; @@ -117,11 +125,14 @@ void binauraliserNF_create pData->codecStatus = CODEC_STATUS_NOT_INITIALISED; pData->procStatus = PROC_STATUS_NOT_ONGOING; pData->reInitHRTFsAndGainTables = 1; - for(ch=0; chrecalc_hrtf_interpFLAG[ch] = 1; + pData->recalc_dvfCoeffFLAG[ch] = 1; pData->src_gains[ch] = 1.f; } pData->recalc_M_rotFLAG = 1; + + pData->src_dirs_cur = pData->src_dirs_deg; } void binauraliserNF_destroy @@ -144,8 +155,6 @@ void binauraliserNF_destroy free(pData->outframeTD); free(pData->inputframeTF); free(pData->outputframeTF); - free(pData->binsrcsTD); - free(pData->binauralTF); free(pData->hrtf_vbap_gtableComp); free(pData->hrtf_vbap_gtableIdx); free(pData->hrtf_fb); @@ -197,8 +206,9 @@ void binauraliserNF_initCodec pData->progressBar0_1 = 0.0f; /* check if TFT needs to be reinitialised */ - binauraliserNF_initTFT(hBin); - + // use binauraliser_initTFT() Freq-domain DVF, binauraliserNF_initTFT() for Time-domain DVF + binauraliser_initTFT(hBin); + /* reinit HRTFs and interpolation tables */ if(pData->reInitHRTFsAndGainTables){ binauraliser_initHRTFsAndGainTables(hBin); @@ -209,10 +219,9 @@ void binauraliserNF_initCodec strcpy(pData->progressBarText,"Done!"); pData->progressBar0_1 = 1.0f; pData->codecStatus = CODEC_STATUS_INITIALISED; - } -void binauraliserNF_process +void binauraliserNF_process /* FREQ DOMAIN version */ ( void * const hBin, const float *const * inputs, @@ -223,17 +232,12 @@ void binauraliserNF_process ) { binauraliserNF_data *pData = (binauraliserNF_data*)(hBin); - int ch, ear, i, band, nSources, srci; - float src_dirs[MAX_NUM_INPUTS][2], src_dists[MAX_NUM_INPUTS], Rxyz[3][3]; - float hypotxy, headRadiusRecip, fs, ffThresh; - int enableRotation; - float thetaLR[2] = { 0.0, 0.0 }; - float rho, wzL, wzR; + int ch, ear, i, band, nSources, enableRotation; + float hypotxy, headRadiusRecip, fs, ffThresh, rho; + float Rxyz[3][3]; + float alphaLR[2] = { 0.0, 0.0 }; /* copy user parameters to local variables */ - memcpy(src_dirs, pData->src_dirs_deg, MAX_NUM_INPUTS * sizeof(float)*2); - memcpy(src_dists, pData->src_dists_m, MAX_NUM_INPUTS * sizeof(float)); - nSources = pData->nSources; enableRotation = pData->enableRotation; headRadiusRecip = pData->head_radius_recip; @@ -244,10 +248,25 @@ void binauraliserNF_process if ((nSamples == BINAURALISER_FRAME_SIZE) && (pData->hrtf_fb!=NULL) && (pData->codecStatus==CODEC_STATUS_INITIALISED) ) { pData->procStatus = PROC_STATUS_ONGOING; + /* Load time-domain data */ + for (i = 0; i < SAF_MIN(nSources, nInputs); i++) + utility_svvcopy(inputs[i], BINAURALISER_FRAME_SIZE, pData->inputFrameTD[i]); + for (; i < nSources; i++) + memset(pData->inputFrameTD[i], 0, BINAURALISER_FRAME_SIZE * sizeof(float)); + + /* Apply source gains */ + for (ch = 0; ch < nSources; ch++) { + if(fabsf(pData->src_gains[ch] - 1.f) > 1e-6f) + utility_svsmul(pData->inputFrameTD[ch], &(pData->src_gains[ch]), BINAURALISER_FRAME_SIZE, NULL); + } + + /* Apply time-frequency transform (TFT) */ + afSTFT_forward_knownDimensions(pData->hSTFT, pData->inputFrameTD, BINAURALISER_FRAME_SIZE, MAX_NUM_INPUTS, TIME_SLOTS, pData->inputframeTF); + /* Rotate source directions */ if (enableRotation && pData->recalc_M_rotFLAG) { yawPitchRoll2Rzyx (pData->yaw, pData->pitch, pData->roll, pData->useRollPitchYawFlag, Rxyz); - for(i=0; isrc_dirs_xyz[i][0] = cosf(DEG2RAD(pData->src_dirs_deg[i][1])) * cosf(DEG2RAD(pData->src_dirs_deg[i][0])); pData->src_dirs_xyz[i][1] = cosf(DEG2RAD(pData->src_dirs_deg[i][1])) * sinf(DEG2RAD(pData->src_dirs_deg[i][0])); pData->src_dirs_xyz[i][2] = sinf(DEG2RAD(pData->src_dirs_deg[i][1])); @@ -265,91 +284,92 @@ void binauraliserNF_process pData->recalc_M_rotFLAG = 0; } - /* Load time-domain data */ - for (i = 0; i < SAF_MIN(nSources, nInputs); i++) - utility_svvcopy(inputs[i], BINAURALISER_FRAME_SIZE, pData->inputFrameTD[i]); - for (; i < nSources; i++) - memset(pData->inputFrameTD[i], 0, BINAURALISER_FRAME_SIZE * sizeof(float)); - - /* Apply source gains */ - for (ch = 0; ch < nSources; ch++) { - if(fabsf(pData->src_gains[ch] - 1.f) > 1e-6f) - utility_svsmul(pData->inputFrameTD[ch], &(pData->src_gains[ch]), BINAURALISER_FRAME_SIZE, NULL); - } + /* Interpolate and apply HRTFs, apply DVF magnitude filter */ + /* Zero out TF summing bus */ + memset(FLATTEN3D(pData->outputframeTF), 0, HYBRID_BANDS*NUM_EARS*TIME_SLOTS * sizeof(float_complex)); - /* Zero out busses - not needed; overwritten, not summeed over */ - // memset(FLATTEN2D(pData->binsrcsTD), 0, NUM_EARS * BINAURALISER_FRAME_SIZE * sizeof(float)); - // memset(FLATTEN3D(pData->binauralTF), 0, HYBRID_BANDS * nSources * NUM_EARS * TIME_SLOTS * sizeof(float_complex)); - - /* Apply time-frequency transform (TFT) */ - afSTFT_forward_knownDimensions(pData->hSTFT, pData->inputFrameTD, BINAURALISER_FRAME_SIZE, MAX_NUM_INPUTS, TIME_SLOTS, pData->inputframeTF); - - /* Interpolate and apply HRTFs, maintain individual binaural source streams */ for (ch = 0; ch < nSources; ch++) { - /* Interpolate hrtfs */ + /* Interpolate HRTFs */ if (pData->recalc_hrtf_interpFLAG[ch]) { - if (enableRotation) - binauraliser_interpHRTFs(hBin, pData->interpMode, pData->src_dirs_rot_deg[ch][0], pData->src_dirs_rot_deg[ch][1], pData->hrtf_interp[ch]); - else - binauraliser_interpHRTFs(hBin, pData->interpMode, pData->src_dirs_deg[ch][0], pData->src_dirs_deg[ch][1], pData->hrtf_interp[ch]); + if (enableRotation) { + pData->src_dirs_cur = pData->src_dirs_rot_deg; + } else { + pData->src_dirs_cur = pData->src_dirs_deg; + } + binauraliser_interpHRTFs(hBin, pData->interpMode, pData->src_dirs_cur[ch][0], pData->src_dirs_cur[ch][1], pData->hrtf_interp[ch]); pData->recalc_hrtf_interpFLAG[ch] = 0; + pData->recalc_dvfCoeffFLAG[ch] = 1; + } + /* Update DVF filters with change in direction and/or distance */ + if (pData->recalc_dvfCoeffFLAG[ch]) { + rho = pData->src_dists_m[ch] * headRadiusRecip; + doaToIpsiInteraural(pData->src_dirs_cur[ch][0], pData->src_dirs_cur[ch][1], &alphaLR[0], NULL); + calcDVFCoeffs(alphaLR[0], rho, fs, pData->b_dvf[ch][0], pData->a_dvf[ch][0]); + calcDVFCoeffs(alphaLR[1], rho, fs, pData->b_dvf[ch][1], pData->a_dvf[ch][1]); + pData->recalc_dvfCoeffFLAG[ch] = 0; + + evalIIRTransferFunctionf(pData->b_dvf[ch][0], pData->a_dvf[ch][0], 2, + pData->freqVector, HYBRID_BANDS, fs, 0, + pData->dvfmags[ch][0], + pData->dvfphases[ch][0]); // NULL to return magnitude only + + evalIIRTransferFunctionf(pData->b_dvf[ch][1], pData->a_dvf[ch][1], 2, + pData->freqVector, HYBRID_BANDS, fs, 0, + pData->dvfmags[ch][1], + pData->dvfphases[ch][1]); // NULL to return magnitude only + pData->recalc_dvfCoeffFLAG[ch] = 0; } - /* Expand each source channel to binaural by applying HRTF filter (FD) */ - for (band = 0; band < HYBRID_BANDS; band++) { - for (ear = 0; ear < NUM_EARS; ear++) { - utility_cvsmul(pData->inputframeTF[band][ch], - &pData->hrtf_interp[ch][band][ear], TIME_SLOTS, - &pData->binauralTF[band][NUM_EARS * ch + ear][0]); + /* Convolve this channel with the interpolated HRTF, and add it to the binaural buffer */ + if (pData->src_dists_m[ch] < ffThresh) { + /* Near field: convolve this channel with the HRTF & DVF filter */ + float_complex dvfScale; + for (band = 0; band < HYBRID_BANDS; band++) { + for (ear = 0; ear < NUM_EARS; ear++) { + /* combine mag and phase response of HRTF and DVF */ + /* apply magnitude & phase */ + dvfScale = ccmulf(cmplxf(pData->dvfmags[ch][ear][band], pData->dvfphases[ch][ear][band]), pData->hrtf_interp[ch][band][ear]); + /* apply magnitude only, no phase */ + // dvfScale = crmulf(pData->hrtf_interp[ch][band][ear], pData->dvfmags[ch][ear][band]); + /* bypass dvf */ + // dvfScale = pData->hrtf_interp[ch][band][ear]; + + cblas_caxpy(TIME_SLOTS, &dvfScale, + pData->inputframeTF[band][ch], 1, + pData->outputframeTF[band][ear], 1); + } } + } else { + /* Far field: convolve this channel with the HRTF filter only */ + for (band = 0; band < HYBRID_BANDS; band++) + for (ear = 0; ear < NUM_EARS; ear++) + cblas_caxpy(TIME_SLOTS, &pData->hrtf_interp[ch][band][ear], + pData->inputframeTF[band][ch], 1, + pData->outputframeTF[band][ear], 1); } } - /* Bring sources back to time domain. - * Note, in the case of AFSTFT_TIME_CH_BANDS, args dataFD_nCH & dataFD_nHops are unused */ - afSTFT_backward_knownDimensions(pData->hSTFT, pData->binauralTF, BINAURALISER_FRAME_SIZE, MAX_NUM_INPUTS * NUM_EARS, TIME_SLOTS, pData->binsrcsTD); + /* scale by number of sources */ + cblas_sscal(/*re+im*/2*HYBRID_BANDS*NUM_EARS*TIME_SLOTS, 1.0f/sqrtf((float)nSources), (float*)FLATTEN3D(pData->outputframeTF), 1); - /* Iterate over sources, applying DVF to those in the near field */ - for (ch = 0; ch < nSources; ch++) { - int l = NUM_EARS * ch; - int r = l + 1; - if (src_dists[ch] < ffThresh) { - /* Get previous frame's last sample for DVF IIR filter */ - wzL = pData->binsrcsTD[l][BINAURALISER_FRAME_SIZE - 1]; - wzR = pData->binsrcsTD[r][BINAURALISER_FRAME_SIZE - 1]; - - rho = src_dists[ch] * headRadiusRecip; - convertFrontalDoAToIpsilateral(src_dirs[ch][0], &thetaLR[0]); - applyDVF(thetaLR[0], rho, pData->binsrcsTD[l], BINAURALISER_FRAME_SIZE, fs, &wzL, pData->binsrcsTD[l]); - applyDVF(thetaLR[1], rho, pData->binsrcsTD[r], BINAURALISER_FRAME_SIZE, fs, &wzR, pData->binsrcsTD[r]); - } - } + /* inverse-TFT */ + afSTFT_backward_knownDimensions(pData->hSTFT, pData->outputframeTF, BINAURALISER_FRAME_SIZE, NUM_EARS, TIME_SLOTS, pData->outframeTD); - /* Zero out all plugin output channels, cblas_saxpy below accumulates signals into the first 2 */ - for (ch = 0; ch < nOutputs; ch++) - memset(outputs[ch], 0, BINAURALISER_FRAME_SIZE * sizeof(float)); - - /* Sum binaural sources */ - // Note, it appears SAF doesn't support in-place vector addition, e.g. ippsAdd_32f_I - for (srci = 0; srci < nSources; srci++) { - for (ear = 0; ear < NUM_EARS; ear++) { - cblas_saxpy(BINAURALISER_FRAME_SIZE, 1.f, - pData->binsrcsTD[srci * NUM_EARS + ear], 1, outputs[ear], 1); - } - } - /* Scale binaural output by number of sources */ - for (ear = 0; ear < NUM_EARS; ear++) { - cblas_sscal(BINAURALISER_FRAME_SIZE, 1.0f / sqrtf((float)nSources), outputs[ear], 1); - } + /* Copy to output buffer */ + for (ch = 0; ch < SAF_MIN(NUM_EARS, nOutputs); ch++) + utility_svvcopy(pData->outframeTD[ch], BINAURALISER_FRAME_SIZE, outputs[ch]); + for (; ch < nOutputs; ch++) + memset(outputs[ch], 0, BINAURALISER_FRAME_SIZE*sizeof(float)); } - else { - for (ch = 0; ch < nOutputs; ch++) - memset(outputs[ch], 0, BINAURALISER_FRAME_SIZE * sizeof(float)); + else{ + for (ch=0; ch < nOutputs; ch++) + memset(outputs[ch],0, BINAURALISER_FRAME_SIZE*sizeof(float)); } pData->procStatus = PROC_STATUS_NOT_ONGOING; } + /* Set Functions */ void binauraliserNF_setSourceDist_m(void* const hBin, int index, float newDist_m) @@ -358,6 +378,7 @@ void binauraliserNF_setSourceDist_m(void* const hBin, int index, float newDist_m newDist_m = SAF_MAX(newDist_m, pData->nearfield_limit_m); if(pData->src_dists_m[index] != newDist_m){ pData->src_dists_m[index] = newDist_m; + pData->recalc_dvfCoeffFLAG[index] = 1; } } @@ -372,8 +393,10 @@ void binauraliserNF_setInputConfigPreset(void* const hBin, int newPresetID) if(pData->nSources != pData->new_nSources) binauraliser_setCodecStatus(hBin, CODEC_STATUS_NOT_INITIALISED); - for(ch=0; chrecalc_hrtf_interpFLAG[ch] = 1; + pData->recalc_dvfCoeffFLAG[ch] = 1; + } } float binauraliserNF_getSourceDist_m(void* const hBin, int index) @@ -397,7 +420,6 @@ float binauraliserNF_getFarfieldHeadroom(void* const hBin) return pData->farfield_headroom; } - float binauraliserNF_getNearfieldLimit_m(void* const hBin) { binauraliserNF_data *pData = (binauraliserNF_data*)(hBin); diff --git a/examples/src/binauraliser_nf/binauraliser_nf_internal.h b/examples/src/binauraliser_nf/binauraliser_nf_internal.h index 9c405f44b..5fece8606 100644 --- a/examples/src/binauraliser_nf/binauraliser_nf_internal.h +++ b/examples/src/binauraliser_nf/binauraliser_nf_internal.h @@ -56,12 +56,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 */ @@ -70,7 +71,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 */ @@ -81,19 +82,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] */ @@ -102,7 +103,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 */ @@ -125,23 +126,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 */ /* ========================================================================== */ diff --git a/framework/modules/saf_utilities/saf_utility_dvf.c b/framework/modules/saf_utilities/saf_utility_dvf.c index 11d5427e2..2c20ba9cb 100644 --- a/framework/modules/saf_utilities/saf_utility_dvf.c +++ b/framework/modules/saf_utilities/saf_utility_dvf.c @@ -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,13 +101,13 @@ 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, @@ -115,55 +115,51 @@ void interpHighShelfParams ) { 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; + } } diff --git a/framework/modules/saf_utilities/saf_utility_dvf.h b/framework/modules/saf_utilities/saf_utility_dvf.h index 2391132c6..81d6dba9d 100644 --- a/framework/modules/saf_utilities/saf_utility_dvf.h +++ b/framework/modules/saf_utilities/saf_utility_dvf.h @@ -41,65 +41,65 @@ extern "C" { #include /** - * Apply the Distance Variation function to the input signal, as described in [1]. + * Calculate the Distance Variation Function (DVF) filter coefficients, + * as described in [1]. * * @see [1] S. Spagnol, E. Tavazzi, and F. Avanzini, “Distance rendering and * perception of nearby virtual sound sources with a near-field filter * model,” Applied Acoustics, vol. 115, pp. 61–73, Jan. 2017, * doi: 10.1016/j.apacoust.2016.08.015. * - * @param[in] theta Ipsilateral azimuth, on the inter-aural axis [0..180] - * * (degrees) - * @param[in] rho Source distance, normalized to head radius, >= 1 - * @param[in] in_signal (&) Input signal pointer. - * @param[in] nSamples Number of samples to process - * @param[in] fs Sample rate - * @param[out] wz (&) Filter coefficients to be passed to the next block - * @param[out] out_signal (&) Output signal pointer + * @param[in] alpha Lateral angle, similar the interaural-polar convention, + * but specified as an offset from the interaural axis, + * [0, 180] (deg). See `doaToIpsiInteraural()` + * to convert frontal azimuth/elevation to the expected + * format. + * @param[in] rho Source distance, normalized to head radius, >= 1. + * @param[in] fs Sample rate. + * @param[out] b Numerator coefficients for the DVF shelving filter. + * @param[out] a Denominator coefficients for the DVF shelving filter */ -void applyDVF(/* Input Arguments */ - float theta, - float rho, - float* in_signal, - int nSamples, - float fs, - /* Output Arguments */ - float* wz, - float* out_signal); +void calcDVFCoeffs( + /* Input Arguments */ + float alpha, + float rho, + float fs, + /* Output Arguments */ + float * b, + float * a); /** - * Apply the Distance Variation function to the input signal + * Calculate the shelving filter parameters for the Distance Variation Function + * filter from the source (ipsilateral) azimuth and distance. * - * @test test__dvf_calcHighShelfParams(), test__dvf_interpHighShelfParams() - * - * @param[in] theta Ipsilateral azimuth, on the inter-aural axis [0..180] (deg) + * @param[in] theta Lateral angle, on the inter-aural axis [0..180] (deg) * @param[in] rho Source distance, normalized to head radius, >= 1 - * @param[out] iG0 (&) interpolated DC gain - * @param[out] iGInf (&) interpolated high shelf gain - * @param[out] iFc (&) interpolated high shelf cutoff frequency + * @param[out] iG0 Interpolated DC gain + * @param[out] iGInf Interpolated high shelf gain + * @param[out] iFc Interpolated high shelf cutoff frequency */ -void interpHighShelfParams(/* Input Arguments */ - float theta, - float rho, - /* Output Arguments */ - float* iG0, - float* iGInf, - float* iFc); + +void interpDVFShelfParams( + /* Input Arguments */ + float theta, + float rho, + /* Output Arguments */ + float* iG0, + float* iGInf, + float* iFc); /** - * Apply the Distance Variation function to the input signal - * - * @test test__dvf_calcIIRCoeffs() + * Calculate the DVF filter coefficients from shelving filter parameters. * - * @param[in] g0 High shelf gain at DC [dB] - * @param[in] gInf High shelf gain at Nyquist frequency [dB] - * @param[in] fc Shelf cutoff frequency [Hz] - * @param[out] fs Sample rate - * @param[out] b0 (&) IIR numerator coefficient 1 - * @param[out] b1 (&) IIR numerator coefficient 2 - * @param[out] a1 (&) IIR denominator coefficient 2 + * @param[in] g0 High shelf gain at DC [dB]. + * @param[in] gInf High shelf gain at Nyquist frequency [dB]. + * @param[in] fc Shelf cutoff frequency [Hz]. + * @param[out] fs Sample rate. + * @param[out] b0 Numerator coefficient 1. + * @param[out] b1 Numerator coefficient 2. + * @param[out] a1 Denominator coefficient 2. */ -void calcIIRCoeffs(/* Input Arguments */ +void dvfShelfCoeffs(/* Input Arguments */ float g0, float gInf, float fc, @@ -110,36 +110,46 @@ void calcIIRCoeffs(/* Input Arguments */ float* a1); /** - * Calculate the High shelf gains and cutoff parameters, given a azimuth index - * 'i' and distance 'rho'. This will be called twice per parameter change and + * Calculate the high shelf gains and cutoff parameters, given a azimuth index + * `i` and distance `rho`. This will be called twice per parameter change and * the shelf filter parameters will be linearly interpolated according to the - * precise azimuth. + * azimuth. * - * @param[in] i Coefficient table row index () + * @param[in] i Coefficient table row index * @param[in] rho Normalized source distance (1 = head radius) * @param[out] g0 High shelf gain at DC [dB] * @param[out] gInf High shelf gain at Nyquist frequency [dB] * @param[out] fc Shelf cutoff frequency [Hz] */ -void calcHighShelfParams(int i, - float rho, - float* g0, - float* gInf, - float* fc); +void calcDVFShelfParams(int i, + float rho, + float* g0, + float* gInf, + float* fc); /** - * Convert user-supplied DoA theta (0˚ forward) to a theta value understood by - * the filter designer (DoA on the interaural axis, 0˚ is ipsilateral). Note - * because it's based on a spherical head, the sign is always positive. Note - * input is expected to be [-180, 180] + * Convert a frontal azimuth/elevation to a modified Interaural-Polar + * coordinate. Whereas Interaural-Polar coordinates are with reference to the + * median plane, alpha [0, 90], beta [0, 180] this modification is with + * reference to the transverse plane (ipsilateral ear direction), alpha + * [0, 180], beta [0, 90]. This is intended for the input to + * `interpDVFShelfParams()` for calculating DVF filter parameters, which + * are framed as an offset from the interaural axis, and based on a spherical + * head model (i.e. elevation translates to a change in lateral angle). * - * @param[in] thetaFront Source DoA, 0˚ is forward-facing, positive angles move - * counter-clockwise [deg, (-180, 180)] - * @param[out] ipsiDoaLR High shelf gain at DC [dB] + * @param[in] azimuth Source DoA, 0˚ is forward-facing, angle increases + * counter-clockwise (deg, [-360, 360]). + * @param[in] elevation Source elevation, angles increase upward from the + * horizon (deg, [-180, 180]). + * @param[out] alphaLR 2-element array of lateral angle alpha for left and + * right ear (deg, [0,180]]). + * @param[out] betaLR 2-element array of vertal angle beta for left and + * right ear (deg, [0,90]]). */ -void convertFrontalDoAToIpsilateral(float thetaFront, - float* ipsiDoaLR); - +void doaToIpsiInteraural(float azimuth, + float elevation, + float* alphaLR, + float* betaLR); #ifdef __cplusplus }/* extern "C" */ diff --git a/framework/modules/saf_utilities/saf_utility_filters.c b/framework/modules/saf_utilities/saf_utility_filters.c index cd47fdffb..55c0b004a 100644 --- a/framework/modules/saf_utilities/saf_utility_filters.c +++ b/framework/modules/saf_utilities/saf_utility_filters.c @@ -592,7 +592,7 @@ void evalBiQuadTransferFunction /* substituting Euler, z = e^(-jwn) = cos(wn) + j*sin(wn), into: * H(z) = (b0 + b1*z^(-1) + b2*z^(-2)) / (1 + a1*z^(-1) + a2*z^(-2)): */ denom_real = 1.0f + a[1]*cosf(w) + a[2]*cosf(2.0f*w); - denom_imag = a[1]*sinf(w) + a[2]*sinf(2*w); + denom_imag = a[1]*sinf(w) + a[2]*sinf(2.0f*w); num_real = b[0] + b[1]*cosf(w) + b[2]*cosf(2.0f*w); num_imag = b[1]*sinf(w) + b[2]*sinf(2.0f*w); @@ -606,6 +606,136 @@ void evalBiQuadTransferFunction } } +void evalIIRTransferFunctionf + ( + float* b_coeff, /* Note coeffs are *floats* */ + float* a_coeff, + int nCoeffs, + float* freqs, + int nFreqs, + float fs, + int mag2dB, + float* magnitude, + float* phase_rad +) +{ + int ff; + float w, x, a, b, c, d, h_re, h_im, cosx, sinx; + float norm_frq = -2.0 * SAF_PI / fs; // -1 factored in from negated 'n' below + + /* + * substitute Euler's z = e^(jwn) = cos(wn) + j*sin(wn) + * into + * [b0*z^(0) + b1*z^(-1) ... + bn*z^(-n)] + * H(z) = -------------------------------------- + * [a0*z^(0) + a1*z^(-1) ... + an*z^(-n)] + */ + + for(ff = 0; ff < nFreqs; ff++){ + w = freqs[ff] * norm_frq; + + // n = 0; real terms are b[0] and a[0], imag terms are 0 + a = b_coeff[0]; // num_re; b_coeff[0] * cos(x); x = -1.f * n * w; + b = 0.0; // num_imag; b_coeff[0] * sin(x); + c = a_coeff[0]; // den_re + d = 0.0; // den_imag + + /* sum over remaining numerator and denominator terms */ + for(int n = 1; n < nCoeffs; n++){ + x = (n * w); // n * -1.f * w (-1.f applied in norm_frq) + cosx = cosf(x); // alt: cosx = 1 - 2.f * powf(sinf(x/2.f), 2.f); + sinx = sinf(x); + a += b_coeff[n] * cosx; // 'a' + b += b_coeff[n] * sinx; // 'b' + c += a_coeff[n] * cosx; // 'c' + d += a_coeff[n] * sinx; // 'd' + } + + /* 1 / (c^2 + d^2 + eps) */ + double dvsr = 1.0 / (powf(c, 2.f) + powf(d, 2.f) + 2.23e-7f); + + if(magnitude!=NULL){ + /* sqrt((a^2 + b^2) / (c^2 + d^2)) */ + magnitude[ff] = (float)sqrt( (powf(a, 2.0f) + powf(b, 2.0f)) * dvsr ); + if(mag2dB) + magnitude[ff] = 20.0f*log10f(magnitude[ff]); + } + + if(phase_rad!=NULL) { + /* complex division */ + h_re = (a*c + b*d) * dvsr; + h_im = (b*c - a*d) * dvsr; + phase_rad[ff] = (float)atan2(h_im, h_re); + } + } +} + +void evalIIRTransferFunction + ( + double* b_coeff, /* Note coeffs are *doubles* */ + double* a_coeff, + int nCoeffs, + float* freqs, + int nFreqs, + float fs, + int mag2dB, + float* magnitude, + float* phase_rad +) +{ + int ff; + float w; + double x, a, b, c, d, h_re, h_im, cosx, sinx; + float norm_frq = -2.0 * SAF_PI / fs; // -1 factored in from negated 'n' below + + /* + * substitute Euler's z = e^(jwn) = cos(wn) + j*sin(wn) + * into + * [b0*z^(0) + b1*z^(-1) ... + bn*z^(-n)] + * H(z) = -------------------------------------- + * [a0*z^(0) + a1*z^(-1) ... + an*z^(-n)] + */ + + for(ff = 0; ff < nFreqs; ff++){ + w = freqs[ff] * norm_frq; + + // n = 0; real terms are b[0] and a[0], imag terms are 0 + a = b_coeff[0]; // num_re; b_coeff[0] * cos(x); x = -1.f * n * w; + b = 0.0; // num_imag; b_coeff[0] * sin(x); + c = a_coeff[0]; // den_re + d = 0.0; // den_imag + + /* sum over remaining numerator and denominator terms */ + for(int n = 1; n < nCoeffs; n++){ + x = (double)(n * w); // n * -1.f * w (-1.f applied in norm_frq) + cosx = 1 - 2.f * pow(sin(x/2.f), 2.f); // cos(x) = 1 - 2sin^2(x/2) + sinx = sin(x); + a += b_coeff[n] * cosx; // 'a' + b += b_coeff[n] * sinx; // 'b' + c += a_coeff[n] * cosx; // 'c' + d += a_coeff[n] * sinx; // 'd' + } + + /* 1 / (c^2 + d^2 + eps) */ + double dvsr = 1.0 / (pow(c, 2.f) + pow(d, 2.f) + 2.23e-17f); + + if(magnitude!=NULL){ + /* sqrt((a^2 + b^2) / (c^2 + d^2)) */ + magnitude[ff] = (float)sqrt( (pow(a, 2.0f) + pow(b, 2.0f)) * dvsr ); + if(mag2dB) + magnitude[ff] = 20.0f*log10f(magnitude[ff]); + } + + if(phase_rad!=NULL) { + /* complex division */ + h_re = (a*c + b*d) * dvsr; + h_im = (b*c - a*d) * dvsr; + phase_rad[ff] = (float)atan2(h_im, h_re); + } + } +} + + void applyIIR ( float* in_signal, diff --git a/framework/modules/saf_utilities/saf_utility_filters.h b/framework/modules/saf_utilities/saf_utility_filters.h index cf0c28830..1fc875368 100644 --- a/framework/modules/saf_utilities/saf_utility_filters.h +++ b/framework/modules/saf_utilities/saf_utility_filters.h @@ -271,6 +271,84 @@ void evalBiQuadTransferFunction(/* Input arguments */ float* magnitude, float* phase_rad); +/** + * Computes magnitude and phase response of an IIR filter from its coefficients + * at user-specified frequencies (Hz). + * + * The function optionally returns magnitude and/or phase. The function is + * tested against Matlab's 'freqz' function. + * + * @note This function operates on filter coefficients of type **double**, for + * more accuracy for higher order filters, in particular at lower frequencies. + * See `evalIIRTransferFunctionf()` (operating on **float** type coefficients) + * for a more efficient implementation (e.g. for lower order filters). + * + * @test test__evalIIRTransferFunction() + * + * @param[in] b Filter coefficients for the numerator; nCoeffs x 1 + * @param[in] a Filter coefficients for the denominator; nCoeffs x 1 + * @param[in] nCoeffs Number of filter coefficients + * @param[in] freqs Frequencies at which to evaluate the + * magnitude/phase response (Hz); nFreqs x 1 + * @param[in] nFreqs Number of frequencies + * @param[in] fs Sampling rate (Hz) + * @param[in] mag2dB 0: 'magnitude' returned in linear scale, 1: dB scale + * @param[out] magnitude Magnitudes at each frequency (set to NULL of not + * wanted); nFreqs x 1 + * @param[out] phase_rad Phases at each frequency (radians; set to NULL of + * not wanted); nFreqs x 1 + */ +void evalIIRTransferFunction(/* Input arguments */ + double* b, + double* a, + int nCoeffs, + float* freqs, + int nFreqs, + float fs, + int mag2dB, + /* Output arguments */ + float* magnitude, + float* phase_rad); + +/** + * Computes magnitude and phase response of an IIR filter from its coefficients + * (floats) at user-specified frequencies (Hz). + * + * @note This function operates on filter coefficients of type **float**, and is + * suitable for low order filters, or otherwise less accuracy at low frequencies + * for higher order filters. For higher-order filters, especially when evaluated + * at low frequencies, `evalIIRTransferFunction()` is recommended. + * + * The function optionally returns magnitude and/or phase. The function is + * tested against Matlab's 'freqz' function. + * + * @test test__evalIIRTransferFunction() + * + * @param[in] b Filter coefficients for the numerator; nCoeffs x 1 + * @param[in] a Filter coefficients for the denominator; nCoeffs x 1 + * @param[in] nCoeffs Number of filter coefficients + * @param[in] freqs Frequencies at which to evaluate the + * magnitude/phase response (Hz); nFreqs x 1 + * @param[in] nFreqs Number of frequencies + * @param[in] fs Sampling rate (Hz) + * @param[in] mag2dB 0: 'magnitude' returned in linear scale, 1: dB scale + * @param[out] magnitude Magnitudes at each frequency (set to NULL of not + * wanted); nFreqs x 1 + * @param[out] phase_rad Phases at each frequency (radians; set to NULL of + * not wanted); nFreqs x 1 + */ +void evalIIRTransferFunctionf(/* Input arguments */ + float* b, + float* a, + int nCoeffs, + float* freqs, + int nFreqs, + float fs, + int mag2dB, + /* Output arguments */ + float* magnitude, + float* phase_rad); + /** * Applies an IIR filter to a time-domain signal (using the direct form II * difference equation) diff --git a/test/include/saf_test.h b/test/include/saf_test.h index 64ddeac0a..8f2b2db81 100644 --- a/test/include/saf_test.h +++ b/test/include/saf_test.h @@ -158,6 +158,10 @@ void test__latticeDecorrelator(void); * Testing that the coefficients computed with butterCoeffs() are numerically * similar to the "butter" function in Matlab */ void test__butterCoeffs(void); +/** + * Testing that the magnitudes and phases returned by evalIIRTransferFunction() + * are numerically similar to the "freqz" function in Matlab */ +void test__evalIIRTransferFunction(void); /** * Testing that the faf_IIRFilterbank can reconstruct the original signal power */ @@ -169,15 +173,15 @@ void test__gexpm(void); /** * Calculate high shelf parameters, g0, gInf, fc, from the lookup table * coefficients (10 degree steps) */ -void test__dvf_calcHighShelfParams(void); +void test__dvf_calcDVFShelfParams(void); /** * Test the interpolation of high shelf parameters based on distance and * incidence angle parameters */ -void test__dvf_interpHighShelfParams(void); +void test__dvf_interpDVFShelfParams(void); /** - * Test the generation of high shelf coeffs based on shelf gains and fc + * Test the generation of high shelf coeffs based on shelf gain and fc * parameters */ -void test__dvf_calcIIRCoeffs(void); +void test__dvf_dvfShelfCoeffs(void); /* ========================================================================== */ diff --git a/test/src/saf_test.c b/test/src/saf_test.c index c69d3051f..570bd5521 100644 --- a/test/src/saf_test.c +++ b/test/src/saf_test.c @@ -105,11 +105,12 @@ int main_test(void) { RUN_TEST(test__unique_i); RUN_TEST(test__latticeDecorrelator); RUN_TEST(test__butterCoeffs); + RUN_TEST(test__evalIIRTransferFunction); RUN_TEST(test__faf_IIRFilterbank); RUN_TEST(test__gexpm); - RUN_TEST(test__dvf_calcHighShelfParams); - RUN_TEST(test__dvf_interpHighShelfParams); - RUN_TEST(test__dvf_calcIIRCoeffs); + RUN_TEST(test__dvf_calcDVFShelfParams); + RUN_TEST(test__dvf_interpDVFShelfParams); + RUN_TEST(test__dvf_dvfShelfCoeffs); /* SAF cdf4sap module unit tests */ RUN_TEST(test__formulate_M_and_Cr); @@ -172,7 +173,7 @@ int main_test(void) { RUN_TEST(test__saf_example_ambi_bin); RUN_TEST(test__saf_example_ambi_dec); RUN_TEST(test__saf_example_ambi_enc); - RUN_TEST(test__saf_example_array2sh); + RUN_TEST(test__saf_example_array2sh); RUN_TEST(test__saf_example_rotator); RUN_TEST(test__saf_example_spreader); #endif /* SAF_ENABLE_EXAMPLES_TESTS */ diff --git a/test/src/test__utilities_module.c b/test/src/test__utilities_module.c index a2b103a8b..e1fd62993 100644 --- a/test/src/test__utilities_module.c +++ b/test/src/test__utilities_module.c @@ -982,6 +982,216 @@ void test__butterCoeffs(void){ } } +/* Test evalIIRTransferFunction() + * Coefficients for the first 7 tests below are taken from the butterworth tests + * above, so can be compared against the mag/phase results given by MATLAB's + * freqz function. The 8th test loop evaluates results for 1st order shelving + * filters with coefficients generated by the DVF filter functions. The last + * test loop runs those same tests, but with the floating point version of + * evalIIRTransferFunctionf() (valid for low order filters) + */ +void test__evalIIRTransferFunction(void) { + int i, nCoeffs; + float phase_ref, mag_ref, mag_ref_db; + + /* Config */ + const double phaseTolerance = 0.0174533f * 5.f; /* ~ 1 degree * mul */ + const double magToleranceDb = 0.1f; /* tolerance in dB, for a target magnitude of 0dB */ + const float errScale = 2.f / 120.f; /* tolerance grows for lower dB target: toleranceLevel/atLevel. + * e.g. 2/120 = 2dB tolerance for -120 target dB */ + const int nFreqs = 10; + const float freqs[10] = {147.21423,270.49564,411.40091,687.90202,1395.3072,2024.3936,3696.9416,6784.4745,9798.67,17594.058}; + float mag[10], phase[10]; + float fs = 44.1e3f; + + /* evalIIRTransferFunction(): coeffs of type double */ + + /* Test 1 * 1st order Low-pass filter */ + nCoeffs = 2; + const double a_t1[2] = {1,-0.6681786}; + const double b_t1[2] = {0.1659107,0.1659107}; + const double mag_ref1[10] = {0.99861294,0.99533929,0.98931332,0.97092312,0.89393904,0.80765661,0.59366549,0.35440473,0.23070415,0.06521195}; + const double phase_ref1[10] = {-0.052676048,-0.096585083,-0.14632679,-0.24173907,-0.46473796,-0.63062928,-0.93519006,-1.2085189,-1.337995,-1.5055381}; + evalIIRTransferFunction((double*)b_t1, (double*)a_t1, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); + for(i = 0; i < nFreqs; i++) { + phase_ref = phase_ref1[i]; + mag_ref = mag_ref1[i]; + mag_ref_db = 20*log10(mag_ref); + TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i])); + TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]); + } + + /* Test 2 * 2nd order Low-pass filter */ + nCoeffs = 3; + const double a_t2[3] = {1,-0,0.1715729}; + const double b_t2[3] = {0.2928932,0.5857864,0.2928932}; + const float mag_ref2[10] = {0.99999991,0.99999985,0.99999955,0.99999702,0.99995046,0.99977761,0.99736787,0.96409579,0.81776268,0.1073164}; + const float phase_ref2[10] = {-0.014832279,-0.027258003,-0.041470589,-0.069414188,-0.14150066,-0.20679977,-0.39012494,-0.7974393,-1.3261562,-2.6614056}; + evalIIRTransferFunction((double*)b_t2, (double*)a_t2, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); + for(i = 0; i < nFreqs; i++) { + phase_ref = phase_ref2[i]; + mag_ref = mag_ref2[i]; + mag_ref_db = 20*log10(mag_ref); + TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i])); + TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]); + } + + /* Test 3 * 3rd order Low-pass filter */ + nCoeffs = 4; + const double a_t3[4] = {1,-2.9476416,2.896645,-0.9489859}; + const double b_t3[4] = {2.2e-06,6.6e-06,6.6e-06,2.2e-06}; + const float mag_ref3[10] = {0.8954383,0.3011618,0.0892913,0.0191409,0.0022769,0.0007374,0.0001152,1.56e-05,3.8e-06,1e-07}; + const float phase_ref3[10] = {-1.8249735,3.0678618,2.4995092,2.1114704,1.8340934,1.751328,1.6679375,1.6206872,1.6020054,1.579398}; + evalIIRTransferFunction((double*)b_t3, (double*)a_t3, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); + for(i = 0; i < nFreqs; i++) { + phase_ref = phase_ref3[i]; + mag_ref = mag_ref3[i]; + mag_ref_db = 20*log10(mag_ref); + TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i])); + TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]); + } + + /* Test 4 * 6th order Low-pass filter */ + nCoeffs = 7; + const double a_t4[7] = {1,-5.49431292177096,12.5978414666894,-15.4285267903275,10.6436770055305,-3.92144696766748,0.6027721469713}; + const double b_t4[7] = {6.15535184628202e-08,3.69321110776921e-07,9.23302776942303e-07,1.2310703692564e-06,9.23302776942303e-07,3.69321110776921e-07,6.15535184628202e-08}; + const float mag_ref4[10] = {0.9999999834907868,0.9999997831836054,0.9999679556869572,0.9849426248859378,0.08033081621985783,0.008452216914022819,0.0002063542729228268,3.793812554381118e-06,2.274031694371124e-07,9.970589432354785e-11}; + const float phase_ref4[10] = {-0.6201852189230334,-1.148525513374147,-1.774695369143539,3.109543344373707,-0.4296773811384472,-1.349824316530828,-2.195405632723407,-2.65814688739603,-2.839508904295157,-3.058387834019209}; + evalIIRTransferFunction((double*)b_t4, (double*)a_t4, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); + for(i = 0; i < nFreqs; i++) { + float phase_ref = phase_ref4[i]; + float mag_ref = mag_ref4[i]; + float mag_ref_db = 20*log10(mag_ref); + TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i])); + TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]); + } + + /* Test 5 * 3rd order High-pass filter */ + nCoeffs = 4; + const double a_t5[4] = {1,-2.2191686,1.7151178,-0.4535459}; + const double b_t5[4] = {0.673479,-2.0204371,2.0204371,-0.673479}; + const float mag_ref5[10] = {0.0001466,0.0009096,0.0032014,0.0149875,0.125037,0.362653,0.927991,0.9985214,0.9999112,0.9999999}; + const float phase_ref5[10] = {-1.6762949,-1.7648759,-1.866651,-2.0692621,-2.6256366,3.0800183,1.6530258,0.7789431,0.4789307,0.1307956}; + evalIIRTransferFunction((double*)b_t5, (double*)a_t5, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); + for(i = 0; i < nFreqs; i++) { + float phase_ref = phase_ref5[i]; + float mag_ref = mag_ref5[i]; + float mag_ref_db = 20*log10(mag_ref); + TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i])); + TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]); + } + + /* Test 6 * 4th order High-pass filter + * 400 Hz cut (differs from butterworth test above) */ + nCoeffs = 5; + const double a_t6[5] = {1,-3.863184622426,5.598835456747838,-3.607752453919942,0.872108645089876}; + const double b_t6[5] = {4.3909323578772e-07,1.75637294315089e-06,2.63455941472633e-06,1.75637294315089e-06,4.3909323578772e-07}; + const float mag_ref6[10] = {0.9996691528983467,0.9595570109649983,0.5370184819357747,0.08100263003740536,0.004753436194609436,0.001057169058757887,8.896712774518116e-05,6.197328265811134e-06,9.491865964914827e-07,5.478157027512644e-09}; + const float phase_ref6[10] = {-1.072517623166929,-2.13344694428915,2.732267641095127,1.462991201859678,0.6929733816699927,0.4733493046075806,0.2541184532330854,0.130425028023503,0.08157492611996242,0.02248140228360206}; + evalIIRTransferFunction((double*)b_t6, (double*)a_t6, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); + for(i = 0; i < nFreqs; i++) { + float phase_ref = phase_ref6[i]; + float mag_ref = mag_ref6[i]; + float mag_ref_db = 20*log10(mag_ref); + TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i])); + TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]); + } + + /* Test 7 * 2nd order Band-pass filter */ + nCoeffs = 5; + const double a_t7[5] = {1,-3.9431258,5.832267,-3.8351187,0.9459779}; + const double b_t7[5] = {0.0003751,0,-0.0007501,0,0.0003751}; + const float mag_ref7[10] = {0.7829909,0.9051549,0.5636772,0.1816557,0.0400635,0.0185759,0.0053305,0.0014022,0.0005484,4.16e-05}; + const float phase_ref7[10] = {0.4017825,-0.7852502,-1.8127451,-2.4983166,-2.8544848,-2.9475768,-3.0381483,-3.0886103,-3.1084696,-3.1324667}; + evalIIRTransferFunction((double*)b_t7, (double*)a_t7, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); for(i=0; i