Spaces:
Runtime error
Runtime error
| """Contains functions for generating and using equal-loudness contours for | |
| side-presented steady pure tones according to the ISO/IEC 226 standard. | |
| Code from https://gist.github.com/sammosummo/777debf946d0356acada and Cédric Colas | |
| """ | |
| __author__ = 'Sam Mathias' | |
| __version__ = 1.0 | |
| import numpy as np | |
| from scipy.interpolate import interp1d | |
| a_freq = 440 | |
| a_pitch = 69 | |
| upright_piano_dynamic_range = 60 # in db | |
| a_440_db_at_max_velocity = 75 | |
| a_440_amplitude_at_max_velocity = 10 ** (a_440_db_at_max_velocity / 20) | |
| def iso266(phon, return_freq=False): | |
| """Returns an equal-loudness contour evaluated at 29 frequencies between | |
| 20 Hz and 12.5 kHz according to the ISO/IEC 226 standard [1]_. | |
| Parameters | |
| ---------- | |
| phon : int or float | |
| The phon value represented by the equal-loudness contour, where a value | |
| of :math:`x` phon is the loudness of 1-KHz steady pure tone presented | |
| at :math:`x` dB SPL. Must be between 0 and 90. | |
| return_freq : bool, optional | |
| If True, the function returns the frequency values as well as the SPL | |
| values of the contour. Default is False. | |
| Returns | |
| ------- | |
| array_like | |
| Either a 1-D or a 2-D numpy array, depending on `return_freq`. | |
| Reference | |
| --------- | |
| .. [1] ISO/IEC (2003). ISO/IEC 226:2003 Acoustics -- Normal equal-loudness- | |
| level contours. | |
| http://www.iso.org/iso/catalogue_detail.htm?csnumber=34222. | |
| Example | |
| ------- | |
| elc = iso266(60, return_freq=True) | |
| print elc | |
| [[ 20. 25. 31.5 40. 50. | |
| 63. 80. 100. 125. 160. | |
| 200. 250. 315. 400. 500. | |
| 630. 800. 1000. 1250. 1600. | |
| 2000. 2500. 3150. 4000. 5000. | |
| 6300. 8000. 10000. 12500. ] | |
| [ 109.51132227 104.22789784 99.07786826 94.17731862 | |
| 89.96345731 85.94342131 82.05340072 78.65461863 | |
| 75.56345314 72.4743448 69.86431929 67.53483532 | |
| 65.39173983 63.45099627 62.0511792 60.81495942 | |
| 59.88668375 60.011588 62.1549143 63.18935604 | |
| 59.96161453 57.25515019 56.42385863 57.56993838 | |
| 60.8882125 66.36125056 71.66396598 73.15510401 | |
| 68.63077045]] | |
| """ | |
| if not 0 <= phon <= 90: | |
| raise ValueError('Cannot calculate for this value.') | |
| f = np.array([ | |
| 20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500, | |
| 630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000, 5000, 6300, 8000, | |
| 10000, 12500 | |
| ]) | |
| af = np.array([ | |
| 0.532, 0.506, 0.480, 0.455, 0.432, 0.409, 0.387, 0.367, 0.349, 0.330, | |
| 0.315, 0.301, 0.288, 0.276, 0.267, 0.259, 0.253, 0.250, 0.246, 0.244, | |
| 0.243, 0.243, 0.243, 0.242, 0.242, 0.245, 0.254, 0.271, 0.301 | |
| ]) | |
| Lu = np.array([ | |
| -31.6, -27.2, -23.0, -19.1, -15.9, -13.0, -10.3, -8.1, -6.2, -4.5, | |
| -3.1, -2.0, -1.1, -0.4, 0.0, 0.3, 0.5, 0.0, -2.7, -4.1, -1.0, 1.7, | |
| 2.5, 1.2, -2.1, -7.1, -11.2, -10.7, -3.1 | |
| ]) | |
| Tf = np.array([ | |
| 78.5, 68.7, 59.5, 51.1, 44.0, 37.5, 31.5, 26.5, 22.1, 17.9, 14.4, 11.4, | |
| 8.6, 6.2, 4.4, 3.0, 2.2, 2.4, 3.5, 1.7, -1.3, -4.2, -6.0, -5.4, -1.5, | |
| 6.0, 12.6, 13.9, 12.3 | |
| ]) | |
| Ln = phon | |
| Af = 4.47e-3 * (10 ** (.025 * Ln) - 1.15) \ | |
| + (.4 * 10 ** (((Tf + Lu) / 10.) - 9)) ** af | |
| Lp = ((10 / af) * np.log10(Af)) - Lu + 94 | |
| spl = Lp | |
| freq = f | |
| if return_freq is True: | |
| return np.array([freq, spl]) | |
| else: | |
| return spl | |
| def equal_loudness(phon, freqs, return_freq=False): | |
| """Returns equal-loudness levels for any frequencies between 20 Hz and | |
| 12.5 kHz according to the ISO/IEC 226 standard [1]_. | |
| Parameters | |
| ---------- | |
| phon : number | |
| The phon value represented by the equal-loudness contour, where a value | |
| of :math:`x` phon is the loudness of 1-KHz steady pure tone presented | |
| at :math:`x` dB SPL. Must be between 0 and 90. | |
| freqs : number or array_like | |
| Frequency or frequencies in Hz to be evaluated. Must be between 20 and | |
| 12500. | |
| return_freq : bool, optional | |
| If True, the function returns the frequency values as well as the SPL | |
| values of the contour. Default is False. | |
| Returns | |
| ------- | |
| array_like | |
| Either a 1-D or a 2-D numpy array, depending on `return_freq`. | |
| Reference | |
| --------- | |
| .. [1] ISO/IEC (2003). ISO/IEC 226:2003 Acoustics -- Normal equal-loudness- | |
| level contours. | |
| http://www.iso.org/iso/catalogue_detail.htm?csnumber=34222. | |
| Example | |
| ------- | |
| >>> el = equal_loudness(60, [500, 1000, 2000], return_freq=True) | |
| >>> print el | |
| [[ 500. 1000. 2000. ] | |
| [ 62.0511792 60.011588 59.96161453]] | |
| """ | |
| f = interp1d(*iso266(phon, True), kind='cubic') | |
| if return_freq is True: | |
| return np.array([freqs, f(freqs)]) | |
| else: | |
| return f(freqs) | |
| def get_loudness(spl, freq): | |
| """Returns the approximate loudness level in phons for a side-presented | |
| steady pure tone according to the ISO/IEC 226 standard [1]_. | |
| This function generates a range of equal-loudness contours and interpolates | |
| between them. Therefore it is more efficient to pass many level and | |
| frequency values to one function call than it is to make many function | |
| calls. | |
| Parameters | |
| ---------- | |
| spl : number or array_like | |
| Sound pressure level or levels in dB SPL. | |
| freq : number or array_like | |
| Frequency or frequencies in Hz. | |
| Returns | |
| ------- | |
| number or array_like | |
| Phon value(s). | |
| Reference | |
| --------- | |
| .. [1] ISO/IEC (2003). ISO/IEC 226:2003 Acoustics -- Normal equal-loudness- | |
| level contours. | |
| http://www.iso.org/iso/catalogue_detail.htm?csnumber=34222. | |
| Example | |
| ------- | |
| phons = get_loudness([50, 60, 70] [500, 500, 500]) | |
| print phons | |
| [ 47.3 57.8 68.4] | |
| """ | |
| phons = np.arange(0, 90.1, .1) | |
| freqs = np.arange(20, 12501) | |
| spls = np.empty((len(phons), len(freqs))) | |
| for i, phon in enumerate(phons): | |
| spls[i] = equal_loudness(phon, freqs) | |
| if not hasattr(spl, '__iter__'): | |
| spl = [spl] | |
| if not hasattr(freq, '__iter__'): | |
| freq = [freq] | |
| spls = spls.T | |
| results = [] | |
| for _spl, _freq in zip(spl, freq): | |
| ix = (np.abs(freqs - _freq)).argmin() | |
| iy = (np.abs(spls[ix] - _spl)).argmin() | |
| results.append(phons[iy]) | |
| if len(results) == 1: | |
| return results[0] | |
| else: | |
| return np.array(results) | |
| def pitch2freq(pitch): | |
| # from https://music.arts.uci.edu/dobrian/maxcookbook/pitch-and-loudness-formulae | |
| relative_pitch = pitch - 69 | |
| factor = 2 ** (relative_pitch / 12) | |
| freq = a_freq * factor | |
| return freq | |
| def velocity2amplitude(velocity): | |
| # from https://www.cs.cmu.edu/~rbd/papers/velocity-icmc2006.pdf | |
| r = 10 ** (upright_piano_dynamic_range / 20) | |
| b = 127 / (126 * np.sqrt(r)) - (1 / 126) | |
| m = (1 - b) / 127 | |
| a = (m * velocity + b) ** 2 | |
| a *= a_440_amplitude_at_max_velocity # scale amplitudes to get realistic perceived loudness | |
| return a | |
| def amplitude2db(amplitude): | |
| power_db = 20 * np.log10(amplitude) | |
| return power_db | |
| def get_db_of_equivalent_loudness_at_440hz(freqs, db): | |
| phons = get_loudness(db, freqs) | |
| equal_dbs = [] | |
| for p in phons: | |
| equal_dbs.append(equal_loudness(p, [440])[0]) | |
| return np.array(equal_dbs) | |
| def compute_total_loudness(eq_amplitudes_440hz, onsets, offsets): | |
| # Compute the instantaneous amplitude, turn it back to dbs, then to perceived loudness with unique freq 440 Hz | |
| # model amplitude as square function, loudness = peak amplitude from onset to offset, 0 afterwards. | |
| # an exponential model might be better | |
| assert all([len(values) == len(onsets) for values in [eq_amplitudes_440hz, offsets]]) | |
| timepoints = np.array(sorted(onsets + offsets)) | |
| amplitudes_per_time = np.zeros(len(timepoints)) | |
| # on each segment, compute the total amplitude | |
| # amplitudes are not just summed: p1+p2 = sqrt(p1**2 + p2**2) | |
| # ref: https://erlend-viggen.no/acoustic-quantities-1/ | |
| for i_n in range(len(onsets)): | |
| indexes = np.where(np.logical_and(timepoints >= onsets[i_n], timepoints < offsets[i_n])) | |
| amplitudes_per_time[indexes] += eq_amplitudes_440hz[i_n] ** 2 | |
| for i in range(len(amplitudes_per_time)): | |
| amplitudes_per_time[i] = np.sqrt(amplitudes_per_time[i]) | |
| # compute power | |
| power_per_time = amplitude2db(amplitudes_per_time) | |
| power_per_time[np.where(power_per_time == -np.inf)] = 0 | |
| # compute loudness | |
| loudness_per_time = get_loudness(power_per_time, [440] * len(power_per_time)) # amplitudes at 440hz, they were computed to make same loudness as original amplitudes at original F. | |
| # now integrate | |
| total_loudness = 0 | |
| for i_t in range(len(timepoints) - 1): | |
| total_loudness += loudness_per_time[i_t] * (timepoints[i_t + 1] - timepoints[i_t]) | |
| return total_loudness / (timepoints[-1] - timepoints[0]), np.std(loudness_per_time) | |
| if __name__ == '__main__': | |
| pass |