e@1: # -*- coding: utf-8 -*-
e@1: """
e@1: Created on Mon Mar  7 13:03:01 2016
e@1: 
e@1: @author: Emmanouil Theofanis Chourdakis
e@1: """
e@1: 
e@1: #from numpy import *
e@1: from scipy.optimize import *
e@1: from numpy import sqrt,array,polyval,exp,log,log10,mean,sum,cos,pi,arange,var,loadtxt
e@1: 
e@1: ga = sqrt(2)
e@1: g1_min = 0.001
e@1: g1_max = 0.999
e@1: 
e@1: d1_min = 0.010
e@1: d1_max = 0.900
e@1: 
e@1: da_min = 0.006
e@1: da_max = 0.012
e@1: 
e@1: gc_min = 0.001
e@1: gc_max = 0.999
e@1: 
e@1: G_min = 0.001
e@1: G_max = 0.999
e@1: 
e@1: # Limits are defined by the variables minima and maxima
e@1: 
e@1: T60_min = 0.2
e@1: T60_max = 4
e@1: #
e@1: D_min = 1000
e@1: D_max = 10000
e@1: #
e@1: C_min = -20
e@1: C_max = 10
e@1: #
e@1: Tc_min = 0.01 
e@1: Tc_max = 2.0
e@1: #
e@1: SC_min = 200
e@1: SC_max = 11025
e@1: def normalize(f,fmin,fmax):
e@1:     return (f - fmin)/(fmax - fmin) 
e@1: 
e@1: def normalize_param(p):
e@1:     return array([
e@1:             normalize(p[0],T60_min,T60_max),
e@1:             normalize(p[1],D_min,D_max),
e@1:             normalize(p[2],C_min,C_max),
e@1:             normalize(p[3],Tc_min,Tc_max)
e@1:             ])    
e@1: 
e@1: def get_high_param_norm(x):
e@1:     return array([
e@1:             normalize(T60(x),T60_min,T60_max),
e@1:             normalize(D(x),D_min,D_max),
e@1:             normalize(C(x),C_min,C_max),
e@1:             normalize(Tc(x),Tc_min,Tc_max)
e@1:             ])
e@1: def T601d(d1,g1,gc,G):
e@1:     g1 = abs(g1)
e@1:     return -d1/log(g1)*(log(G)+log(1-gc)+log(ga)+3*log(10))
e@1:     
e@1:     
e@1: 
e@1: def D1d(d1,da,t=0.1):
e@1:     kS = 20.78125
e@1:     
e@1:     return kS*t/d1/da
e@1:         
e@1:         
e@1: def C1d(g1,gc, G):
e@1:     s_ = 0
e@1:     for k in range(1, 7):
e@1:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@1:         
e@1:     return -10 * log10(G**2*(1-gc)/(1+gc)*s_)      
e@1:     
e@1: def Tc1d(d1,da,g1):
e@1:     gk = lambda k: g1**(1.5)**(1-k)
e@1:     dk = lambda k: d1*1.5**(1-k)
e@1: 
e@1:     sum1 = 0 
e@1:     sum2 = 0
e@1:     
e@1:     for k in range(1, 7):
e@1:         s1 = gk(k)**2/(1-gk(k)**2)
e@1:         sum2 += s1
e@1:         sum1 += dk(k)*s1/(1-gk(k)**2)
e@1:         
e@1:     return sum1/sum2 + da        
e@1:     
e@1: def Tc_g1_d1_D(d1,g1,D, t=0.1):
e@1:     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:     S2 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@1:     Tc = S1/S2 + 20.78125*t/(D*d1)
e@1:     return Tc
e@1: 
e@1: def T60_d1_g1_C(d1,g1,gc,C):
e@1: #    S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@1:     g1 = abs(g1)
e@1:     s_ = 0
e@1:     for k in range(1, 7):
e@1:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))    
e@1:       
e@1:     T60 = d1*log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))/log(g1)
e@1:     
e@1:     return T60 
e@1:     
e@1: def d1_g1_T60_C(g1,gc,T60,C):
e@1:     g1 = abs(g1)
e@1:     s_ = 0
e@1:     for k in range(1, 7):
e@1:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))      
e@1:     #S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@1:     
e@1:     
e@1:     d1 = T60*log(g1)/log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))
e@1:     return d1
e@1:     
e@1: def G_g1_C(g1,gc, C):
e@1:     
e@1:     s_ = 0
e@1:     for k in range(1, 7):
e@1:         s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
e@1:     
e@1:     G = sqrt(exp(-C/10*log(10))/(1-gc)*(1+gc)/s_)
e@1:     return G
e@1:     
e@1: def da_d1_D(d1, D,t=0.1):
e@1:     return 20.78125*t/d1/D
e@1:     
e@1:     
e@1: def Tc_g1_T60_D_C(g1,gc, T60, D, C,SR=44100.0,t=0.1):
e@1:     S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
e@1:     L1 = log(0.001*exp(C*log(10)/20)/(ga*sqrt((gc + 1)/((-gc + 1)*S1))*(-gc + 1)))
e@1:     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:     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:     
e@1:     return Tc
e@1: 
e@1: def sc_csv_to_c():
e@1:     T = loadtxt('scmapping.csv')
e@1:     SampleRates = T[:,0]
e@1:     P = T[:,1:]
e@1:     
e@1:     str_ = ""
e@1:     
e@1:     for SR in SampleRates:
e@1:         p = P[SampleRates == SR][0]
e@1:         str_ += "case %d:\n" % SR
e@1:         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:         str_ += "break;\n" 
e@1:     
e@1:     return str_
e@1: 
e@1: def f_gc(sc,SR=44100):
e@1:     T = loadtxt('scmapping.csv')
e@1:     SampleRates = T[:,0]
e@1:     P = T[:,1:]
e@1:     p = P[SampleRates == SR][0]
e@1:     return polyval(p, sc)
e@1:     
e@1:     
e@1: def errorfunc(g1, Target, gc, weights=(1.0, 1.0, 1.0, 1.0),SR=44100.0):
e@1:     
e@1:     (kT60, kD, kC, kTc) = Target
e@1:     G = G_g1_C(g1, gc, kC)
e@1: #     print "G1: ", g1
e@1: #     print "G: ", G
e@1:     
e@1:     
e@1:     # Check bounds
e@1: 
e@1:     if G>G_max:
e@1:         G = G_max
e@1:     if G<G_min:
e@1:         G = G_min
e@1:         
e@1: 
e@1:     d1 = d1_g1_T60_C(g1,gc, kT60, kC)
e@1:     
e@1: #     print "D1: ", d1
e@1:     
e@1:     if d1>d1_max:
e@1:         d1 = d1_max
e@1:     if d1<d1_min:
e@1:         d1 = d1_min
e@1:         
e@1:     da = da_d1_D(d1, kD)
e@1:     
e@1: #     print "DA: ", da
e@1:     
e@1:     if da>da_max:
e@1:         da = da_max
e@1:     if da<da_min:
e@1:        da = da_min
e@1:     
e@1:     T_60_est = T601d(d1,g1,gc,G)
e@1:     D_est = D1d(d1,da)
e@1:     C_est = C1d(g1,gc,G)
e@1:     Tc_est = Tc1d(d1,da,g1)
e@1:   #  print "TCEST:", Tc_est
e@1:   
e@1: #     print "T_60_est", T_60_est
e@1: #     print "D_est", D_est
e@1: #     print "C_est", C_est
e@1: #     print "Tc_est", Tc_est
e@1: #     
e@1: 
e@1:     est_vect = normalize_param([T_60_est, D_est, C_est, Tc_est])
e@1:     orig_vect = normalize_param([kT60, kD, kC, kTc])
e@1:     
e@1:     errvec = abs(est_vect-orig_vect)
e@1: #     print "ESTVEC: ", est_vect
e@1: #     print "ORIVEC: ", orig_vect
e@1: # 
e@1: #     print "ERRVEC: ", errvec
e@1:     w1 = weights[0]
e@1:     w2 = weights[1]
e@1:     w3 = weights[2]
e@1:     w4 = weights[3]
e@1: 
e@1:     errvec = array([w1*errvec[0],w2*errvec[1],w3*errvec[2],w4*errvec[3]])
e@1:     
e@1:     return mean(errvec**2) + var(errvec**2)**2 #+ skew(abs(est_vect-orig_vect)**2)
e@1:     
e@1: def SC1d(gc,SR=44100):
e@1:     
e@1:     lut = 1-2*gc*cos(2*pi*arange(0,int(SR/2)+1)/SR)+gc**2
e@1:     
e@1:     fgcn = lambda  n: 1/lut[n]
e@1: 
e@1:     sum1 = 0
e@1:     sum2 = 0
e@1:     
e@1:     for n in range(0, int(SR/2)):
e@1:         fgcn_ = fgcn(n)
e@1:         sum1 += n*fgcn_
e@1:         sum2 += fgcn_
e@1:         
e@1:     return sum1/sum2    
e@1: 
e@1: 
e@1: def forward_mapping(d1, da, g1, gc, G):
e@1:     SC = SC1d(gc)    
e@1:     T60 = T601d(d1,g1,gc,G)
e@1:     D = D1d(d1,da)
e@1:     C = C1d(g1,gc,G)
e@1:     Tc = Tc1d(d1,da,g1)
e@1:     
e@1:     return (T60, D, C, Tc, SC)
e@1:     
e@1: def inverse_mapping(T60, D, C, Tc, SC, SR=44100):
e@1:     """ returns a vector (d1, da, g1, gc, G) of filter coefficients """
e@1:         
e@1:     #1. Estimate gc from spectral centroid 
e@1:     gc = f_gc(SC, SR)
e@1: #     print "GC: ",gc
e@1:     
e@1:     #2. Estimate g1
e@1:     g1 = minimize_scalar(lambda x: errorfunc(x, (T60,D,C,Tc), gc), bounds=(g1_min,g1_max), method='bounded').x
e@1: #     print "G1: ", g1
e@1:     
e@1:     #4. a. Estimate G from G_g1_C
e@1:     G = G_g1_C(g1, gc, C)
e@1:     if G > G_max:
e@1:         G = G_max
e@1:     elif G < G_min:
e@1:         G = G_min
e@1: 
e@1: #     print "G: ", G
e@1: 
e@1:     #4. b. Estimate d1 from d1_g1_T60_C
e@1:     d1 = d1_g1_T60_C(g1,gc,T60,C)
e@1:     if d1 > d1_max:
e@1:         d1 = d1_max
e@1:     elif d1 < d1_min:
e@1:         d1 = d1_min
e@1:         
e@1: #     print "D1: ", d1
e@1:     #4. c. Estimate da from da_d1_D
e@1:     da = da_d1_D(d1,gc, D)
e@1:     if da > da_max:
e@1:         da = da_max
e@1:     elif da < da_min:
e@1:         da = da_min
e@1: 
e@1: #     print "DA: ", da
e@1:     return (d1, da, g1, gc, G)
e@1: 
e@1: #print forward_mapping(0.1, 0.006, 0.5, 0.5, 0.5)
e@1: 
e@1: #print forward_mapping(0.05, 0.006, 0.5, 0.5, 0.5)
e@1: 
e@1: #print Tc_g1_d1_D(0.1,0.5,1000, t=0.1)
e@1: #print T60_d1_g1_C(0.01, 0.9, 0.9, -5.0)
e@1: #print d1_g1_T60_C(0.9, 0.8, 4.2, -10.0)
e@1: #0.447907
e@1: #0.0788184
e@1: #0.945694
e@1: #print G_g1_C(0.2, 0.4, -2.0)
e@1: #print da_d1_D(0.05, 0.1, 100)
e@1: #print f_gc(11025, 44100)
e@1: #0.415625
e@1: #print sc_csv_to_c()
e@1: #-0.000462180326584
e@1: 
e@1: #0.980471
e@1: 
e@1: print "--"
e@1: #print
e@1: 
e@1: v= inverse_mapping(0.8465784284662089, 3463.5416666666665, -0.10779739192724623, 0.09213541733941949, 6419.2075162127485)
e@1: #v =  inverse_mapping(0.7, 3000, -5, 0.9, 5000)
e@1: print v
e@1: print forward_mapping(*v)
e@1: #print forward_mapping(0.15437302257477994, 0.012, 0.22280866837217597, 0.63966852775153527, 0.999)
e@1: # EST_VEC: 
e@1: # 0.170152218017423446 
e@1: # 0.27372685185185186 
e@1: # 0.663073420269091773 
e@1: # 0.0412740845527812192 
e@1: # 
e@1: # ORI_VEC: 
e@1: # 0.17015221801742339 
e@1: # 0.27372685185185186 
e@1: # 0.663073420269091773 
e@1: # 0.0412740790650349201 
e@1: # 
e@1: # ERR_VEC: 
e@1: # 5.5511151231257827e-17 
e@1: #           0 
e@1: #           0 
e@1: # 5.48774629904880129e-09 
e@1: 
e@1: # inverse mapping:
e@1: # 0.100257463495416962 
e@1: # 0.00299229593030461835 
e@1: # 0.498999999999999999 
e@1: # 0.499829818494971656 
e@1: # 0.500748281620112645 
e@1: 
e@1: #print "---"
e@1: #print errorfunc(0.4, (0.8465784284662089, 3463.5416666666665, -0.10779739192724623, 0.09213541733941949), 0.3)
e@1: #print T601d(0.30876416353,0.1,0.8,1.9511185262)
e@1: