-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathiso226.py
127 lines (110 loc) · 4.82 KB
/
iso226.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#!/usr/bin/env python3
"""
A python-ported version from Jeff Tackett Matlab script:
https://uk.mathworks.com/matlabcentral/fileexchange/
7028-iso-226-equal-loudness-level-contour-signal
"""
# function [spl, freq] = iso226(phon);
# %
# % Generates an Equal Loudness Contour as described in ISO 226
# %
# % Usage: [SPL FREQ] = ISO226(PHON);
# %
# % PHON is the phon value in dB SPL that you want the equal
# % loudness curve to represent. (1phon = 1dB @ 1kHz)
# % SPL is the Sound Pressure Level amplitude returned for
# % each of the 29 frequencies evaluated by ISO226.
# % FREQ is the returned vector of frequencies that ISO226
# % evaluates to generate the contour.
# %
# % Desc: This function will return the equal loudness contour for
# % your desired phon level. The frequencies evaulated in this
# % function only span from 20Hz - 12.5kHz, and only 29 selective
# % frequencies are covered. This is the limitation of the ISO
# % standard.
# %
# % In addition the valid phon range should be 0 - 90 dB SPL.
# % Values outside this range do not have experimental values
# % and their contours should be treated as inaccurate.
# %
# % If more samples are required you should be able to easily
# % interpolate these values using spline().
# %
# % Author: Jeff Tackett 03/01/05
#
#
#
# % /---------------------------------------\
# %%%%%%%%%%%%%%%%% TABLES FROM ISO226 %%%%%%%%%%%%%%%%%
# % \---------------------------------------/
# f = [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 = [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 = [-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 = [ 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];
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# %Error Trapping
# if((phon < 0) | (phon > 90))
# disp('Phon value out of bounds!')
# spl = 0;
# freq = 0;
# else
# %Setup user-defined values for equation
# Ln = phon;
#
# %Deriving sound pressure level from loudness level (iso226 sect 4.1)
# Af=4.47E-3 * (10.^(0.025*Ln) - 1.15) + (0.4*10.^(((Tf+Lu)/10)-9 )).^af;
# Lp=((10./af).*log10(Af)) - Lu + 94;
#
# %Return user data
# spl = Lp;
# freq = f;
# end
# ADDITIONAL INFO
# A more recent routine ISO 226:2003 has been published here:
# https://github.com/IoSR-Surrey/MatlabToolbox/blob/master/%2Biosr/%2Bauditory/iso226.m
# But the curves are pretty similar.
import numpy as np
FREQS = 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] )
def iso226(phon):
if phon < 0 or phon > 90:
raise ValueError('Phon value out of bounds!')
# Setup user-defined values for equation
Ln = phon;
# Deriving sound pressure level from loudness level (iso226 sect 4.1)
Af = 4.47E-3 * (10**(0.025 * Ln) - 1.15) + ( 0.4 * 10**(((Tf+Lu)/10)-9) ) ** af
Lp = ((10 / af) * np.log10(Af)) - Lu + 94
# Return user data
spl = Lp
return FREQS, spl
def make_eq_ld_curves():
""" equal loudness curves array from 0 phon to 90 phon
"""
eq_ld_curves = np.zeros( (91, len(FREQS)) )
for phon in np.arange(0, 91, 1):
_, mag = iso226(phon)
eq_ld_curves[phon] = mag
return eq_ld_curves
EQ_LD_CURVES = make_eq_ld_curves()