|
e@1
|
1 # -*- coding: utf-8 -*-
|
|
e@1
|
2 """
|
|
e@1
|
3 Created on Mon Mar 7 13:03:01 2016
|
|
e@1
|
4
|
|
e@1
|
5 @author: Emmanouil Theofanis Chourdakis
|
|
e@1
|
6 """
|
|
e@1
|
7
|
|
e@1
|
8 #from numpy import *
|
|
e@1
|
9 from scipy.optimize import *
|
|
e@1
|
10 from numpy import sqrt,array,polyval,exp,log,log10,mean,sum,cos,pi,arange,var,loadtxt
|
|
e@1
|
11
|
|
e@1
|
12 ga = sqrt(2)
|
|
e@1
|
13 g1_min = 0.001
|
|
e@1
|
14 g1_max = 0.999
|
|
e@1
|
15
|
|
e@1
|
16 d1_min = 0.010
|
|
e@1
|
17 d1_max = 0.900
|
|
e@1
|
18
|
|
e@1
|
19 da_min = 0.006
|
|
e@1
|
20 da_max = 0.012
|
|
e@1
|
21
|
|
e@1
|
22 gc_min = 0.001
|
|
e@1
|
23 gc_max = 0.999
|
|
e@1
|
24
|
|
e@1
|
25 G_min = 0.001
|
|
e@1
|
26 G_max = 0.999
|
|
e@1
|
27
|
|
e@1
|
28 # Limits are defined by the variables minima and maxima
|
|
e@1
|
29
|
|
e@1
|
30 T60_min = 0.2
|
|
e@1
|
31 T60_max = 4
|
|
e@1
|
32 #
|
|
e@1
|
33 D_min = 1000
|
|
e@1
|
34 D_max = 10000
|
|
e@1
|
35 #
|
|
e@1
|
36 C_min = -20
|
|
e@1
|
37 C_max = 10
|
|
e@1
|
38 #
|
|
e@1
|
39 Tc_min = 0.01
|
|
e@1
|
40 Tc_max = 2.0
|
|
e@1
|
41 #
|
|
e@1
|
42 SC_min = 200
|
|
e@1
|
43 SC_max = 11025
|
|
e@1
|
44 def normalize(f,fmin,fmax):
|
|
e@1
|
45 return (f - fmin)/(fmax - fmin)
|
|
e@1
|
46
|
|
e@1
|
47 def normalize_param(p):
|
|
e@1
|
48 return array([
|
|
e@1
|
49 normalize(p[0],T60_min,T60_max),
|
|
e@1
|
50 normalize(p[1],D_min,D_max),
|
|
e@1
|
51 normalize(p[2],C_min,C_max),
|
|
e@1
|
52 normalize(p[3],Tc_min,Tc_max)
|
|
e@1
|
53 ])
|
|
e@1
|
54
|
|
e@1
|
55 def get_high_param_norm(x):
|
|
e@1
|
56 return array([
|
|
e@1
|
57 normalize(T60(x),T60_min,T60_max),
|
|
e@1
|
58 normalize(D(x),D_min,D_max),
|
|
e@1
|
59 normalize(C(x),C_min,C_max),
|
|
e@1
|
60 normalize(Tc(x),Tc_min,Tc_max)
|
|
e@1
|
61 ])
|
|
e@1
|
62 def T601d(d1,g1,gc,G):
|
|
e@1
|
63 g1 = abs(g1)
|
|
e@1
|
64 return -d1/log(g1)*(log(G)+log(1-gc)+log(ga)+3*log(10))
|
|
e@1
|
65
|
|
e@1
|
66
|
|
e@1
|
67
|
|
e@1
|
68 def D1d(d1,da,t=0.1):
|
|
e@1
|
69 kS = 20.78125
|
|
e@1
|
70
|
|
e@1
|
71 return kS*t/d1/da
|
|
e@1
|
72
|
|
e@1
|
73
|
|
e@1
|
74 def C1d(g1,gc, G):
|
|
e@1
|
75 s_ = 0
|
|
e@1
|
76 for k in range(1, 7):
|
|
e@1
|
77 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
|
|
e@1
|
78
|
|
e@1
|
79 return -10 * log10(G**2*(1-gc)/(1+gc)*s_)
|
|
e@1
|
80
|
|
e@1
|
81 def Tc1d(d1,da,g1):
|
|
e@1
|
82 gk = lambda k: g1**(1.5)**(1-k)
|
|
e@1
|
83 dk = lambda k: d1*1.5**(1-k)
|
|
e@1
|
84
|
|
e@1
|
85 sum1 = 0
|
|
e@1
|
86 sum2 = 0
|
|
e@1
|
87
|
|
e@1
|
88 for k in range(1, 7):
|
|
e@1
|
89 s1 = gk(k)**2/(1-gk(k)**2)
|
|
e@1
|
90 sum2 += s1
|
|
e@1
|
91 sum1 += dk(k)*s1/(1-gk(k)**2)
|
|
e@1
|
92
|
|
e@1
|
93 return sum1/sum2 + da
|
|
e@1
|
94
|
|
e@1
|
95 def Tc_g1_d1_D(d1,g1,D, t=0.1):
|
|
e@1
|
96 S1 = sum([1.5**(-k + 1)*d1*g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1)**2 for k in range(1, 7)])
|
|
e@1
|
97 S2 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
|
|
e@1
|
98 Tc = S1/S2 + 20.78125*t/(D*d1)
|
|
e@1
|
99 return Tc
|
|
e@1
|
100
|
|
e@1
|
101 def T60_d1_g1_C(d1,g1,gc,C):
|
|
e@1
|
102 # S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
|
|
e@1
|
103 g1 = abs(g1)
|
|
e@1
|
104 s_ = 0
|
|
e@1
|
105 for k in range(1, 7):
|
|
e@1
|
106 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
|
|
e@1
|
107
|
|
e@1
|
108 T60 = d1*log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))/log(g1)
|
|
e@1
|
109
|
|
e@1
|
110 return T60
|
|
e@1
|
111
|
|
e@1
|
112 def d1_g1_T60_C(g1,gc,T60,C):
|
|
e@1
|
113 g1 = abs(g1)
|
|
e@1
|
114 s_ = 0
|
|
e@1
|
115 for k in range(1, 7):
|
|
e@1
|
116 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
|
|
e@1
|
117 #S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
|
|
e@1
|
118
|
|
e@1
|
119
|
|
e@1
|
120 d1 = T60*log(g1)/log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))
|
|
e@1
|
121 return d1
|
|
e@1
|
122
|
|
e@1
|
123 def G_g1_C(g1,gc, C):
|
|
e@1
|
124
|
|
e@1
|
125 s_ = 0
|
|
e@1
|
126 for k in range(1, 7):
|
|
e@1
|
127 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
|
|
e@1
|
128
|
|
e@1
|
129 G = sqrt(exp(-C/10*log(10))/(1-gc)*(1+gc)/s_)
|
|
e@1
|
130 return G
|
|
e@1
|
131
|
|
e@1
|
132 def da_d1_D(d1, D,t=0.1):
|
|
e@1
|
133 return 20.78125*t/d1/D
|
|
e@1
|
134
|
|
e@1
|
135
|
|
e@1
|
136 def Tc_g1_T60_D_C(g1,gc, T60, D, C,SR=44100.0,t=0.1):
|
|
e@1
|
137 S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
|
|
e@1
|
138 L1 = log(0.001*exp(C*log(10)/20)/(ga*sqrt((gc + 1)/((-gc + 1)*S1))*(-gc + 1)))
|
|
e@1
|
139 S2 = sum([1.5**(-k + 1)*T60*g1**(2*1.5**(-k + 1))*log(g1)/((-g1**(2*1.5**(-k + 1)) + 1)**2*L1) for k in range(1, 7)])
|
|
e@1
|
140 Tc = S2/S1 + 20.78125*t*log(0.001*exp(C*log(10)/20)/(ga*sqrt((gc + 1)/((-gc + 1)*S1))*(-gc + 1)))/(D*T60*log(g1))
|
|
e@1
|
141
|
|
e@1
|
142 return Tc
|
|
e@1
|
143
|
|
e@1
|
144 def sc_csv_to_c():
|
|
e@1
|
145 T = loadtxt('scmapping.csv')
|
|
e@1
|
146 SampleRates = T[:,0]
|
|
e@1
|
147 P = T[:,1:]
|
|
e@1
|
148
|
|
e@1
|
149 str_ = ""
|
|
e@1
|
150
|
|
e@1
|
151 for SR in SampleRates:
|
|
e@1
|
152 p = P[SampleRates == SR][0]
|
|
e@1
|
153 str_ += "case %d:\n" % SR
|
|
e@1
|
154 str_ += "\tp0 = %1.18e;\n\tp1 = %1.18e;\n\tp2 = %1.18e;\n\tp3 = %1.18e;\n\tp4 = %1.18e;\n" % (p[0], p[1], p[2], p[3], p[4])
|
|
e@1
|
155 str_ += "break;\n"
|
|
e@1
|
156
|
|
e@1
|
157 return str_
|
|
e@1
|
158
|
|
e@1
|
159 def f_gc(sc,SR=44100):
|
|
e@1
|
160 T = loadtxt('scmapping.csv')
|
|
e@1
|
161 SampleRates = T[:,0]
|
|
e@1
|
162 P = T[:,1:]
|
|
e@1
|
163 p = P[SampleRates == SR][0]
|
|
e@1
|
164 return polyval(p, sc)
|
|
e@1
|
165
|
|
e@1
|
166
|
|
e@1
|
167 def errorfunc(g1, Target, gc, weights=(1.0, 1.0, 1.0, 1.0),SR=44100.0):
|
|
e@1
|
168
|
|
e@1
|
169 (kT60, kD, kC, kTc) = Target
|
|
e@1
|
170 G = G_g1_C(g1, gc, kC)
|
|
e@1
|
171 # print "G1: ", g1
|
|
e@1
|
172 # print "G: ", G
|
|
e@1
|
173
|
|
e@1
|
174
|
|
e@1
|
175 # Check bounds
|
|
e@1
|
176
|
|
e@1
|
177 if G>G_max:
|
|
e@1
|
178 G = G_max
|
|
e@1
|
179 if G<G_min:
|
|
e@1
|
180 G = G_min
|
|
e@1
|
181
|
|
e@1
|
182
|
|
e@1
|
183 d1 = d1_g1_T60_C(g1,gc, kT60, kC)
|
|
e@1
|
184
|
|
e@1
|
185 # print "D1: ", d1
|
|
e@1
|
186
|
|
e@1
|
187 if d1>d1_max:
|
|
e@1
|
188 d1 = d1_max
|
|
e@1
|
189 if d1<d1_min:
|
|
e@1
|
190 d1 = d1_min
|
|
e@1
|
191
|
|
e@1
|
192 da = da_d1_D(d1, kD)
|
|
e@1
|
193
|
|
e@1
|
194 # print "DA: ", da
|
|
e@1
|
195
|
|
e@1
|
196 if da>da_max:
|
|
e@1
|
197 da = da_max
|
|
e@1
|
198 if da<da_min:
|
|
e@1
|
199 da = da_min
|
|
e@1
|
200
|
|
e@1
|
201 T_60_est = T601d(d1,g1,gc,G)
|
|
e@1
|
202 D_est = D1d(d1,da)
|
|
e@1
|
203 C_est = C1d(g1,gc,G)
|
|
e@1
|
204 Tc_est = Tc1d(d1,da,g1)
|
|
e@1
|
205 # print "TCEST:", Tc_est
|
|
e@1
|
206
|
|
e@1
|
207 # print "T_60_est", T_60_est
|
|
e@1
|
208 # print "D_est", D_est
|
|
e@1
|
209 # print "C_est", C_est
|
|
e@1
|
210 # print "Tc_est", Tc_est
|
|
e@1
|
211 #
|
|
e@1
|
212
|
|
e@1
|
213 est_vect = normalize_param([T_60_est, D_est, C_est, Tc_est])
|
|
e@1
|
214 orig_vect = normalize_param([kT60, kD, kC, kTc])
|
|
e@1
|
215
|
|
e@1
|
216 errvec = abs(est_vect-orig_vect)
|
|
e@1
|
217 # print "ESTVEC: ", est_vect
|
|
e@1
|
218 # print "ORIVEC: ", orig_vect
|
|
e@1
|
219 #
|
|
e@1
|
220 # print "ERRVEC: ", errvec
|
|
e@1
|
221 w1 = weights[0]
|
|
e@1
|
222 w2 = weights[1]
|
|
e@1
|
223 w3 = weights[2]
|
|
e@1
|
224 w4 = weights[3]
|
|
e@1
|
225
|
|
e@1
|
226 errvec = array([w1*errvec[0],w2*errvec[1],w3*errvec[2],w4*errvec[3]])
|
|
e@1
|
227
|
|
e@1
|
228 return mean(errvec**2) + var(errvec**2)**2 #+ skew(abs(est_vect-orig_vect)**2)
|
|
e@1
|
229
|
|
e@1
|
230 def SC1d(gc,SR=44100):
|
|
e@1
|
231
|
|
e@1
|
232 lut = 1-2*gc*cos(2*pi*arange(0,int(SR/2)+1)/SR)+gc**2
|
|
e@1
|
233
|
|
e@1
|
234 fgcn = lambda n: 1/lut[n]
|
|
e@1
|
235
|
|
e@1
|
236 sum1 = 0
|
|
e@1
|
237 sum2 = 0
|
|
e@1
|
238
|
|
e@1
|
239 for n in range(0, int(SR/2)):
|
|
e@1
|
240 fgcn_ = fgcn(n)
|
|
e@1
|
241 sum1 += n*fgcn_
|
|
e@1
|
242 sum2 += fgcn_
|
|
e@1
|
243
|
|
e@1
|
244 return sum1/sum2
|
|
e@1
|
245
|
|
e@1
|
246
|
|
e@1
|
247 def forward_mapping(d1, da, g1, gc, G):
|
|
e@1
|
248 SC = SC1d(gc)
|
|
e@1
|
249 T60 = T601d(d1,g1,gc,G)
|
|
e@1
|
250 D = D1d(d1,da)
|
|
e@1
|
251 C = C1d(g1,gc,G)
|
|
e@1
|
252 Tc = Tc1d(d1,da,g1)
|
|
e@1
|
253
|
|
e@1
|
254 return (T60, D, C, Tc, SC)
|
|
e@1
|
255
|
|
e@1
|
256 def inverse_mapping(T60, D, C, Tc, SC, SR=44100):
|
|
e@1
|
257 """ returns a vector (d1, da, g1, gc, G) of filter coefficients """
|
|
e@1
|
258
|
|
e@1
|
259 #1. Estimate gc from spectral centroid
|
|
e@1
|
260 gc = f_gc(SC, SR)
|
|
e@1
|
261 # print "GC: ",gc
|
|
e@1
|
262
|
|
e@1
|
263 #2. Estimate g1
|
|
e@1
|
264 g1 = minimize_scalar(lambda x: errorfunc(x, (T60,D,C,Tc), gc), bounds=(g1_min,g1_max), method='bounded').x
|
|
e@1
|
265 # print "G1: ", g1
|
|
e@1
|
266
|
|
e@1
|
267 #4. a. Estimate G from G_g1_C
|
|
e@1
|
268 G = G_g1_C(g1, gc, C)
|
|
e@1
|
269 if G > G_max:
|
|
e@1
|
270 G = G_max
|
|
e@1
|
271 elif G < G_min:
|
|
e@1
|
272 G = G_min
|
|
e@1
|
273
|
|
e@1
|
274 # print "G: ", G
|
|
e@1
|
275
|
|
e@1
|
276 #4. b. Estimate d1 from d1_g1_T60_C
|
|
e@1
|
277 d1 = d1_g1_T60_C(g1,gc,T60,C)
|
|
e@1
|
278 if d1 > d1_max:
|
|
e@1
|
279 d1 = d1_max
|
|
e@1
|
280 elif d1 < d1_min:
|
|
e@1
|
281 d1 = d1_min
|
|
e@1
|
282
|
|
e@1
|
283 # print "D1: ", d1
|
|
e@1
|
284 #4. c. Estimate da from da_d1_D
|
|
e@1
|
285 da = da_d1_D(d1,gc, D)
|
|
e@1
|
286 if da > da_max:
|
|
e@1
|
287 da = da_max
|
|
e@1
|
288 elif da < da_min:
|
|
e@1
|
289 da = da_min
|
|
e@1
|
290
|
|
e@1
|
291 # print "DA: ", da
|
|
e@1
|
292 return (d1, da, g1, gc, G)
|
|
e@1
|
293
|
|
e@1
|
294 #print forward_mapping(0.1, 0.006, 0.5, 0.5, 0.5)
|
|
e@1
|
295
|
|
e@1
|
296 #print forward_mapping(0.05, 0.006, 0.5, 0.5, 0.5)
|
|
e@1
|
297
|
|
e@1
|
298 #print Tc_g1_d1_D(0.1,0.5,1000, t=0.1)
|
|
e@1
|
299 #print T60_d1_g1_C(0.01, 0.9, 0.9, -5.0)
|
|
e@1
|
300 #print d1_g1_T60_C(0.9, 0.8, 4.2, -10.0)
|
|
e@1
|
301 #0.447907
|
|
e@1
|
302 #0.0788184
|
|
e@1
|
303 #0.945694
|
|
e@1
|
304 #print G_g1_C(0.2, 0.4, -2.0)
|
|
e@1
|
305 #print da_d1_D(0.05, 0.1, 100)
|
|
e@1
|
306 #print f_gc(11025, 44100)
|
|
e@1
|
307 #0.415625
|
|
e@1
|
308 #print sc_csv_to_c()
|
|
e@1
|
309 #-0.000462180326584
|
|
e@1
|
310
|
|
e@1
|
311 #0.980471
|
|
e@1
|
312
|
|
e@1
|
313 print "--"
|
|
e@1
|
314 #print
|
|
e@1
|
315
|
|
e@1
|
316 v= inverse_mapping(0.8465784284662089, 3463.5416666666665, -0.10779739192724623, 0.09213541733941949, 6419.2075162127485)
|
|
e@1
|
317 #v = inverse_mapping(0.7, 3000, -5, 0.9, 5000)
|
|
e@1
|
318 print v
|
|
e@1
|
319 print forward_mapping(*v)
|
|
e@1
|
320 #print forward_mapping(0.15437302257477994, 0.012, 0.22280866837217597, 0.63966852775153527, 0.999)
|
|
e@1
|
321 # EST_VEC:
|
|
e@1
|
322 # 0.170152218017423446
|
|
e@1
|
323 # 0.27372685185185186
|
|
e@1
|
324 # 0.663073420269091773
|
|
e@1
|
325 # 0.0412740845527812192
|
|
e@1
|
326 #
|
|
e@1
|
327 # ORI_VEC:
|
|
e@1
|
328 # 0.17015221801742339
|
|
e@1
|
329 # 0.27372685185185186
|
|
e@1
|
330 # 0.663073420269091773
|
|
e@1
|
331 # 0.0412740790650349201
|
|
e@1
|
332 #
|
|
e@1
|
333 # ERR_VEC:
|
|
e@1
|
334 # 5.5511151231257827e-17
|
|
e@1
|
335 # 0
|
|
e@1
|
336 # 0
|
|
e@1
|
337 # 5.48774629904880129e-09
|
|
e@1
|
338
|
|
e@1
|
339 # inverse mapping:
|
|
e@1
|
340 # 0.100257463495416962
|
|
e@1
|
341 # 0.00299229593030461835
|
|
e@1
|
342 # 0.498999999999999999
|
|
e@1
|
343 # 0.499829818494971656
|
|
e@1
|
344 # 0.500748281620112645
|
|
e@1
|
345
|
|
e@1
|
346 #print "---"
|
|
e@1
|
347 #print errorfunc(0.4, (0.8465784284662089, 3463.5416666666665, -0.10779739192724623, 0.09213541733941949), 0.3)
|
|
e@1
|
348 #print T601d(0.30876416353,0.1,0.8,1.9511185262)
|
|
e@1
|
349
|
