view experiment-reverb/code/optim.py @ 2:c87a9505f294 tip

Added LICENSE for code, removed .wav files

# -*- 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()