Mercurial > hg > chourdakisreiss2016

annotate experiment-reverb/code/mapping.py @ 2:c87a9505f294 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added LICENSE for code, removed .wav files
author Emmanouil Theofanis Chourdakis <e.t.chourdakis@qmul.ac.uk>
date Sat, 30 Sep 2017 13:25:50 +0100
parents 246d5546657c
children
rev   line source
e@0 1 # -*- coding: utf-8 -*-
e@0 2 """
e@0 3 Created on Mon Mar 7 13:03:01 2016
e@0 4
e@0 5 @author: Emmanouil Theofanis Chourdakis
e@0 6 """
e@0 7
e@0 8 from numpy import *
e@0 9 from scipy.optimize import *
e@0 10
e@0 11
e@0 12 ga = sqrt(2)
e@0 13 g1_min = 0.136
e@0 14 g1_max = 0.999
e@0 15
e@0 16 d1_min = 0.010
e@0 17 d1_max = 0.900
e@0 18
e@0 19 da_min = 0.006
e@0 20 da_max = 0.012
e@0 21
e@0 22 gc_min = 0.001
e@0 23 gc_max = 0.999
e@0 24
e@0 25 G_min = 0.001
e@0 26 G_max = 0.999
e@0 27
e@0 28 # Limits are defined by the variables minima and maxima
e@0 29
e@0 30 T60_min = 0.02
e@0 31 T60_max = 5
e@0 32 #
e@0 33 D_min = 200
e@0 34 D_max = 10000
e@0 35 #
e@0 36 C_min = -20
e@0 37 C_max = 10
e@0 38 #
e@0 39 Tc_min = 0.001
e@0 40 Tc_max = 5.0
e@0 41 #
e@0 42 SC_min = 200
e@0 43 SC_max = 11025
e@0 44 def normalize(f,fmin,fmax):
e@0 45 return (f - fmin)/(fmax - fmin)
e@0 46
e@0 47 def normalize_param(p):
e@0 48 return array([
e@0 49 normalize(p[0],T60_min,T60_max),
e@0 50 normalize(p[1],D_min,D_max),
e@0 51 normalize(p[2],C_min,C_max),
e@0 52 normalize(p[3],Tc_min,Tc_max)
e@0 53 ])
e@0 54
e@0 55 def get_high_param_norm(x):
e@0 56 return array([
e@0 57 normalize(T60(x),T60_min,T60_max),
e@0 58 normalize(D(x),D_min,D_max),
e@0 59 normalize(C(x),C_min,C_max),
e@0 60 normalize(Tc(x),Tc_min,Tc_max)
e@0 61 ])
e@0 62 def T601d(d1,g1,gc,G):
e@0 63 g1 = abs(g1)
e@0 64 return -d1/log(g1)*(log(G)+log(1-gc)+log(ga)+3*log(10))
e@0 65
e@0 66
e@0 67
e@0 68 def D1d(d1,da,t=0.1):
e@0 69 kS = 20.78125
e@0 70
e@0 71 return kS*t/d1/da
e@0 72
e@0 73
e@0 74 def C1d(g1,gc, G):
e@0 75 s_ = 0
e@0 76 for k in range(1, 7):
e@0 77 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@0 78
e@0 79 return -10 * log10(G**2*(1-gc)/(1+gc)*s_)
e@0 80
e@0 81 def Tc1d(d1,da,g1):
e@0 82 gk = lambda k: g1**(1.5)**(1-k)
e@0 83 dk = lambda k: d1*1.5**(1-k)
e@0 84
e@0 85 sum1 = 0
e@0 86 sum2 = 0
e@0 87
e@0 88 for k in range(1, 7):
e@0 89 s1 = gk(k)**2/(1-gk(k)**2)
e@0 90 sum2 += s1
e@0 91 sum1 += dk(k)*s1/(1-gk(k)**2)
e@0 92
e@0 93 return sum1/sum2 + da
e@0 94
e@0 95 def Tc_g1_d1_D(d1,g1,D, t=0.1):
e@0 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@0 97 S2 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0 98 Tc = S1/S2 + 20.78125*t/(D*d1)
e@0 99 return Tc
e@0 100
e@0 101 def T60_d1_g1_C(d1,g1,gc,C):
e@0 102 # S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0 103 g1 = abs(g1)
e@0 104 s_ = 0
e@0 105 for k in range(1, 7):
e@0 106 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@0 107
e@0 108 T60 = d1*log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))/log(g1)
e@0 109
e@0 110 return T60
e@0 111
e@0 112 def d1_g1_T60_C(g1,gc,T60,C):
e@0 113 g1 = abs(g1)
e@0 114 s_ = 0
e@0 115 for k in range(1, 7):
e@0 116 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@0 117 #S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0 118
e@0 119
e@0 120 d1 = T60*log(g1)/log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))
e@0 121 return d1
e@0 122
e@0 123 def G_g1_C(g1,gc, C):
e@0 124
e@0 125 s_ = 0
e@0 126 for k in range(1, 7):
e@0 127 s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@0 128
e@0 129 G = sqrt(exp(-C/10*log(10))/(1-gc)*(1+gc)/s_)
e@0 130 return G
e@0 131
e@0 132 def da_d1_D(d1,gc, D,t=0.1):
e@0 133 return 20.78125*t/d1/D
e@0 134
e@0 135
e@0 136 def Tc_g1_T60_D_C(g1,gc, T60, D, C):
e@0 137 S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0 138 L1 = log(0.001*exp(C*log(10)/20)/(ga*sqrt((gc + 1)/((-gc + 1)*S1))*(-gc + 1)))
e@0 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@0 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@0 141
e@0 142 return Tc
e@0 143
e@0 144 def f_gc(sc,SR=44100):
e@0 145 T = loadtxt('scmapping.csv')
e@0 146 SampleRates = T[:,0]
e@0 147 P = T[:,1:]
e@0 148 p = P[SampleRates == SR][0]
e@0 149 return polyval(p, sc)
e@0 150
e@0 151 def errorfunc(g1, Target, gc, weights=(1.0, 1.0, 1.0, 1.0)):
e@0 152
e@0 153 (kT60, kD, kC, kTc) = Target
e@0 154 G = G_g1_C(g1, gc, kC)
e@0 155
e@0 156
e@0 157 # Check bounds
e@0 158
e@0 159 if G>G_max:
e@0 160 G = G_max
e@0 161 if G<G_min:
e@0 162 G = G_min
e@0 163
e@0 164
e@0 165 d1 = d1_g1_T60_C(g1,gc, kT60, kC)
e@0 166
e@0 167 if d1>d1_max:
e@0 168 d1 = d1_max
e@0 169 if d1<d1_min:
e@0 170 d1 = d1_min
e@0 171
e@0 172 da = da_d1_D(d1, gc, kD)
e@0 173
e@0 174 if da>da_max:
e@0 175 da = da_max
e@0 176 if da<da_min:
e@0 177 da = da_min
e@0 178
e@0 179 T_60_est = T601d(d1,g1,gc,G)
e@0 180 D_est = D1d(d1,da)
e@0 181 C_est = C1d(g1,gc,G)
e@0 182 Tc_est = Tc1d(d1,da,g1)
e@0 183
e@0 184 est_vect = normalize_param(array([T_60_est, D_est, C_est, Tc_est]))
e@0 185 orig_vect = normalize_param([kT60, kD, kC, kTc])
e@0 186
e@0 187 errvec = abs(est_vect-orig_vect)
e@0 188 w1 = weights[0]
e@0 189 w2 = weights[1]
e@0 190 w3 = weights[2]
e@0 191 w4 = weights[3]
e@0 192
e@0 193 errvec = array([w1*errvec[0],w2*errvec[1],w3*errvec[2],w4*errvec[3]])
e@0 194
e@0 195 return mean(errvec**2) + var(errvec**2)**2 #+ skew(abs(est_vect-orig_vect)**2)
e@0 196
e@0 197 def SC1d(gc,SR=44100):
e@0 198
e@0 199 lut = 1-2*gc*cos(2*pi*arange(0,int(SR/2)+1)/SR)+gc**2
e@0 200
e@0 201 fgcn = lambda n: 1/lut[n]
e@0 202
e@0 203 sum1 = 0
e@0 204 sum2 = 0
e@0 205
e@0 206 for n in range(0, int(SR/2)):
e@0 207 fgcn_ = fgcn(n)
e@0 208 sum1 += n*fgcn_
e@0 209 sum2 += fgcn_
e@0 210
e@0 211 return sum1/sum2
e@0 212
e@0 213
e@0 214 def forward_mapping(d1, da, g1, gc, G):
e@0 215 SC = SC1d(gc)
e@0 216 T60 = T601d(d1,g1,gc,G)
e@0 217 D = D1d(d1,da)
e@0 218 C = C1d(g1,gc,G)
e@0 219 Tc = Tc1d(d1,da,g1)
e@0 220
e@0 221 return (T60, D, C, Tc, SC)
e@0 222
e@0 223 def inverse_mapping(T60, D, C, Tc, SC, SR=44100):
e@0 224 """ returns a vector (d1, da, g1, gc, G) of filter coefficients """
e@0 225
e@0 226 #1. Estimate gc from spectral centroid
e@0 227 gc = f_gc(SC, SR)
e@0 228
e@0 229 #2. Estimate g1
e@0 230 g1 = minimize_scalar(lambda x: errorfunc(x, (T60,D,C,Tc), gc), bounds=(g1_min,g1_max), method='bounded').x
e@0 231
e@0 232 #4. a. Estimate G from G_g1_C
e@0 233 G = G_g1_C(g1, gc, C)
e@0 234 if G > G_max:
e@0 235 G = G_max
e@0 236 elif G < G_min:
e@0 237 G = G_min
e@0 238
e@0 239 #4. b. Estimate d1 from d1_g1_T60_C
e@0 240 d1 = d1_g1_T60_C(g1,gc,T60,C)
e@0 241 if d1 > d1_max:
e@0 242 d1 = d1_max
e@0 243 elif d1 < d1_min:
e@0 244 d1 = d1_min
e@0 245
e@0 246 #4. c. Estimate da from da_d1_D
e@0 247 da = da_d1_D(d1,gc, D)
e@0 248 if da > da_max:
e@0 249 da = da_max
e@0 250 elif da < da_min:
e@0 251 da = da_min
e@0 252
e@0 253 return (d1, da, g1, gc, G)
e@0 254