e@0: # -*- coding: utf-8 -*-
e@0: """
e@0: Created on Mon Mar  7 13:03:01 2016
e@0: 
e@0: @author: Emmanouil Theofanis Chourdakis
e@0: """
e@0: 
e@0: from numpy import *
e@0: from scipy.optimize import *
e@0: 
e@0: 
e@0: ga = sqrt(2)
e@0: g1_min = 0.136
e@0: g1_max = 0.999
e@0: 
e@0: d1_min = 0.010
e@0: d1_max = 0.900
e@0: 
e@0: da_min = 0.006
e@0: da_max = 0.012
e@0: 
e@0: gc_min = 0.001
e@0: gc_max = 0.999
e@0: 
e@0: G_min = 0.001
e@0: G_max = 0.999
e@0: 
e@0: # Limits are defined by the variables minima and maxima
e@0: 
e@0: T60_min = 0.02
e@0: T60_max = 5
e@0: #
e@0: D_min = 200
e@0: D_max = 10000
e@0: #
e@0: C_min = -20
e@0: C_max = 10
e@0: #
e@0: Tc_min = 0.001 
e@0: Tc_max = 5.0
e@0: #
e@0: SC_min = 200
e@0: SC_max = 11025
e@0: def normalize(f,fmin,fmax):
e@0:     return (f - fmin)/(fmax - fmin) 
e@0: 
e@0: def normalize_param(p):
e@0:     return array([
e@0:             normalize(p[0],T60_min,T60_max),
e@0:             normalize(p[1],D_min,D_max),
e@0:             normalize(p[2],C_min,C_max),
e@0:             normalize(p[3],Tc_min,Tc_max)
e@0:             ])    
e@0: 
e@0: def get_high_param_norm(x):
e@0:     return array([
e@0:             normalize(T60(x),T60_min,T60_max),
e@0:             normalize(D(x),D_min,D_max),
e@0:             normalize(C(x),C_min,C_max),
e@0:             normalize(Tc(x),Tc_min,Tc_max)
e@0:             ])
e@0: def T601d(d1,g1,gc,G):
e@0:     g1 = abs(g1)
e@0:     return -d1/log(g1)*(log(G)+log(1-gc)+log(ga)+3*log(10))
e@0:     
e@0:     
e@0: 
e@0: def D1d(d1,da,t=0.1):
e@0:     kS = 20.78125
e@0:     
e@0:     return kS*t/d1/da
e@0:         
e@0:         
e@0: def C1d(g1,gc, G):
e@0:     s_ = 0
e@0:     for k in range(1, 7):
e@0:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@0:         
e@0:     return -10 * log10(G**2*(1-gc)/(1+gc)*s_)      
e@0:     
e@0: def Tc1d(d1,da,g1):
e@0:     gk = lambda k: g1**(1.5)**(1-k)
e@0:     dk = lambda k: d1*1.5**(1-k)
e@0: 
e@0:     sum1 = 0 
e@0:     sum2 = 0
e@0:     
e@0:     for k in range(1, 7):
e@0:         s1 = gk(k)**2/(1-gk(k)**2)
e@0:         sum2 += s1
e@0:         sum1 += dk(k)*s1/(1-gk(k)**2)
e@0:         
e@0:     return sum1/sum2 + da        
e@0:     
e@0: def Tc_g1_d1_D(d1,g1,D, t=0.1):
e@0:     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:     S2 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0:     Tc = S1/S2 + 20.78125*t/(D*d1)
e@0:     return Tc
e@0: 
e@0: def T60_d1_g1_C(d1,g1,gc,C):
e@0: #    S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0:     g1 = abs(g1)
e@0:     s_ = 0
e@0:     for k in range(1, 7):
e@0:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))    
e@0:       
e@0:     T60 = d1*log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))/log(g1)
e@0:     
e@0:     return T60 
e@0:     
e@0: def d1_g1_T60_C(g1,gc,T60,C):
e@0:     g1 = abs(g1)
e@0:     s_ = 0
e@0:     for k in range(1, 7):
e@0:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))      
e@0:     #S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0:     
e@0:     
e@0:     d1 = T60*log(g1)/log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))
e@0:     return d1
e@0:     
e@0: def G_g1_C(g1,gc, C):
e@0:     
e@0:     s_ = 0
e@0:     for k in range(1, 7):
e@0:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@0:     
e@0:     G = sqrt(exp(-C/10*log(10))/(1-gc)*(1+gc)/s_)
e@0:     return G
e@0:     
e@0: def da_d1_D(d1,gc, D,t=0.1):
e@0:     return 20.78125*t/d1/D
e@0:     
e@0:     
e@0: def Tc_g1_T60_D_C(g1,gc, T60, D, C):
e@0:     S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@0:     L1 = log(0.001*exp(C*log(10)/20)/(ga*sqrt((gc + 1)/((-gc + 1)*S1))*(-gc + 1)))
e@0:     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:     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:     
e@0:     return Tc
e@0: 
e@0: def f_gc(sc,SR=44100):
e@0:     T = loadtxt('scmapping.csv')
e@0:     SampleRates = T[:,0]
e@0:     P = T[:,1:]
e@0:     p = P[SampleRates == SR][0]
e@0:     return polyval(p, sc)
e@0: 
e@0: def errorfunc(g1, Target, gc, weights=(1.0, 1.0, 1.0, 1.0)):
e@0:     
e@0:     (kT60, kD, kC, kTc) = Target
e@0:     G = G_g1_C(g1, gc, kC)
e@0:     
e@0:     
e@0:     # Check bounds
e@0: 
e@0:     if G>G_max:
e@0:         G = G_max
e@0:     if G<G_min:
e@0:         G = G_min
e@0:         
e@0: 
e@0:     d1 = d1_g1_T60_C(g1,gc, kT60, kC)
e@0:     
e@0:     if d1>d1_max:
e@0:         d1 = d1_max
e@0:     if d1<d1_min:
e@0:         d1 = d1_min
e@0:         
e@0:     da = da_d1_D(d1, gc, kD)
e@0:     
e@0:     if da>da_max:
e@0:         da = da_max
e@0:     if da<da_min:
e@0:        da = da_min
e@0:     
e@0:     T_60_est = T601d(d1,g1,gc,G)
e@0:     D_est = D1d(d1,da)
e@0:     C_est = C1d(g1,gc,G)
e@0:     Tc_est = Tc1d(d1,da,g1)
e@0: 
e@0:     est_vect = normalize_param(array([T_60_est, D_est, C_est, Tc_est]))
e@0:     orig_vect = normalize_param([kT60, kD, kC, kTc])
e@0:     
e@0:     errvec = abs(est_vect-orig_vect)
e@0:     w1 = weights[0]
e@0:     w2 = weights[1]
e@0:     w3 = weights[2]
e@0:     w4 = weights[3]
e@0: 
e@0:     errvec = array([w1*errvec[0],w2*errvec[1],w3*errvec[2],w4*errvec[3]])
e@0:     
e@0:     return mean(errvec**2) + var(errvec**2)**2 #+ skew(abs(est_vect-orig_vect)**2)
e@0:     
e@0: def SC1d(gc,SR=44100):
e@0:     
e@0:     lut = 1-2*gc*cos(2*pi*arange(0,int(SR/2)+1)/SR)+gc**2
e@0:     
e@0:     fgcn = lambda  n: 1/lut[n]
e@0: 
e@0:     sum1 = 0
e@0:     sum2 = 0
e@0:     
e@0:     for n in range(0, int(SR/2)):
e@0:         fgcn_ = fgcn(n)
e@0:         sum1 += n*fgcn_
e@0:         sum2 += fgcn_
e@0:         
e@0:     return sum1/sum2    
e@0: 
e@0: 
e@0: def forward_mapping(d1, da, g1, gc, G):
e@0:     SC = SC1d(gc)    
e@0:     T60 = T601d(d1,g1,gc,G)
e@0:     D = D1d(d1,da)
e@0:     C = C1d(g1,gc,G)
e@0:     Tc = Tc1d(d1,da,g1)
e@0:     
e@0:     return (T60, D, C, Tc, SC)
e@0:     
e@0: def inverse_mapping(T60, D, C, Tc, SC, SR=44100):
e@0:     """ returns a vector (d1, da, g1, gc, G) of filter coefficients """
e@0:         
e@0:     #1. Estimate gc from spectral centroid 
e@0:     gc = f_gc(SC, SR)
e@0:     
e@0:     #2. Estimate g1
e@0:     g1 = minimize_scalar(lambda x: errorfunc(x, (T60,D,C,Tc), gc), bounds=(g1_min,g1_max), method='bounded').x
e@0:     
e@0:     #4. a. Estimate G from G_g1_C
e@0:     G = G_g1_C(g1, gc, C)
e@0:     if G > G_max:
e@0:         G = G_max
e@0:     elif G < G_min:
e@0:         G = G_min
e@0: 
e@0:     #4. b. Estimate d1 from d1_g1_T60_C
e@0:     d1 = d1_g1_T60_C(g1,gc,T60,C)
e@0:     if d1 > d1_max:
e@0:         d1 = d1_max
e@0:     elif d1 < d1_min:
e@0:         d1 = d1_min
e@0:         
e@0:     #4. c. Estimate da from da_d1_D
e@0:     da = da_d1_D(d1,gc, D)
e@0:     if da > da_max:
e@0:         da = da_max
e@0:     elif da < da_min:
e@0:         da = da_min
e@0: 
e@0:     return (d1, da, g1, gc, G)
e@0: