e@0: # -*- coding: utf-8 -*-
e@0: """
e@0: Created on Tue Jun 16 15:23:28 2015
e@0: 
e@0: @author: mmxgn
e@0: """
e@0: 
e@0: 
e@0: from matplotlib.pyplot import *
e@0: 
e@0: from numpy.core._internal import _gcd as gcd
e@0: 
e@0: from numpy import *
e@0: 
e@0: def calculate_parameters(d1,g1):
e@0:     
e@0:     d2 = int(round((1.5)**(-1)*d1))
e@0: 
e@0:     while gcd(d2,d1) != 1:
e@0:         d2 += 1
e@0:     
e@0:     d3 = int(round((1.5)**(-2)*d1))
e@0:     
e@0:     while gcd(d3, d2) != 1 or gcd(d3, d1) != 1:
e@0:         d3 += 1
e@0:         
e@0:     d4 = int(round((1.5)**(-3)*d1))
e@0:     
e@0:     while gcd(d4, d3) != 1 or gcd(d4, d2) != 1  or gcd(d4, d1) != 1:
e@0:         d4 += 1
e@0: 
e@0:         
e@0:     d5 = int(round((1.5)**(-4)*d1))
e@0: 
e@0:     while gcd(d5, d4) != 1 or gcd(d5, d3) != 1  or gcd(d5, d2) != 1 or gcd(d5, d1) != 1:
e@0:         d5 += 1        
e@0:     
e@0:     d6 = int(round((1.5)**(-5)*d1))
e@0:     while gcd(d6, d5) != 1 or gcd(d6, d4) != 1  or gcd(d6, d3) != 1 or gcd(d6, d2) != 1 or gcd(d6, d1) != 1:
e@0:         d6 += 1          
e@0:     g2 = g1**(1.5)**(-1)*g1
e@0:     g3 = g1**(1.5)**(-2)*g1
e@0:     g4 = g1**(1.5)**(-3)*g1     
e@0:     g5 = g1**(1.5)**(-4)*g1 
e@0:     g6 = g1**(1.5)**(-5)*g1        
e@0:     
e@0:     return (d1, d2, d3, d4, d5, d6, g1, g2, g3, g4, g5, g6)        
e@0: def estimate_SC(d1, g1, da, G, gc, SR):
e@0:     n = arange(0, SR/2+1)
e@0:     return sum(n/(1+gc**2-2*gc*cos(2*pi*n/SR)))/sum(1/(1+gc**2-2*gc*cos(2*pi*n/SR)))
e@0:     
e@0:     
e@0: def estimate_T60(d1, g1, da, G, gc, SR):
e@0:     ga = 1/sqrt(2)
e@0:     return d1/SR/log(g1)*log(10**-3/ga/(1-gc)/G)
e@0:         
e@0: def estimate_D(d1,g1, da, G, gc, SR):
e@0:     Dm = zeros((6,10))
e@0:     delays = zeros((6,))
e@0:     gains = zeros((6,1))        
e@0:     (delays[0],delays[1],delays[2],delays[3],delays[4],delays[5],gains[0],gains[1],gains[2],gains[3],gains[4],gains[5]) = calculate_parameters(d1,g1)    
e@0:     for k in range(1, 7):
e@0:         for m in range(1, 11):
e@0:             Dm[k-1,m-1] = max(0.1-m*delays[k-1]/SR,0)        
e@0:     
e@0:     return sum(Dm)
e@0: 
e@0:     
e@0: def estimate_C(d1, g1, da, G, gc, SR):
e@0:     g2 = g1**(1.5)**(-1)*g1
e@0:     g3 = g1**(1.5)**(-2)*g1
e@0:     g4 = g1**(1.5)**(-3)*g1     
e@0:     g5 = g1**(1.5)**(-4)*g1 
e@0:     g6 = g1**(1.5)**(-5)*g1         
e@0:     gains = zeros((6,1))
e@0:     gains[0] = g1
e@0:     gains[1] = g2
e@0:     gains[2] = g3
e@0:     gains[3] = g4
e@0:     gains[4] = g5
e@0:     gains[5] = g6
e@0:     
e@0:     return -10*log10(G**2*(1-gc)/(1+gc)*sum(1/(1-gains**2)))
e@0: 
e@0:     
e@0: def estimate_Tc(d1, g1, da, G, gc, SR):
e@0:     delays = zeros((6,))
e@0:     gains = zeros((6,1))        
e@0:     (delays[0],delays[1],delays[2],delays[3],delays[4],delays[5],gains[0],gains[1],gains[2],gains[3],gains[4],gains[5]) = calculate_parameters(d1,g1) 
e@0:     return sum(delays/SR*gains**2/(1-gains**2)**2)/sum(gains**2/(1-gains**2)) + da/SR        
e@0:         
e@0:     
e@0: GC = arange(0,1,0.01)
e@0: 
e@0: SR = 44100
e@0: SCy = array([estimate_SC(0,0,0,0,g,SR) for g in GC])
e@0: 
e@0: z = polyfit(GC,SCy,9)
e@0: 
e@0: 
e@0: #plot(GC,SCy)
e@0: estSC = poly1d(z)
e@0: #plot(GC, estSC(GC))
e@0: 
e@0: 
e@0: 
e@0: print mean((SCy-estSC(GC))**2)
e@0: #show()
e@0: from scipy.optimize import minimize_scalar
e@0: SCp = 8000
e@0: f = lambda x: abs(estSC(x) - SCp)
e@0: res = minimize_scalar(f, bounds=(0,1), method='bounded')
e@0: 
e@0: def utility_function(d1t,g1,dat,gc,G,SR,T60,T60min,T60max,D,Dmin,Dmax,C,Cmin,Cmax,Tc,Tcmin,Tcmax,SC,SCmin,SCmax):
e@0:     d1 = int(d1t*SR)
e@0:     da = int(dat*SR)
e@0:     T60p = estimate_T60(d1,g1,da,G,gc,SR)
e@0:     Dp = estimate_D(d1,g1,da,G,gc,SR)
e@0:     Cp = estimate_C(d1,g1,da,G,gc,SR)
e@0:     
e@0:     
e@0:     Tcp = estimate_Tc(d1,g1,da,G,gc,SR)
e@0:     SCp = estimate_SC(d1,g1,da,G,gc,SR)
e@0:     
e@0:     T60err = abs(T60p-T60)/(T60max-T60min)
e@0:     Derr = abs(Dp-D)/(Dmax-Dmin)
e@0:     Cerr = abs(Cp-C)/(Cmax-Cmin)
e@0:     Tcerr = abs(Tcp-Tc)/(Tcmax-Tcmin)
e@0:     Scerr = abs(SCp-SC)/(SCmax-SCmin)
e@0: 
e@0: 
e@0:     return T60err + Derr + Cerr + Tcerr + Scerr    
e@0:     
e@0: print utility_function(0.05, 0.5, 0.006, 0.5, 0.5, 44100, 0.373, 0.001, 120, 2841.75, 1, 8000, 2.181, -10, 20, 0.167, 0.02, 0.7, 3056.711, 10, SR/4)
e@0:     
e@0: 
e@0: g1x = arange(0.01,0.99,0.01)
e@0: #g1_utility = [utility_function(0.05, g1, 0.006, 0.5, 0.5, 44100, 0.373, 0.001, 120, 2841.75, 1, 8000, 2.181, -10, 20, 0.167, 0.02, 0.7, 3056.711, 10, SR/4) for g1 in arange(0.01,0.99,0.01)]
e@0: 
e@0: d1x = arange(0.01,0.1,0.001)*SR
e@0: #d1_utility = [utility_function(d1, 0.5, 0.006, 0.5, 0.5, 44100, 0.373, 0.001, 120, 2841.75, 1, 8000, 2.181, -10, 20, 0.167, 0.02, 0.7, 3056.711, 10, SR/4) for d1 in d1x]
e@0: 
e@0: dax = arange(0.006, 0.012, 0.001)*SR
e@0: #da_utility = [utility_function(0.05, 0.5, da, 0.5, 0.5, 44100, 0.373, 0.001, 120, 2841.75, 1, 8000, 2.181, -10, 20, 0.167, 0.02, 0.7, 3056.711, 10, SR/4) for da in dax]
e@0: 
e@0: gcx = arange(0.01, 0.99, 0.01)
e@0: #gc_utility = [utility_function(0.05, 0.5, 0.006, gc, 0.5, 44100, 0.373, 0.001, 120, 2841.75, 1, 8000, 2.181, -10, 20, 0.167, 0.02, 0.7, 3056.711, 10, SR/4) for gc in gcx]
e@0: 
e@0: 
e@0: flist = [estimate_T60, estimate_D, estimate_C, estimate_Tc, estimate_SC]
e@0: ix = 0
e@0: for func in flist:
e@0:     p1 = d1x
e@0:     p2 = g1x
e@0:     
e@0:     image = zeros((len(p1), len(p2)))
e@0:     for i in range(0, len(p1)):
e@0:         for j in range(0, len(p2)):
e@0:             image[i,j] = func(int(p1[i]), p2[j], int(0.012*SR), 0.5, 0.5, SR)
e@0:             figure()
e@0:             imshow(image)
e@0:             ylabel('g1')
e@0:             xlabel('d1')
e@0:             imsave(arr=image, fname='%d-g1-d1.pgf' % ix)
e@0:             
e@0: show()
e@0:            
e@0:         
e@0:             
e@0: 
e@0: 
e@0:     
e@0: 
e@0: