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

Added LICENSE for code, removed .wav files

# -*- coding: utf-8 -*-
"""
Created on Mon Mar  7 13:03:01 2016

@author: Emmanouil Theofanis Chourdakis
"""

from numpy import *
from scipy.optimize import *


ga = sqrt(2)
g1_min = 0.136
g1_max = 0.999

d1_min = 0.010
d1_max = 0.900

da_min = 0.006
da_max = 0.012

gc_min = 0.001
gc_max = 0.999

G_min = 0.001
G_max = 0.999

# Limits are defined by the variables minima and maxima

T60_min = 0.02
T60_max = 5
#
D_min = 200
D_max = 10000
#
C_min = -20
C_max = 10
#
Tc_min = 0.001 
Tc_max = 5.0
#
SC_min = 200
SC_max = 11025
def normalize(f,fmin,fmax):
    return (f - fmin)/(fmax - fmin) 

def normalize_param(p):
    return array([
            normalize(p[0],T60_min,T60_max),
            normalize(p[1],D_min,D_max),
            normalize(p[2],C_min,C_max),
            normalize(p[3],Tc_min,Tc_max)
            ])    

def get_high_param_norm(x):
    return array([
            normalize(T60(x),T60_min,T60_max),
            normalize(D(x),D_min,D_max),
            normalize(C(x),C_min,C_max),
            normalize(Tc(x),Tc_min,Tc_max)
            ])
def T601d(d1,g1,gc,G):
    g1 = abs(g1)
    return -d1/log(g1)*(log(G)+log(1-gc)+log(ga)+3*log(10))
    
    

def D1d(d1,da,t=0.1):
    kS = 20.78125
    
    return kS*t/d1/da
        
        
def C1d(g1,gc, G):
    s_ = 0
    for k in range(1, 7):
        s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
        
    return -10 * log10(G**2*(1-gc)/(1+gc)*s_)      
    
def Tc1d(d1,da,g1):
    gk = lambda k: g1**(1.5)**(1-k)
    dk = lambda k: d1*1.5**(1-k)

    sum1 = 0 
    sum2 = 0
    
    for k in range(1, 7):
        s1 = gk(k)**2/(1-gk(k)**2)
        sum2 += s1
        sum1 += dk(k)*s1/(1-gk(k)**2)
        
    return sum1/sum2 + da        
    
def Tc_g1_d1_D(d1,g1,D, t=0.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)])
    S2 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
    Tc = S1/S2 + 20.78125*t/(D*d1)
    return Tc

def T60_d1_g1_C(d1,g1,gc,C):
#    S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
    g1 = abs(g1)
    s_ = 0
    for k in range(1, 7):
        s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))    
      
    T60 = d1*log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))/log(g1)
    
    return T60 
    
def d1_g1_T60_C(g1,gc,T60,C):
    g1 = abs(g1)
    s_ = 0
    for k in range(1, 7):
        s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))      
    #S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
    
    
    d1 = T60*log(g1)/log(0.001/ga*exp(C*log(10)/20)/(sqrt((gc + 1)/((-gc + 1)*s_))*(-gc + 1)))
    return d1
    
def G_g1_C(g1,gc, C):
    
    s_ = 0
    for k in range(1, 7):
        s_ += g1**(2*1.5**(1-k))/(1-g1**(2*1.5**(1-k)))
    
    G = sqrt(exp(-C/10*log(10))/(1-gc)*(1+gc)/s_)
    return G
    
def da_d1_D(d1,gc, D,t=0.1):
    return 20.78125*t/d1/D
    
    
def Tc_g1_T60_D_C(g1,gc, T60, D, C):
    S1 = sum([g1**(2*1.5**(-k + 1))/(-g1**(2*1.5**(-k + 1)) + 1) for k in range(1, 7)])
    L1 = log(0.001*exp(C*log(10)/20)/(ga*sqrt((gc + 1)/((-gc + 1)*S1))*(-gc + 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)])
    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))
    
    return Tc

def f_gc(sc,SR=44100):
    T = loadtxt('scmapping.csv')
    SampleRates = T[:,0]
    P = T[:,1:]
    p = P[SampleRates == SR][0]
    return polyval(p, sc)

def errorfunc(g1, Target, gc, weights=(1.0, 1.0, 1.0, 1.0)):
    
    (kT60, kD, kC, kTc) = Target
    G = G_g1_C(g1, gc, kC)
    
    
    # Check bounds

    if G>G_max:
        G = G_max
    if G<G_min:
        G = G_min
        

    d1 = d1_g1_T60_C(g1,gc, kT60, kC)
    
    if d1>d1_max:
        d1 = d1_max
    if d1<d1_min:
        d1 = d1_min
        
    da = da_d1_D(d1, gc, kD)
    
    if da>da_max:
        da = da_max
    if da<da_min:
       da = da_min
    
    T_60_est = T601d(d1,g1,gc,G)
    D_est = D1d(d1,da)
    C_est = C1d(g1,gc,G)
    Tc_est = Tc1d(d1,da,g1)

    est_vect = normalize_param(array([T_60_est, D_est, C_est, Tc_est]))
    orig_vect = normalize_param([kT60, kD, kC, kTc])
    
    errvec = abs(est_vect-orig_vect)
    w1 = weights[0]
    w2 = weights[1]
    w3 = weights[2]
    w4 = weights[3]

    errvec = array([w1*errvec[0],w2*errvec[1],w3*errvec[2],w4*errvec[3]])
    
    return mean(errvec**2) + var(errvec**2)**2 #+ skew(abs(est_vect-orig_vect)**2)
    
def SC1d(gc,SR=44100):
    
    lut = 1-2*gc*cos(2*pi*arange(0,int(SR/2)+1)/SR)+gc**2
    
    fgcn = lambda  n: 1/lut[n]

    sum1 = 0
    sum2 = 0
    
    for n in range(0, int(SR/2)):
        fgcn_ = fgcn(n)
        sum1 += n*fgcn_
        sum2 += fgcn_
        
    return sum1/sum2    


def forward_mapping(d1, da, g1, gc, G):
    SC = SC1d(gc)    
    T60 = T601d(d1,g1,gc,G)
    D = D1d(d1,da)
    C = C1d(g1,gc,G)
    Tc = Tc1d(d1,da,g1)
    
    return (T60, D, C, Tc, SC)
    
def inverse_mapping(T60, D, C, Tc, SC, SR=44100):
    """ returns a vector (d1, da, g1, gc, G) of filter coefficients """
        
    #1. Estimate gc from spectral centroid 
    gc = f_gc(SC, SR)
    
    #2. Estimate g1
    g1 = minimize_scalar(lambda x: errorfunc(x, (T60,D,C,Tc), gc), bounds=(g1_min,g1_max), method='bounded').x
    
    #4. a. Estimate G from G_g1_C
    G = G_g1_C(g1, gc, C)
    if G > G_max:
        G = G_max
    elif G < G_min:
        G = G_min

    #4. b. Estimate d1 from d1_g1_T60_C
    d1 = d1_g1_T60_C(g1,gc,T60,C)
    if d1 > d1_max:
        d1 = d1_max
    elif d1 < d1_min:
        d1 = d1_min
        
    #4. c. Estimate da from da_d1_D
    da = da_d1_D(d1,gc, D)
    if da > da_max:
        da = da_max
    elif da < da_min:
        da = da_min

    return (d1, da, g1, gc, G)