Loudness according to Annex F of ECMA-74

[mglesser]

Thanks to @djcaminero work, we have a first draft of ECMA-74 annex F hearing model to compute loudness in MOSQITO. However the model doesn't pass the validation stage yet (see figure below with the loudness contour from ISO 226 in dotted line and the loudness contour from MOSQITO hearing model in solid lines).

Any help from the community to review @djcaminero code will be welcome. You can reproduce the figure above by running mosqito/validation/hearing_model/hearing_model_validation.py script. That's a good starting point for a review.

If you would like to contribute but are not familiar with github, please contact me. I will be happy to help.

[mglesser]

I'm begining a step by step validation of the implementation. Step 1 is the outer and middle/inner ear filtering (section F 3.2 of the ECMA). By running mosqito/validations/hearing_model/valid_ear_filter_design.py, i get a response compliant to the standard, see figure below (ECMA-74 on the left, MOSQITO on the right)

[mglesser]

Second step is the gammatone filters design. I have created a new 'gammatone' function from 'afb_coefficients' and part of 'afb_ppal_parameters' the objective is to be able to easily compare MOSQITO implementation of the gammatone filter with the one of scipy. The function is not used in the main code yet, but it should be eventually.

The filters from mosqito match the one from scipy enough to consider that it is, at least, not totally wrong considering that it is not necessarily the same auditory model that is implemented (run mosqito/validations/hearing_model/valid_gammatone.py to get the plots):

[mglesser]

I checked the validation script and realized that it was not giving a correct evaluation of the implementation. I've programmed an alternative validation script. The principle is the following : for each frequency, the Loudness is calculated for 4 tones with an SPL amplitude determined by the isophone curves 20, 40, 60 and 80. A satisfying implementation should return values in phon close to 20, 40, 60 and 80 respectively. The result of the current implementation is the following (theoretical isophon in dotted lines, results from MOSQITO in solid lines):

Update 13/08/2021
I've programmed an alternative validation script. The principle is the following : for each frequency, the Loudness is calculated for 4 tones with an SPL amplitude determined by the isophone curves 20, 40, 60 and 80. A satisfying implementation should return values in phon close to 20, 40, 60 and 80 for each frequency (dotted lines in the figure below). The current implementation is in solid lines.

[kmh1201]

I changed variable centre_freq_array in the call of function annex_f_validation() in comp_loudness() to variable centre_freq for working diagnostics: (comp_loudness(signal, validation=True)).

More detailed explanation:
If the validation option of function comp_loudness(signal, validation=False) is set to True, function annex_f_validation() is called at the end with one of the arguments being centre_freq_array. In the current version, variable centre_freq_array is an empty list, thus, interrupting the validation process. The appropriate variable to finish the validation is called centre_freq:

[...]
if validation:
annex_f_validation(
fs,
signal,
signal_filtered,
sos_ear,
centre_freq,
f_bandwidth_array,
t_delay_array,
d_coefficients_array,
sb_array,
sh_array,
am_mod_coefficient_array,
bm_mod_coefficient_array,
band_pass_signal_array,
n_array,
t_array,
)

[mglesser]

I've been reviewing the whole code. I took this opportunity to implement an alternative version of the comp_loudness function. Some of the steps have been simplified to allow for easyer debugging. However, the results remains the same.

[mglesser]

I noticed that if I do not apply the Outer and middle/inner ear filtering (section 5.1.2 of the ECMA-418-2), I get a overall trend of the hearing model that is closer to the theoretical equal loudness contours. It suggests that the filter is not applied correctly to the signal.

[mglesser]

I've finally spotted two one errors in the "outer and midle ear filtering" section of the code.

1. The sosfiltfilt function was used to apply the filter.

signal_filtered = sp.signal.sosfiltfilt(sos_ear, signal, axis=0)

As stated in the doc : "This function applies a linear digital filter twice, once forward and once backwards. The combined filter has zero phase and a filter order twice that of the original." The filter order was then twice too high.

2. The filter predict the response of the ear. In order to "compensate" this response and go from physical amplitude in dB to physiological amplitude in sone, one should apply the inverse filter:

We are now getting dangerously close to the theoretical equal loudness contours. Next mistakes are to be search in a portion of the code that affect differently the low and high amplitudes. I would guess the "Non-linerarity" section.

[mglesser]

Solved by PR #33