Mercurial > hg > chourdakisreiss2016

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added LICENSE for code, removed .wav files
author Emmanouil Theofanis Chourdakis <e.t.chourdakis@qmul.ac.uk>
date Sat, 30 Sep 2017 13:25:50 +0100
parents 246d5546657c
children
rev   line source
e@0 1 # -*- coding: utf-8 -*-
e@0 2 """
e@0 3 Created on Tue Jun 16 15:23:28 2015
e@0 4
e@0 5 @author: mmxgn
e@0 6 """
e@0 7
e@0 8
e@0 9 from matplotlib.pyplot import *
e@0 10
e@0 11 from numpy.core._internal import _gcd as gcd
e@0 12
e@0 13 from numpy import *
e@0 14
e@0 15 def calculate_parameters(d1,g1):
e@0 16
e@0 17 d2 = int(round((1.5)**(-1)*d1))
e@0 18
e@0 19 while gcd(d2,d1) != 1:
e@0 20 d2 += 1
e@0 21
e@0 22 d3 = int(round((1.5)**(-2)*d1))
e@0 23
e@0 24 while gcd(d3, d2) != 1 or gcd(d3, d1) != 1:
e@0 25 d3 += 1
e@0 26
e@0 27 d4 = int(round((1.5)**(-3)*d1))
e@0 28
e@0 29 while gcd(d4, d3) != 1 or gcd(d4, d2) != 1 or gcd(d4, d1) != 1:
e@0 30 d4 += 1
e@0 31
e@0 32
e@0 33 d5 = int(round((1.5)**(-4)*d1))
e@0 34
e@0 35 while gcd(d5, d4) != 1 or gcd(d5, d3) != 1 or gcd(d5, d2) != 1 or gcd(d5, d1) != 1:
e@0 36 d5 += 1
e@0 37
e@0 38 d6 = int(round((1.5)**(-5)*d1))
e@0 39 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 40 d6 += 1
e@0 41 g2 = g1**(1.5)**(-1)*g1
e@0 42 g3 = g1**(1.5)**(-2)*g1
e@0 43 g4 = g1**(1.5)**(-3)*g1
e@0 44 g5 = g1**(1.5)**(-4)*g1
e@0 45 g6 = g1**(1.5)**(-5)*g1
e@0 46
e@0 47 return (d1, d2, d3, d4, d5, d6, g1, g2, g3, g4, g5, g6)
e@0 48 def estimate_SC(d1, g1, da, G, gc, SR):
e@0 49 n = arange(0, SR/2+1)
e@0 50 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 51
e@0 52
e@0 53 def estimate_T60(d1, g1, da, G, gc, SR):
e@0 54 ga = 1/sqrt(2)
e@0 55 return d1/SR/log(g1)*log(10**-3/ga/(1-gc)/G)
e@0 56
e@0 57 def estimate_D(d1,g1, da, G, gc, SR):
e@0 58 Dm = zeros((6,10))
e@0 59 delays = zeros((6,))
e@0 60 gains = zeros((6,1))
e@0 61 (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 62 for k in range(1, 7):
e@0 63 for m in range(1, 11):
e@0 64 Dm[k-1,m-1] = max(0.1-m*delays[k-1]/SR,0)
e@0 65
e@0 66 return sum(Dm)
e@0 67
e@0 68
e@0 69 def estimate_C(d1, g1, da, G, gc, SR):
e@0 70 g2 = g1**(1.5)**(-1)*g1
e@0 71 g3 = g1**(1.5)**(-2)*g1
e@0 72 g4 = g1**(1.5)**(-3)*g1
e@0 73 g5 = g1**(1.5)**(-4)*g1
e@0 74 g6 = g1**(1.5)**(-5)*g1
e@0 75 gains = zeros((6,1))
e@0 76 gains[0] = g1
e@0 77 gains[1] = g2
e@0 78 gains[2] = g3
e@0 79 gains[3] = g4
e@0 80 gains[4] = g5
e@0 81 gains[5] = g6
e@0 82
e@0 83 return -10*log10(G**2*(1-gc)/(1+gc)*sum(1/(1-gains**2)))
e@0 84
e@0 85
e@0 86 def estimate_Tc(d1, g1, da, G, gc, SR):
e@0 87 delays = zeros((6,))
e@0 88 gains = zeros((6,1))
e@0 89 (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 90 return sum(delays/SR*gains**2/(1-gains**2)**2)/sum(gains**2/(1-gains**2)) + da/SR
e@0 91
e@0 92
e@0 93 GC = arange(0,1,0.01)
e@0 94
e@0 95 SR = 44100
e@0 96 SCy = array([estimate_SC(0,0,0,0,g,SR) for g in GC])
e@0 97
e@0 98 z = polyfit(GC,SCy,9)
e@0 99
e@0 100
e@0 101 #plot(GC,SCy)
e@0 102 estSC = poly1d(z)
e@0 103 #plot(GC, estSC(GC))
e@0 104
e@0 105
e@0 106
e@0 107 print mean((SCy-estSC(GC))**2)
e@0 108 #show()
e@0 109 from scipy.optimize import minimize_scalar
e@0 110 SCp = 8000
e@0 111 f = lambda x: abs(estSC(x) - SCp)
e@0 112 res = minimize_scalar(f, bounds=(0,1), method='bounded')
e@0 113
e@0 114 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 115 d1 = int(d1t*SR)
e@0 116 da = int(dat*SR)
e@0 117 T60p = estimate_T60(d1,g1,da,G,gc,SR)
e@0 118 Dp = estimate_D(d1,g1,da,G,gc,SR)
e@0 119 Cp = estimate_C(d1,g1,da,G,gc,SR)
e@0 120
e@0 121
e@0 122 Tcp = estimate_Tc(d1,g1,da,G,gc,SR)
e@0 123 SCp = estimate_SC(d1,g1,da,G,gc,SR)
e@0 124
e@0 125 T60err = abs(T60p-T60)/(T60max-T60min)
e@0 126 Derr = abs(Dp-D)/(Dmax-Dmin)
e@0 127 Cerr = abs(Cp-C)/(Cmax-Cmin)
e@0 128 Tcerr = abs(Tcp-Tc)/(Tcmax-Tcmin)
e@0 129 Scerr = abs(SCp-SC)/(SCmax-SCmin)
e@0 130
e@0 131
e@0 132 return T60err + Derr + Cerr + Tcerr + Scerr
e@0 133
e@0 134 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 135
e@0 136
e@0 137 g1x = arange(0.01,0.99,0.01)
e@0 138 #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 139
e@0 140 d1x = arange(0.01,0.1,0.001)*SR
e@0 141 #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 142
e@0 143 dax = arange(0.006, 0.012, 0.001)*SR
e@0 144 #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 145
e@0 146 gcx = arange(0.01, 0.99, 0.01)
e@0 147 #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 148
e@0 149
e@0 150 flist = [estimate_T60, estimate_D, estimate_C, estimate_Tc, estimate_SC]
e@0 151 ix = 0
e@0 152 for func in flist:
e@0 153 p1 = d1x
e@0 154 p2 = g1x
e@0 155
e@0 156 image = zeros((len(p1), len(p2)))
e@0 157 for i in range(0, len(p1)):
e@0 158 for j in range(0, len(p2)):
e@0 159 image[i,j] = func(int(p1[i]), p2[j], int(0.012*SR), 0.5, 0.5, SR)
e@0 160 figure()
e@0 161 imshow(image)
e@0 162 ylabel('g1')
e@0 163 xlabel('d1')
e@0 164 imsave(arr=image, fname='%d-g1-d1.pgf' % ix)
e@0 165
e@0 166 show()
e@0 167
e@0 168
e@0 169
e@0 170
e@0 171
e@0 172
e@0 173
e@0 174