Mercurial > hg > chourdakisreiss2016

comparison cpp/chourdakisreiss2016cpp/Parameter.cpp @ 1:144fbd1d29c3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added c++ class Parameter
author Emmanouil Theofanis Chourdakis <e.t.chourdakis@qmul.ac.uk>
date Wed, 28 Dec 2016 17:22:49 +0000
parents
children
comparison
equal deleted inserted replaced
0:246d5546657c 1:144fbd1d29c3
1 /*
2 * Parameter.cpp
3 *
4 * Created on: Dec 14, 2016
5 * Author: mmxgn
6 */
7
8 #include <cmath>
9 #include "Parameter.h"
10 //#include "csv.h"
11
12
13 /*
14 * Given a k, return g_k and d_k based on
15 * g_1 and d_1 respectively.
16 *
17 */
18
19 #define GK(k) powf(G1, powf(1.5, 1-k))
20 #define _GK(k) powf(_G1, powf(1.5, 1-k))
21
22 #define DK(k) (D1*powf(1.5, 1-k))
23 #define _DK(k) (_D1*powf(1.5, 1-k))
24
25 /*
26 * The parameter class
27 */
28
29 Parameter::Parameter(double _SR,
30 double _D1,
31 double _DA,
32 double _G1,
33 double _GC,
34 double _G,
35 double _T60,
36 double _ED,
37 double _C,
38 double _TC,
39 double _SC)
40 :
41 D1(_D1),
42 DA(_DA),
43 G1(_G1),
44 GC(_GC),
45 G(_G),
46 T60(_T60),
47 ED(_ED),
48 C(_C),
49 TC(_TC),
50 SC(_SC),
51 SC_MAX(_SR/4.0),
52 SampleRate(_SR)
53 {
54 // pass
55 }
56
57 double Parameter::getC() const {
58 return C;
59 }
60
61 void Parameter::setC(double c) {
62 C = c;
63 }
64
65 double Parameter::getD1() const {
66 return D1;
67 }
68
69 void Parameter::setD1(double d1) {
70 /* Dependency for:
71 * T60, ED, TC
72 */
73
74 T60 = compute_T60(d1, G1, GC, G);
75 ED = compute_ED(d1, DA);
76 TC = compute_TC(d1, G1, DA);
77
78 }
79
80 double Parameter::getDa() const {
81 return DA;
82
83 }
84
85 void Parameter::setDa(double da) {
86 /* Dependency for:
87 * ED, TC
88 */
89
90 ED = compute_ED(D1, da);
91 TC = compute_TC(D1, G1, da);
92
93 }
94
95 double Parameter::getEd() const {
96 return ED;
97 }
98
99 void Parameter::setEd(double ed) {
100 matrix <double,5,1> est_param = inverse_mapping(T60, ed, C, TC, SC, SampleRate,G1);
101 D1 = est_param(0);
102 DA = est_param(1);
103 G1 = est_param(2);
104 GC = est_param(3);
105 G = est_param(4);
106 T60 = compute_T60(D1, G1, GC, G);
107 ED = compute_ED(D1, DA);
108 C = compute_C(G1, GC, G);
109 TC = compute_TC(D1, G1, DA);
110 SC = compute_SC(GC, SampleRate);
111 }
112
113 double Parameter::getG() const {
114 return G;
115 }
116
117 void Parameter::setG(double g) {
118 /*
119 * Dependency for:
120 * T60, C
121 */
122
123 T60 = compute_T60(D1, G1, GC, g);
124 C = compute_C(G1, GC, g);
125
126 }
127
128 double Parameter::getG1() const {
129 return G1;
130 }
131
132 void Parameter::setG1(double g1) {
133 /*
134 * Dependency for:
135 * T60, C, TC
136 */
137
138 T60 = compute_T60(D1, g1, GC, G);
139 C = compute_C(g1, GC, G);
140 TC = compute_TC(D1, g1, DA);
141
142 }
143
144 double Parameter::getGc() const {
145 return GC;
146 }
147
148 void Parameter::setGc(double gc) {
149 /*
150 * Dependency for:
151 * SC
152 */
153
154 SC = compute_SC(gc, SampleRate);
155 }
156
157 double Parameter::getSampleRate() const {
158 return SampleRate;
159 }
160
161 void Parameter::setSampleRate(double sampleRate) {
162 SampleRate = sampleRate;
163 set_low_and_compute(D1, DA, G1, GC, G);
164 matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, TC, SC, SampleRate,G1);
165 D1 = est_param(0);
166 DA = est_param(1);
167 G1 = est_param(2);
168 GC = est_param(3);
169 G = est_param(4);
170 T60 = compute_T60(D1, G1, GC, G);
171 ED = compute_ED(D1, DA);
172 C = compute_C(G1, GC, G);
173 TC = compute_TC(D1, G1, DA);
174 SC = compute_SC(GC, SampleRate);
175 }
176
177 double Parameter::getSc() const {
178 return SC;
179 }
180
181 void Parameter::setSc(double sc) {
182 matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, TC, sc, SampleRate,G1);
183 D1 = est_param(0);
184 DA = est_param(1);
185 G1 = est_param(2);
186 GC = est_param(3);
187 G = est_param(4);
188 T60 = compute_T60(D1, G1, GC, G);
189 ED = compute_ED(D1, DA);
190 C = compute_C(G1, GC, G);
191 TC = compute_TC(D1, G1, DA);
192 SC = compute_SC(GC, SampleRate);
193 }
194
195 double Parameter::getT60() const {
196 return T60;
197 }
198
199 void Parameter::setT60(double t60) {
200 matrix <double,5,1> est_param = inverse_mapping(t60, ED, C, TC, SC, SampleRate,G1);
201 D1 = est_param(0);
202 DA = est_param(1);
203 G1 = est_param(2);
204 GC = est_param(3);
205 G = est_param(4);
206 T60 = compute_T60(D1, G1, GC, G);
207 ED = compute_ED(D1, DA);
208 C = compute_C(G1, GC, G);
209 TC = compute_TC(D1, G1, DA);
210 SC = compute_SC(GC, SampleRate);
211 }
212
213 double Parameter::getTc() const {
214 return TC;
215 }
216
217 void Parameter::setTc(double tc) {
218 matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, tc, SC, SampleRate,G1);
219 D1 = est_param(0);
220 DA = est_param(1);
221 G1 = est_param(2);
222 GC = est_param(3);
223 G = est_param(4);
224 T60 = compute_T60(D1, G1, GC, G);
225 ED = compute_ED(D1, DA);
226 C = compute_C(G1, GC, G);
227 TC = compute_TC(D1, G1, DA);
228 SC = compute_SC(GC, SampleRate);
229 }
230
231 Parameter::~Parameter() {
232 // TODO Auto-generated destructor stub
233 }
234
235 std::ostream& operator<<(std::ostream &s, const Parameter &p) {
236 return s << "Parameter({High:<" <<
237 p.T60 << ", " <<
238 p.ED << ", " <<
239 p.C << ", " <<
240 p.TC << ", " <<
241 p.SC << ">}, " <<
242 "{Low: <" <<
243 p.D1 << ", " <<
244 p.DA << ", " <<
245 p.G1 << ", " <<
246 p.GC << ", " <<
247 p.G << ">}).";
248
249 }
250
251 int Parameter::set_low_and_compute(double D1, double DA, double G1, double GC,
252 double G) {
253
254 /* Temporary values */
255
256 double _T60, _ED, _C, _TC, _SC;
257
258 _T60 = compute_T60(D1, G1, GC, G);
259 if (_T60 < T60_MIN || _T60 > T60_MAX)
260 return -1;
261
262 _ED = compute_ED(D1, DA);
263 if (_ED < ED_MIN || _ED > ED_MAX)
264 return -1;
265
266 _C = compute_C(G1, GC, G);
267
268 if (_C < C_MIN || _C > C_MAX)
269 return -1;
270
271 _TC = compute_TC(D1, G1, DA);
272
273 if (_TC < TC_MIN || _TC > TC_MAX)
274 return -1;
275
276 _SC = compute_SC(GC, SampleRate);
277
278 if (_SC < SC_MIN || _SC > SC_MAX)
279 return -1;
280
281 T60 = _T60;
282 ED = _ED;
283 C = _C;
284 TC = _TC;
285 SC = _SC;
286
287 return 0;
288
289 }
290
291 double Parameter::compute_T60(double _D1, double _G1, double _GC, double _G) {
292 return -_D1/log(_G1)*(log(_G)+log(1-_GC)+log(GA)+3*log(10));
293 }
294
295 double Parameter::compute_ED(double _D1, double _DA) {
296 return 2.078125/_D1/_DA;
297 }
298
299 double Parameter::compute_C(double _G1, double _GC, double _G) {
300 double sum1 = 0;
301 for (int k=1; k<7; k++){
302 sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1, 2*powf(1.5, 1-k)));
303 }
304 return -10*log10(_G*_G*(1-_GC)/(1+_GC)*sum1);
305 }
306
307 double Parameter::compute_TC(double _D1, double _G1, double _DA) {
308 double sum1 = 0, sum2 = 0, s1 = 0;
309 for (int k=1; k<7; k++){
310 s1 = powf(_GK(k), 2.0)/(1-powf(_GK(k), 2.0));
311 sum2 += s1;
312 sum1 += _DK(k)*s1/(1-powf(_GK(k),2.0));
313 }
314
315 return sum1/sum2 + _DA;
316 }
317
318 double Parameter::compute_SC(double _GC, double _SR) {
319 double sum1 = 0, sum2 = 0;
320
321 for (int k=0; k < int(_SR/2.0); k++){
322 double lut = 1.0-2.0*_GC*cos(2.0*M_PI*k/_SR)+_GC*_GC;
323 sum1 += k/lut;
324 sum2 += 1/lut;
325 }
326
327 return sum1/sum2;
328 }
329
330 int Parameter::get_low_and_computer() {
331 return 0;
332 }
333
334 int Parameter::set_high_and_compute(double T60, double ED, double C, double TC,
335 double SC) {
336 return 0;
337 }
338
339 double Parameter::TC_G1_D1_ED(double _D1, double _G1, double _ED, double t) {
340 double sum1 = 0, sum2 = 0;
341
342 for (int k=1; k<7; k++){
343 sum1 += powf(1.5, -k+1)*_D1*
344 powf(_G1, 2*powf(1.5, -k+1))/
345 powf((-powf(_G1, 2*powf(1.5, -k+1))+1),2);
346 sum2 += powf(_G1,2*powf(1.5,-k+1))/
347 (-powf(_G1, 2*powf(1.5,-k+1))+1);
348 }
349
350 return sum1/sum2 + 20.78125*t/(_ED*_D1);
351 }
352
353 double Parameter::T60_D1_G1_C(double _D1, double _G1, double _GC, double _C) {
354 double sum1 = 0;
355
356 for (int k=1; k<7; k++){
357 sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1,2*powf(1.5,1-k)));
358 }
359
360 return _D1*log(0.001/GA*exp(_C*log(10)/20)/
361 (sqrt((1+_GC)/((1-_GC)*sum1))*(1-_GC)))/log(_G1);
362 }
363
364 double Parameter::D1_G1_T60_C(double _G1, double _GC, double _T60, double _C) {
365 double sum1 = 0;
366
367 for (int k=1; k<7; k++){
368 sum1 += powf(_G1,(2*powf(1.5,1-k)))/(1-powf(_G1,2*powf(1.5,1-k)));
369 }
370
371 return _T60*log(_G1)/log(0.001/GA*exp(_C*log(10)/20)/(sqrt((1+_GC)/((1-_GC)*sum1))*(1-_GC)));
372 }
373
374 double Parameter::G_G1_C(double _G1, double _GC, double _C) {
375 double sum1 = 0;
376 for (int k=1; k<7; k++){
377 sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1, 2*powf(1.5, 1-k)));
378 }
379
380 return sqrt(exp(-_C/10*log(10))/(1-_GC)*(1+_GC)/sum1);
381 }
382
383 double Parameter::DA_D1_ED(double _D1, double _ED, double t) {
384 return 20.78125*t/_D1/_ED;
385 }
386
387
388 double Parameter::GC_SC(double _SC, double SR) {
389 double p0 = 1.0,
390 p1 = 1.0,
391 p2 = 1.0,
392 p3 = 1.0,
393 p4 = 1.0;
394
395 switch (int(SR)){
396 case 8000:
397 p0 = -7.545401867799958007e-15;
398 p1 = 6.605438530163032344e-11;
399 p2 = -2.339787304666883204e-07;
400 p3 = -2.380996920823082359e-04;
401 p4 = 1.003686549974644349e+00;
402 break;
403 case 11025:
404 p0 = -2.090495384250846006e-15;
405 p1 = 2.523110858472787412e-11;
406 p2 = -1.232004017234833354e-07;
407 p3 = -1.727859517219043628e-04;
408 p4 = 1.003773048218970398e+00;
409 break;
410 case 16000:
411 p0 = -4.707928872668820147e-16;
412 p1 = 8.250344830381294683e-12;
413 p2 = -5.847797160444276075e-08;
414 p3 = -1.190647306137786752e-04;
415 p4 = 1.003845141060948354e+00;
416 break;
417 case 21000:
418 p0 = -1.585759862171924958e-16;
419 p1 = 3.648228536317020713e-12;
420 p2 = -3.394367375513113085e-08;
421 p3 = -9.071942141117647239e-05;
422 p4 = 1.003883355464454619e+00;
423 break;
424 case 22050:
425 p0 = -1.304517046624577766e-16;
426 p1 = 3.151372288676199373e-12;
427 p2 = -3.078752354037690531e-08;
428 p3 = -8.639997174340922025e-05;
429 p4 = 1.003889194220798364e+00;
430 break;
431 case 32000:
432 p0 = -2.939615050414309854e-17;
433 p1 = 1.030821612883252373e-12;
434 p2 = -1.461695000370224873e-08;
435 p3 = -5.953732000720860076e-05;
436 p4 = 1.003925597252757207e+00;
437 break;
438 case 37800:
439 p0 = -1.509578131570550335e-17;
440 p1 = 6.253538622497773237e-13;
441 p2 = -1.047516110428303656e-08;
442 p3 = -5.040267062802067523e-05;
443 p4 = 1.003938012276286873e+00;
444 break;
445 case 44056:
446 p0 = -8.179918755257568798e-18;
447 p1 = 3.949666985075376660e-13;
448 p2 = -7.711241596699441441e-09;
449 p3 = -4.324592160336278773e-05;
450 p4 = 1.003947751979878067e+00;
451 break;
452 case 44100:
453 p0 = -8.147316012596786372e-18;
454 p1 = 3.937855238616550567e-13;
455 p2 = -7.695860699409773702e-09;
456 p3 = -4.320277669373632828e-05;
457 p4 = 1.003947810730627443e+00;
458 break;
459 case 47250:
460 p0 = -6.182156102335834534e-18;
461 p1 = 3.201545709180980074e-13;
462 p2 = -6.703887839395126650e-09;
463 p3 = -4.032277597168298445e-05;
464 p4 = 1.003951733380953559e+00;
465 break;
466 case 48000:
467 p0 = -5.804669671441499519e-18;
468 p1 = 3.053789928757690291e-13;
469 p2 = -6.496014866101249794e-09;
470 p3 = -3.969277247952757479e-05;
471 p4 = 1.003952591709977415e+00;
472 break;
473 case 50000:
474 p0 = -4.930039877651698485e-18;
475 p1 = 2.701762164752109105e-13;
476 p2 = -5.986696651033372113e-09;
477 p3 = -3.810515833478767255e-05;
478 p4 = 1.003954755096216900e+00;
479 break;
480 case 50400:
481 p0 = -4.775358854140600092e-18;
482 p1 = 2.637937027466497354e-13;
483 p2 = -5.892041055138477995e-09;
484 p3 = -3.780275477003837485e-05;
485 p4 = 1.003955167234265478e+00;
486 break;
487 case 88200:
488 p0 = -5.090139768480031043e-19;
489 p1 = 4.921424693356908116e-14;
490 p2 = -1.923829523675108346e-09;
491 p3 = -2.160214687014651151e-05;
492 p4 = 1.003977276279834019e+00;
493 break;
494 case 96000:
495 p0 = -3.626648836744412987e-19;
496 p1 = 3.816597595693573777e-14;
497 p2 = -1.623898000952206408e-09;
498 p3 = -1.984703070895103917e-05;
499 p4 = 1.003979674995812754e+00;
500 break;
501 case 176400:
502 p0 = -3.180720104052761169e-20;
503 p1 = 6.151207927460715513e-15;
504 p2 = -4.809399117661735116e-10;
505 p3 = -1.080127073758081756e-05;
506 p4 = 1.003992048707243789e+00;
507 break;
508 case 192000:
509 p0 = -2.266250752087159745e-20;
510 p1 = 4.770337959813220317e-15;
511 p2 = -4.059609201723455372e-10;
512 p3 = -9.923682425709201677e-06;
513 p4 = 1.003993250132436454e+00;
514 break;
515 case 352800:
516 p0 = -1.987755090799793513e-21;
517 p1 = 7.688647436953224411e-16;
518 p2 = -1.202327621003883834e-10;
519 p3 = -5.400685657755166908e-06;
520 p4 = 1.003999444875642189e+00;
521 break;
522 case 2822400:
523 p0 = -4.852498403591686415e-25;
524 p1 = 1.501626536939688520e-18;
525 p2 = -1.878606337281021417e-12;
526 p3 = -6.750912735944337505e-07;
527 p4 = 1.004005921979573834e+00;
528 break;
529 case 5644800:
530 p0 = -3.032792721276874507e-26;
531 p1 = 1.877027577372803603e-19;
532 p2 = -4.696510359527020341e-13;
533 p3 = -3.375458367765195566e-07;
534 p4 = 1.004006384824945863e+00;
535 break;
536 }
537
538 return p0*pow(_SC,4) +
539 p1*pow(_SC,3)+
540 p2*pow(_SC,2)+
541 p3*pow(_SC,1)+
542 p4;
543 }
544
545 int Parameter::get_high_and_compute() {
546 }
547
548 double Parameter::errorfunc(double _G1, double _GC, const matrix<double, 4, 1> &Target) {
549 double kT60 = Target(0),
550 kED = Target(1),
551 kC = Target(2),
552 kTC = Target(3);
553
554 double G = G_G1_C(_G1, _GC, kC);
555 double D1 = D1_G1_T60_C(_G1, _GC, kT60, kC);
556 double DA = DA_D1_ED(D1, kED, 0.1);
557 double T60est = compute_T60(D1, _G1, _GC, G);
558 double EDest = compute_ED(D1, DA);
559 double Cest = compute_C(_G1, _GC, G);
560 double TCest = compute_TC(D1, _G1, DA);
561
562 matrix <double, 4, 1> est_vec, orig_vec, err_vec;
563
564 float w0 = 1.0, w1=1.0, w2=1.0, w3=1.0; // Weights for errors.
565
566 est_vec(0) = w0*_N(T60est, T60_MIN, T60_MAX);
567 est_vec(1) = w1*_N(EDest, ED_MIN, ED_MAX);
568 est_vec(2) = w2*_N(Cest, C_MIN, C_MAX);
569 est_vec(3) = w3*_N(TCest, TC_MIN, TC_MAX);
570
571 orig_vec(0) = w0*_N(Target(0), T60_MIN, T60_MAX);
572 orig_vec(1) = w1*_N(Target(1), ED_MIN, ED_MAX);
573 orig_vec(2) = w2*_N(Target(2), C_MIN, C_MAX);
574 orig_vec(3) = w3*_N(Target(3), TC_MIN, TC_MAX);
575
576 err_vec = dlib::abs(est_vec-orig_vec);
577
578 return mean(pow(err_vec,2)) + variance(pow(err_vec,2));
579
580 }
581
582 matrix<double, 5, 1> Parameter::inverse_mapping(double _T60, double _ED,
583 double _C, double _TC, double _SC, double SampleRate, double _G1) {
584 /*
585 * _G1 is a starting point for the search of G1
586 */
587
588 matrix<double, 5, 1> result;
589 for (int i=0; i<5; i++) result(i) = 0.0;
590
591 double GC = GC_SC(_SC, SampleRate);
592 result(3) = GC;
593
594 matrix<double,1,1> starting_point;
595 matrix<double,4,1> Target;
596
597 Target(0) = _T60;
598 Target(1) = _ED;
599 Target(2) = _C;
600 Target(3) = _TC;
601
602 starting_point = _G1; // Starting point is current G1, for faster convergence
603
604 // Lambda function for objective function to minimize.
605
606 auto func = [&] (const matrix<double, 0, 1> &x) { float g1 = x(0); return errorfunc(g1, GC, Target); };
607
608 find_min_box_constrained(bfgs_search_strategy(),
609 objective_delta_stop_strategy(1e-9),
610 func, derivative(func), starting_point, G1_MIN, G1_MAX);
611
612 double G1 = starting_point(0);
613
614 result(2) = G1;
615
616 double G = G_G1_C(G1,GC,_C);
617
618
619 result(4) = G;
620
621 double D1 = D1_G1_T60_C(G1, GC, _T60, _C);
622
623
624 result(0) = D1;
625
626 double DA = DA_D1_ED(D1, _ED);
627 result(1) = DA;
628
629 return result;
630 }