Mercurial > hg > chourdakisreiss2016

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

@author: Emmanouil Theofanis Chourdakis
"""

#from numpy import *
from scipy.optimize import *
from numpy import sqrt,array,polyval,exp,log,log10,mean,sum,cos,pi,arange,var,loadtxt

ga = sqrt(2)
g1_min = 0.001
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.2
T60_max = 4
#
D_min = 1000
D_max = 10000
#
C_min = -20
C_max = 10
#
Tc_min = 0.01
Tc_max = 2.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, D,t=0.1):
    return 20.78125*t/d1/D


def Tc_g1_T60_D_C(g1,gc, T60, D, C,SR=44100.0,t=0.1):
    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 sc_csv_to_c():
    T = loadtxt('scmapping.csv')
    SampleRates = T[:,0]
    P = T[:,1:]

    str_ = ""

    for SR in SampleRates:
        p = P[SampleRates == SR][0]
        str_ += "case %d:\n" % SR
        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])
        str_ += "break;\n"

    return str_

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),SR=44100.0):

    (kT60, kD, kC, kTc) = Target
    G = G_g1_C(g1, gc, kC)
#     print "G1: ", g1
#     print "G: ", G


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

#     print "D1: ", d1

    if d1>d1_max:
        d1 = d1_max
    if d1<d1_min:
        d1 = d1_min

    da = da_d1_D(d1, kD)

#     print "DA: ", da

    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)
  #  print "TCEST:", Tc_est

#     print "T_60_est", T_60_est
#     print "D_est", D_est
#     print "C_est", C_est
#     print "Tc_est", Tc_est
#

    est_vect = normalize_param([T_60_est, D_est, C_est, Tc_est])
    orig_vect = normalize_param([kT60, kD, kC, kTc])

    errvec = abs(est_vect-orig_vect)
#     print "ESTVEC: ", est_vect
#     print "ORIVEC: ", orig_vect
#
#     print "ERRVEC: ", errvec
    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)
#     print "GC: ",gc

    #2. Estimate g1
    g1 = minimize_scalar(lambda x: errorfunc(x, (T60,D,C,Tc), gc), bounds=(g1_min,g1_max), method='bounded').x
#     print "G1: ", g1

    #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

#     print "G: ", G

    #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

#     print "D1: ", d1
    #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

#     print "DA: ", da
    return (d1, da, g1, gc, G)

#print forward_mapping(0.1, 0.006, 0.5, 0.5, 0.5)

#print forward_mapping(0.05, 0.006, 0.5, 0.5, 0.5)

#print Tc_g1_d1_D(0.1,0.5,1000, t=0.1)
#print T60_d1_g1_C(0.01, 0.9, 0.9, -5.0)
#print d1_g1_T60_C(0.9, 0.8, 4.2, -10.0)
#0.447907
#0.0788184
#0.945694
#print G_g1_C(0.2, 0.4, -2.0)
#print da_d1_D(0.05, 0.1, 100)
#print f_gc(11025, 44100)
#0.415625
#print sc_csv_to_c()
#-0.000462180326584

#0.980471

print "--"
#print

v= inverse_mapping(0.8465784284662089, 3463.5416666666665, -0.10779739192724623, 0.09213541733941949, 6419.2075162127485)
#v =  inverse_mapping(0.7, 3000, -5, 0.9, 5000)
print v
print forward_mapping(*v)
#print forward_mapping(0.15437302257477994, 0.012, 0.22280866837217597, 0.63966852775153527, 0.999)
# EST_VEC:
# 0.170152218017423446
# 0.27372685185185186
# 0.663073420269091773
# 0.0412740845527812192
#
# ORI_VEC:
# 0.17015221801742339
# 0.27372685185185186
# 0.663073420269091773
# 0.0412740790650349201
#
# ERR_VEC:
# 5.5511151231257827e-17
#           0
#           0
# 5.48774629904880129e-09

# inverse mapping:
# 0.100257463495416962
# 0.00299229593030461835
# 0.498999999999999999
# 0.499829818494971656
# 0.500748281620112645

#print "---"
#print errorfunc(0.4, (0.8465784284662089, 3463.5416666666665, -0.10779739192724623, 0.09213541733941949), 0.3)
#print T601d(0.30876416353,0.1,0.8,1.9511185262)
author	Emmanouil Theofanis Chourdakis <e.t.chourdakis@qmul.ac.uk>
date	Sat, 30 Sep 2017 13:25:50 +0100
parents	144fbd1d29c3
children