diff experiment-reverb/code/optim.py @ 0:246d5546657c

initial commit, needs cleanup

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