e@1: /*
e@1:  * Parameter.cpp
e@1:  *
e@1:  *  Created on: Dec 14, 2016
e@1:  *      Author: mmxgn
e@1:  */
e@1: 
e@1: #include <cmath>
e@1: #include "Parameter.h"
e@1: //#include "csv.h"
e@1: 
e@1: 
e@1: /*
e@1:  * Given a k, return g_k and d_k based on
e@1:  * g_1 and d_1 respectively.
e@1:  *
e@1:  */
e@1: 
e@1: #define GK(k) powf(G1, powf(1.5, 1-k))
e@1: #define _GK(k) powf(_G1, powf(1.5, 1-k))
e@1: 
e@1: #define DK(k) (D1*powf(1.5, 1-k))
e@1: #define _DK(k) (_D1*powf(1.5, 1-k))
e@1: 
e@1: /*
e@1:  * The parameter class
e@1:  */
e@1: 
e@1: Parameter::Parameter(double _SR,
e@1: 					 double _D1,
e@1: 					 double _DA,
e@1: 					 double _G1,
e@1: 					 double _GC,
e@1: 					 double _G,
e@1: 					 double _T60,
e@1: 					 double _ED,
e@1: 					 double _C,
e@1: 					 double _TC,
e@1: 					 double _SC)
e@1: :
e@1: 	D1(_D1),
e@1: 	DA(_DA),
e@1: 	G1(_G1),
e@1: 	GC(_GC),
e@1: 	G(_G),
e@1: 	T60(_T60),
e@1: 	ED(_ED),
e@1: 	C(_C),
e@1: 	TC(_TC),
e@1: 	SC(_SC),
e@1: 	SC_MAX(_SR/4.0),
e@1: 	SampleRate(_SR)
e@1: {
e@1: 	// pass
e@1: }
e@1: 
e@1: double Parameter::getC() const {
e@1: 	return C;
e@1: }
e@1: 
e@1: void Parameter::setC(double c) {
e@1: 	C = c;
e@1: }
e@1: 
e@1: double Parameter::getD1() const {
e@1: 	return D1;
e@1: }
e@1: 
e@1: void Parameter::setD1(double d1) {
e@1: 	/* Dependency for:
e@1: 	 * 	T60, ED, TC
e@1: 	 */
e@1: 
e@1: 	T60 = compute_T60(d1, G1, GC, G);
e@1: 	ED = compute_ED(d1, DA);
e@1: 	TC = compute_TC(d1, G1, DA);
e@1: 
e@1: }
e@1: 
e@1: double Parameter::getDa() const {
e@1: 	return DA;
e@1: 
e@1: }
e@1: 
e@1: void Parameter::setDa(double da) {
e@1: 	/* Dependency for:
e@1: 	 * 	ED, TC
e@1: 	 */
e@1: 
e@1: 	ED = compute_ED(D1, da);
e@1: 	TC = compute_TC(D1, G1, da);
e@1: 
e@1: }
e@1: 
e@1: double Parameter::getEd() const {
e@1: 	return ED;
e@1: }
e@1: 
e@1: void Parameter::setEd(double ed) {
e@1: 	matrix <double,5,1> est_param = inverse_mapping(T60, ed, C, TC, SC, SampleRate,G1);
e@1: 	D1 = est_param(0);
e@1: 	DA = est_param(1);
e@1: 	G1 = est_param(2);
e@1: 	GC = est_param(3);
e@1: 	G = est_param(4);
e@1: 	T60 = compute_T60(D1, G1, GC, G);
e@1: 	ED = compute_ED(D1, DA);
e@1: 	C = compute_C(G1, GC, G);
e@1: 	TC = compute_TC(D1, G1, DA);
e@1: 	SC = compute_SC(GC, SampleRate);
e@1: }
e@1: 
e@1: double Parameter::getG() const {
e@1: 	return G;
e@1: }
e@1: 
e@1: void Parameter::setG(double g) {
e@1: 	/*
e@1: 	 * Dependency for:
e@1: 	 * 	T60, C
e@1: 	 */
e@1: 
e@1: 	T60 = compute_T60(D1, G1, GC, g);
e@1: 	C = compute_C(G1, GC, g);
e@1: 
e@1: }
e@1: 
e@1: double Parameter::getG1() const {
e@1: 	return G1;
e@1: }
e@1: 
e@1: void Parameter::setG1(double g1) {
e@1: 	/*
e@1: 	 * Dependency for:
e@1: 	 * 	T60, C, TC
e@1: 	 */
e@1: 
e@1: 	T60 = compute_T60(D1, g1, GC, G);
e@1: 	C = compute_C(g1, GC, G);
e@1: 	TC = compute_TC(D1, g1, DA);
e@1: 
e@1: }
e@1: 
e@1: double Parameter::getGc() const {
e@1: 	return GC;
e@1: }
e@1: 
e@1: void Parameter::setGc(double gc) {
e@1: 	/*
e@1: 	 * Dependency for:
e@1: 	 * 	SC
e@1: 	 */
e@1: 
e@1: 	SC = compute_SC(gc, SampleRate);
e@1: }
e@1: 
e@1: double Parameter::getSampleRate() const {
e@1: 	return SampleRate;
e@1: }
e@1: 
e@1: void Parameter::setSampleRate(double sampleRate) {
e@1: 	SampleRate = sampleRate;
e@1: 	set_low_and_compute(D1, DA, G1, GC, G);
e@1: 	matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, TC, SC, SampleRate,G1);
e@1: 	D1 = est_param(0);
e@1: 	DA = est_param(1);
e@1: 	G1 = est_param(2);
e@1: 	GC = est_param(3);
e@1: 	G = est_param(4);
e@1: 	T60 = compute_T60(D1, G1, GC, G);
e@1: 	ED = compute_ED(D1, DA);
e@1: 	C = compute_C(G1, GC, G);
e@1: 	TC = compute_TC(D1, G1, DA);
e@1: 	SC = compute_SC(GC, SampleRate);
e@1: }
e@1: 
e@1: double Parameter::getSc() const {
e@1: 	return SC;
e@1: }
e@1: 
e@1: void Parameter::setSc(double sc) {
e@1: 	matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, TC, sc, SampleRate,G1);
e@1: 	D1 = est_param(0);
e@1: 	DA = est_param(1);
e@1: 	G1 = est_param(2);
e@1: 	GC = est_param(3);
e@1: 	G = est_param(4);
e@1: 	T60 = compute_T60(D1, G1, GC, G);
e@1: 	ED = compute_ED(D1, DA);
e@1: 	C = compute_C(G1, GC, G);
e@1: 	TC = compute_TC(D1, G1, DA);
e@1: 	SC = compute_SC(GC, SampleRate);
e@1: }
e@1: 
e@1: double Parameter::getT60() const {
e@1: 	return T60;
e@1: }
e@1: 
e@1: void Parameter::setT60(double t60) {
e@1: 	matrix <double,5,1> est_param = inverse_mapping(t60, ED, C, TC, SC, SampleRate,G1);
e@1: 	D1 = est_param(0);
e@1: 	DA = est_param(1);
e@1: 	G1 = est_param(2);
e@1: 	GC = est_param(3);
e@1: 	G = est_param(4);
e@1: 	T60 = compute_T60(D1, G1, GC, G);
e@1: 	ED = compute_ED(D1, DA);
e@1: 	C = compute_C(G1, GC, G);
e@1: 	TC = compute_TC(D1, G1, DA);
e@1: 	SC = compute_SC(GC, SampleRate);
e@1: }
e@1: 
e@1: double Parameter::getTc() const {
e@1: 	return TC;
e@1: }
e@1: 
e@1: void Parameter::setTc(double tc) {
e@1: 	matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, tc, SC, SampleRate,G1);
e@1: 	D1 = est_param(0);
e@1: 	DA = est_param(1);
e@1: 	G1 = est_param(2);
e@1: 	GC = est_param(3);
e@1: 	G = est_param(4);
e@1: 	T60 = compute_T60(D1, G1, GC, G);
e@1: 	ED = compute_ED(D1, DA);
e@1: 	C = compute_C(G1, GC, G);
e@1: 	TC = compute_TC(D1, G1, DA);
e@1: 	SC = compute_SC(GC, SampleRate);
e@1: }
e@1: 
e@1: Parameter::~Parameter() {
e@1: 	// TODO Auto-generated destructor stub
e@1: }
e@1: 
e@1: std::ostream& operator<<(std::ostream &s, const Parameter &p) {
e@1:   return s << "Parameter({High:<" <<
e@1: 		  p.T60 << ", " <<
e@1: 		  p.ED << ", " <<
e@1: 		  p.C << ", " <<
e@1: 		  p.TC << ", " <<
e@1: 		  p.SC << ">}, " <<
e@1: 		  "{Low: <" <<
e@1: 		  p.D1 << ", " <<
e@1: 		  p.DA << ", " <<
e@1: 		  p.G1 << ", " <<
e@1: 		  p.GC << ", " <<
e@1: 		  p.G  << ">}).";
e@1: 
e@1: }
e@1: 
e@1: int Parameter::set_low_and_compute(double D1, double DA, double G1, double GC,
e@1: 		double G) {
e@1: 
e@1: 	/* Temporary values */
e@1: 
e@1: 	double _T60, _ED, _C, _TC, _SC;
e@1: 
e@1: 	_T60 = compute_T60(D1, G1, GC, G);
e@1: 	if (_T60 < T60_MIN || _T60 > T60_MAX)
e@1: 		return -1;
e@1: 
e@1: 	_ED = compute_ED(D1, DA);
e@1: 	if (_ED < ED_MIN || _ED > ED_MAX)
e@1: 		return -1;
e@1: 
e@1: 	_C = compute_C(G1, GC, G);
e@1: 
e@1: 	if (_C < C_MIN || _C > C_MAX)
e@1: 		return -1;
e@1: 
e@1: 	_TC = compute_TC(D1, G1, DA);
e@1: 
e@1: 	if (_TC < TC_MIN || _TC > TC_MAX)
e@1: 		return -1;
e@1: 
e@1: 	_SC = compute_SC(GC, SampleRate);
e@1: 
e@1: 	if (_SC < SC_MIN || _SC > SC_MAX)
e@1: 		return -1;
e@1: 
e@1: 	T60 = _T60;
e@1: 	ED = _ED;
e@1: 	C = _C;
e@1: 	TC = _TC;
e@1: 	SC = _SC;
e@1: 
e@1: 	return 0;
e@1: 
e@1: }
e@1: 
e@1: double Parameter::compute_T60(double _D1, double _G1, double _GC, double _G) {
e@1: 	return -_D1/log(_G1)*(log(_G)+log(1-_GC)+log(GA)+3*log(10));
e@1: }
e@1: 
e@1: double Parameter::compute_ED(double _D1, double _DA) {
e@1: 	return 2.078125/_D1/_DA;
e@1: }
e@1: 
e@1: double Parameter::compute_C(double _G1, double _GC, double _G) {
e@1: 	double sum1 = 0;
e@1: 	for (int k=1; k<7; k++){
e@1: 		sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1, 2*powf(1.5, 1-k)));
e@1: 	}
e@1: 	return -10*log10(_G*_G*(1-_GC)/(1+_GC)*sum1);
e@1: }
e@1: 
e@1: double Parameter::compute_TC(double _D1, double _G1, double _DA) {
e@1: 	double sum1 = 0, sum2 = 0, s1 = 0;
e@1: 	for (int k=1; k<7; k++){
e@1: 		s1 = powf(_GK(k), 2.0)/(1-powf(_GK(k), 2.0));
e@1: 		sum2 += s1;
e@1: 		sum1 += _DK(k)*s1/(1-powf(_GK(k),2.0));
e@1: 	}
e@1: 
e@1: 	return sum1/sum2 + _DA;
e@1: }
e@1: 
e@1: double Parameter::compute_SC(double _GC, double _SR) {
e@1: 	double sum1 = 0, sum2 = 0;
e@1: 
e@1: 	for (int k=0; k < int(_SR/2.0); k++){
e@1: 		double lut = 1.0-2.0*_GC*cos(2.0*M_PI*k/_SR)+_GC*_GC;
e@1: 		sum1 += k/lut;
e@1: 		sum2 += 1/lut;
e@1: 	}
e@1: 
e@1: 	return sum1/sum2;
e@1: }
e@1: 
e@1: int Parameter::get_low_and_computer() {
e@1: 	return 0;
e@1: }
e@1: 
e@1: int Parameter::set_high_and_compute(double T60, double ED, double C, double TC,
e@1: 		double SC) {
e@1: 	return 0;
e@1: }
e@1: 
e@1: double Parameter::TC_G1_D1_ED(double _D1, double _G1, double _ED, double t) {
e@1: 	double sum1 = 0, sum2 = 0;
e@1: 
e@1: 	for (int k=1; k<7; k++){
e@1: 		sum1 += powf(1.5, -k+1)*_D1*
e@1: 				powf(_G1, 2*powf(1.5, -k+1))/
e@1: 				powf((-powf(_G1, 2*powf(1.5, -k+1))+1),2);
e@1: 		sum2 += powf(_G1,2*powf(1.5,-k+1))/
e@1: 				(-powf(_G1, 2*powf(1.5,-k+1))+1);
e@1: 	}
e@1: 
e@1: 	return sum1/sum2 + 20.78125*t/(_ED*_D1);
e@1: }
e@1: 
e@1: double Parameter::T60_D1_G1_C(double _D1, double _G1, double _GC, double _C) {
e@1: 	double sum1 = 0;
e@1: 
e@1: 	for (int k=1; k<7; k++){
e@1: 		sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1,2*powf(1.5,1-k)));
e@1: 	}
e@1: 
e@1: 	return _D1*log(0.001/GA*exp(_C*log(10)/20)/
e@1: 			(sqrt((1+_GC)/((1-_GC)*sum1))*(1-_GC)))/log(_G1);
e@1: }
e@1: 
e@1: double Parameter::D1_G1_T60_C(double _G1, double _GC, double _T60, double _C) {
e@1: 	double sum1 = 0;
e@1: 
e@1: 	for (int k=1; k<7; k++){
e@1: 		sum1 += powf(_G1,(2*powf(1.5,1-k)))/(1-powf(_G1,2*powf(1.5,1-k)));
e@1: 	}
e@1: 
e@1: 	return _T60*log(_G1)/log(0.001/GA*exp(_C*log(10)/20)/(sqrt((1+_GC)/((1-_GC)*sum1))*(1-_GC)));
e@1: }
e@1: 
e@1: double Parameter::G_G1_C(double _G1, double _GC, double _C) {
e@1: 	double sum1 = 0;
e@1: 	for (int k=1; k<7; k++){
e@1: 		sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1, 2*powf(1.5, 1-k)));
e@1: 	}
e@1: 
e@1: 	return sqrt(exp(-_C/10*log(10))/(1-_GC)*(1+_GC)/sum1);
e@1: }
e@1: 
e@1: double Parameter::DA_D1_ED(double _D1, double _ED, double t) {
e@1: 	return 20.78125*t/_D1/_ED;
e@1: }
e@1: 
e@1: 
e@1: double Parameter::GC_SC(double _SC, double SR) {
e@1: 	double p0 = 1.0,
e@1: 			p1 = 1.0,
e@1: 			p2 = 1.0,
e@1: 			p3 = 1.0,
e@1: 			p4 = 1.0;
e@1: 
e@1: 	switch (int(SR)){
e@1: 	case 8000:
e@1: 		p0 = -7.545401867799958007e-15;
e@1: 		p1 = 6.605438530163032344e-11;
e@1: 		p2 = -2.339787304666883204e-07;
e@1: 		p3 = -2.380996920823082359e-04;
e@1: 		p4 = 1.003686549974644349e+00;
e@1: 	break;
e@1: 	case 11025:
e@1: 		p0 = -2.090495384250846006e-15;
e@1: 		p1 = 2.523110858472787412e-11;
e@1: 		p2 = -1.232004017234833354e-07;
e@1: 		p3 = -1.727859517219043628e-04;
e@1: 		p4 = 1.003773048218970398e+00;
e@1: 	break;
e@1: 	case 16000:
e@1: 		p0 = -4.707928872668820147e-16;
e@1: 		p1 = 8.250344830381294683e-12;
e@1: 		p2 = -5.847797160444276075e-08;
e@1: 		p3 = -1.190647306137786752e-04;
e@1: 		p4 = 1.003845141060948354e+00;
e@1: 	break;
e@1: 	case 21000:
e@1: 		p0 = -1.585759862171924958e-16;
e@1: 		p1 = 3.648228536317020713e-12;
e@1: 		p2 = -3.394367375513113085e-08;
e@1: 		p3 = -9.071942141117647239e-05;
e@1: 		p4 = 1.003883355464454619e+00;
e@1: 	break;
e@1: 	case 22050:
e@1: 		p0 = -1.304517046624577766e-16;
e@1: 		p1 = 3.151372288676199373e-12;
e@1: 		p2 = -3.078752354037690531e-08;
e@1: 		p3 = -8.639997174340922025e-05;
e@1: 		p4 = 1.003889194220798364e+00;
e@1: 	break;
e@1: 	case 32000:
e@1: 		p0 = -2.939615050414309854e-17;
e@1: 		p1 = 1.030821612883252373e-12;
e@1: 		p2 = -1.461695000370224873e-08;
e@1: 		p3 = -5.953732000720860076e-05;
e@1: 		p4 = 1.003925597252757207e+00;
e@1: 	break;
e@1: 	case 37800:
e@1: 		p0 = -1.509578131570550335e-17;
e@1: 		p1 = 6.253538622497773237e-13;
e@1: 		p2 = -1.047516110428303656e-08;
e@1: 		p3 = -5.040267062802067523e-05;
e@1: 		p4 = 1.003938012276286873e+00;
e@1: 	break;
e@1: 	case 44056:
e@1: 		p0 = -8.179918755257568798e-18;
e@1: 		p1 = 3.949666985075376660e-13;
e@1: 		p2 = -7.711241596699441441e-09;
e@1: 		p3 = -4.324592160336278773e-05;
e@1: 		p4 = 1.003947751979878067e+00;
e@1: 	break;
e@1: 	case 44100:
e@1: 		p0 = -8.147316012596786372e-18;
e@1: 		p1 = 3.937855238616550567e-13;
e@1: 		p2 = -7.695860699409773702e-09;
e@1: 		p3 = -4.320277669373632828e-05;
e@1: 		p4 = 1.003947810730627443e+00;
e@1: 	break;
e@1: 	case 47250:
e@1: 		p0 = -6.182156102335834534e-18;
e@1: 		p1 = 3.201545709180980074e-13;
e@1: 		p2 = -6.703887839395126650e-09;
e@1: 		p3 = -4.032277597168298445e-05;
e@1: 		p4 = 1.003951733380953559e+00;
e@1: 	break;
e@1: 	case 48000:
e@1: 		p0 = -5.804669671441499519e-18;
e@1: 		p1 = 3.053789928757690291e-13;
e@1: 		p2 = -6.496014866101249794e-09;
e@1: 		p3 = -3.969277247952757479e-05;
e@1: 		p4 = 1.003952591709977415e+00;
e@1: 	break;
e@1: 	case 50000:
e@1: 		p0 = -4.930039877651698485e-18;
e@1: 		p1 = 2.701762164752109105e-13;
e@1: 		p2 = -5.986696651033372113e-09;
e@1: 		p3 = -3.810515833478767255e-05;
e@1: 		p4 = 1.003954755096216900e+00;
e@1: 	break;
e@1: 	case 50400:
e@1: 		p0 = -4.775358854140600092e-18;
e@1: 		p1 = 2.637937027466497354e-13;
e@1: 		p2 = -5.892041055138477995e-09;
e@1: 		p3 = -3.780275477003837485e-05;
e@1: 		p4 = 1.003955167234265478e+00;
e@1: 	break;
e@1: 	case 88200:
e@1: 		p0 = -5.090139768480031043e-19;
e@1: 		p1 = 4.921424693356908116e-14;
e@1: 		p2 = -1.923829523675108346e-09;
e@1: 		p3 = -2.160214687014651151e-05;
e@1: 		p4 = 1.003977276279834019e+00;
e@1: 	break;
e@1: 	case 96000:
e@1: 		p0 = -3.626648836744412987e-19;
e@1: 		p1 = 3.816597595693573777e-14;
e@1: 		p2 = -1.623898000952206408e-09;
e@1: 		p3 = -1.984703070895103917e-05;
e@1: 		p4 = 1.003979674995812754e+00;
e@1: 	break;
e@1: 	case 176400:
e@1: 		p0 = -3.180720104052761169e-20;
e@1: 		p1 = 6.151207927460715513e-15;
e@1: 		p2 = -4.809399117661735116e-10;
e@1: 		p3 = -1.080127073758081756e-05;
e@1: 		p4 = 1.003992048707243789e+00;
e@1: 	break;
e@1: 	case 192000:
e@1: 		p0 = -2.266250752087159745e-20;
e@1: 		p1 = 4.770337959813220317e-15;
e@1: 		p2 = -4.059609201723455372e-10;
e@1: 		p3 = -9.923682425709201677e-06;
e@1: 		p4 = 1.003993250132436454e+00;
e@1: 	break;
e@1: 	case 352800:
e@1: 		p0 = -1.987755090799793513e-21;
e@1: 		p1 = 7.688647436953224411e-16;
e@1: 		p2 = -1.202327621003883834e-10;
e@1: 		p3 = -5.400685657755166908e-06;
e@1: 		p4 = 1.003999444875642189e+00;
e@1: 	break;
e@1: 	case 2822400:
e@1: 		p0 = -4.852498403591686415e-25;
e@1: 		p1 = 1.501626536939688520e-18;
e@1: 		p2 = -1.878606337281021417e-12;
e@1: 		p3 = -6.750912735944337505e-07;
e@1: 		p4 = 1.004005921979573834e+00;
e@1: 	break;
e@1: 	case 5644800:
e@1: 		p0 = -3.032792721276874507e-26;
e@1: 		p1 = 1.877027577372803603e-19;
e@1: 		p2 = -4.696510359527020341e-13;
e@1: 		p3 = -3.375458367765195566e-07;
e@1: 		p4 = 1.004006384824945863e+00;
e@1: 	break;
e@1: 	}
e@1: 
e@1: 	return p0*pow(_SC,4) +
e@1: 			p1*pow(_SC,3)+
e@1: 			p2*pow(_SC,2)+
e@1: 			p3*pow(_SC,1)+
e@1: 			p4;
e@1: }
e@1: 
e@1: int Parameter::get_high_and_compute() {
e@1: }
e@1: 
e@1: double Parameter::errorfunc(double _G1, double _GC, const matrix<double, 4, 1> &Target) {
e@1: 	double kT60 = Target(0),
e@1: 			kED = Target(1),
e@1: 			kC  = Target(2),
e@1: 			kTC = Target(3);
e@1: 
e@1: 	double G = G_G1_C(_G1, _GC, kC);
e@1: 	double D1 = D1_G1_T60_C(_G1, _GC, kT60, kC);
e@1: 	double DA = DA_D1_ED(D1, kED, 0.1);
e@1: 	double T60est = compute_T60(D1, _G1, _GC, G);
e@1: 	double EDest = compute_ED(D1, DA);
e@1: 	double Cest = compute_C(_G1, _GC, G);
e@1: 	double TCest = compute_TC(D1, _G1, DA);
e@1: 
e@1: 	matrix <double, 4, 1> est_vec, orig_vec, err_vec;
e@1: 
e@1: 	float w0 = 1.0, w1=1.0, w2=1.0, w3=1.0; // Weights for errors.
e@1: 
e@1: 	est_vec(0) = w0*_N(T60est, T60_MIN, T60_MAX);
e@1: 	est_vec(1) = w1*_N(EDest, ED_MIN, ED_MAX);
e@1: 	est_vec(2) = w2*_N(Cest, C_MIN, C_MAX);
e@1: 	est_vec(3) = w3*_N(TCest, TC_MIN, TC_MAX);
e@1: 
e@1: 	orig_vec(0) = w0*_N(Target(0), T60_MIN, T60_MAX);
e@1: 	orig_vec(1) = w1*_N(Target(1), ED_MIN, ED_MAX);
e@1: 	orig_vec(2) = w2*_N(Target(2), C_MIN, C_MAX);
e@1: 	orig_vec(3) = w3*_N(Target(3), TC_MIN, TC_MAX);
e@1: 
e@1: 	err_vec = dlib::abs(est_vec-orig_vec);
e@1: 
e@1: 	return mean(pow(err_vec,2)) + variance(pow(err_vec,2));
e@1: 
e@1: }
e@1: 
e@1: matrix<double, 5, 1> Parameter::inverse_mapping(double _T60, double _ED,
e@1: 		double _C, double _TC, double _SC, double SampleRate, double _G1) {
e@1:     /*
e@1:      * _G1 is a starting point for the search of G1
e@1:      */
e@1: 
e@1: 	matrix<double, 5, 1> result;
e@1: 	for (int i=0; i<5; i++) result(i) = 0.0;
e@1: 
e@1: 	double GC = GC_SC(_SC, SampleRate);
e@1: 	result(3) = GC;
e@1: 
e@1: 	matrix<double,1,1> starting_point;
e@1: 	matrix<double,4,1> Target;
e@1: 
e@1: 	Target(0) = _T60;
e@1: 	Target(1) = _ED;
e@1: 	Target(2) = _C;
e@1: 	Target(3) = _TC;
e@1: 
e@1: 	starting_point = _G1; // Starting point is current G1, for faster convergence
e@1: 
e@1: 	// Lambda function for objective function to minimize.
e@1: 
e@1: 	auto func = [&] (const matrix<double, 0, 1> &x) { float g1 = x(0); return errorfunc(g1, GC, Target); };
e@1: 
e@1:     find_min_box_constrained(bfgs_search_strategy(),
e@1:                              objective_delta_stop_strategy(1e-9),
e@1: 							 func, derivative(func), starting_point, G1_MIN, G1_MAX);
e@1: 
e@1:     double G1 = starting_point(0);
e@1: 
e@1:     result(2) = G1;
e@1: 
e@1:     double G = G_G1_C(G1,GC,_C);
e@1: 
e@1: 
e@1:     result(4) = G;
e@1: 
e@1:     double D1 = D1_G1_T60_C(G1, GC, _T60, _C);
e@1: 
e@1: 
e@1:     result(0) = D1;
e@1: 
e@1:     double DA = DA_D1_ED(D1, _ED);
e@1:     result(1) = DA;
e@1: 
e@1: 	return result;
e@1: }