view cpp/chourdakisreiss2016cpp/Parameter.cpp @ 2:c87a9505f294 tip

Added LICENSE for code, removed .wav files

/*
 * Parameter.cpp
 *
 *  Created on: Dec 14, 2016
 *      Author: mmxgn
 */

#include <cmath>
#include "Parameter.h"
//#include "csv.h"


/*
 * Given a k, return g_k and d_k based on
 * g_1 and d_1 respectively.
 *
 */

#define GK(k) powf(G1, powf(1.5, 1-k))
#define _GK(k) powf(_G1, powf(1.5, 1-k))

#define DK(k) (D1*powf(1.5, 1-k))
#define _DK(k) (_D1*powf(1.5, 1-k))

/*
 * The parameter class
 */

Parameter::Parameter(double _SR,
					 double _D1,
					 double _DA,
					 double _G1,
					 double _GC,
					 double _G,
					 double _T60,
					 double _ED,
					 double _C,
					 double _TC,
					 double _SC)
:
	D1(_D1),
	DA(_DA),
	G1(_G1),
	GC(_GC),
	G(_G),
	T60(_T60),
	ED(_ED),
	C(_C),
	TC(_TC),
	SC(_SC),
	SC_MAX(_SR/4.0),
	SampleRate(_SR)
{
	// pass
}

double Parameter::getC() const {
	return C;
}

void Parameter::setC(double c) {
	C = c;
}

double Parameter::getD1() const {
	return D1;
}

void Parameter::setD1(double d1) {
	/* Dependency for:
	 * 	T60, ED, TC
	 */

	T60 = compute_T60(d1, G1, GC, G);
	ED = compute_ED(d1, DA);
	TC = compute_TC(d1, G1, DA);

}

double Parameter::getDa() const {
	return DA;

}

void Parameter::setDa(double da) {
	/* Dependency for:
	 * 	ED, TC
	 */

	ED = compute_ED(D1, da);
	TC = compute_TC(D1, G1, da);

}

double Parameter::getEd() const {
	return ED;
}

void Parameter::setEd(double ed) {
	matrix <double,5,1> est_param = inverse_mapping(T60, ed, C, TC, SC, SampleRate,G1);
	D1 = est_param(0);
	DA = est_param(1);
	G1 = est_param(2);
	GC = est_param(3);
	G = est_param(4);
	T60 = compute_T60(D1, G1, GC, G);
	ED = compute_ED(D1, DA);
	C = compute_C(G1, GC, G);
	TC = compute_TC(D1, G1, DA);
	SC = compute_SC(GC, SampleRate);
}

double Parameter::getG() const {
	return G;
}

void Parameter::setG(double g) {
	/*
	 * Dependency for:
	 * 	T60, C
	 */

	T60 = compute_T60(D1, G1, GC, g);
	C = compute_C(G1, GC, g);

}

double Parameter::getG1() const {
	return G1;
}

void Parameter::setG1(double g1) {
	/*
	 * Dependency for:
	 * 	T60, C, TC
	 */

	T60 = compute_T60(D1, g1, GC, G);
	C = compute_C(g1, GC, G);
	TC = compute_TC(D1, g1, DA);

}

double Parameter::getGc() const {
	return GC;
}

void Parameter::setGc(double gc) {
	/*
	 * Dependency for:
	 * 	SC
	 */

	SC = compute_SC(gc, SampleRate);
}

double Parameter::getSampleRate() const {
	return SampleRate;
}

void Parameter::setSampleRate(double sampleRate) {
	SampleRate = sampleRate;
	set_low_and_compute(D1, DA, G1, GC, G);
	matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, TC, SC, SampleRate,G1);
	D1 = est_param(0);
	DA = est_param(1);
	G1 = est_param(2);
	GC = est_param(3);
	G = est_param(4);
	T60 = compute_T60(D1, G1, GC, G);
	ED = compute_ED(D1, DA);
	C = compute_C(G1, GC, G);
	TC = compute_TC(D1, G1, DA);
	SC = compute_SC(GC, SampleRate);
}

double Parameter::getSc() const {
	return SC;
}

void Parameter::setSc(double sc) {
	matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, TC, sc, SampleRate,G1);
	D1 = est_param(0);
	DA = est_param(1);
	G1 = est_param(2);
	GC = est_param(3);
	G = est_param(4);
	T60 = compute_T60(D1, G1, GC, G);
	ED = compute_ED(D1, DA);
	C = compute_C(G1, GC, G);
	TC = compute_TC(D1, G1, DA);
	SC = compute_SC(GC, SampleRate);
}

double Parameter::getT60() const {
	return T60;
}

void Parameter::setT60(double t60) {
	matrix <double,5,1> est_param = inverse_mapping(t60, ED, C, TC, SC, SampleRate,G1);
	D1 = est_param(0);
	DA = est_param(1);
	G1 = est_param(2);
	GC = est_param(3);
	G = est_param(4);
	T60 = compute_T60(D1, G1, GC, G);
	ED = compute_ED(D1, DA);
	C = compute_C(G1, GC, G);
	TC = compute_TC(D1, G1, DA);
	SC = compute_SC(GC, SampleRate);
}

double Parameter::getTc() const {
	return TC;
}

void Parameter::setTc(double tc) {
	matrix <double,5,1> est_param = inverse_mapping(T60, ED, C, tc, SC, SampleRate,G1);
	D1 = est_param(0);
	DA = est_param(1);
	G1 = est_param(2);
	GC = est_param(3);
	G = est_param(4);
	T60 = compute_T60(D1, G1, GC, G);
	ED = compute_ED(D1, DA);
	C = compute_C(G1, GC, G);
	TC = compute_TC(D1, G1, DA);
	SC = compute_SC(GC, SampleRate);
}

Parameter::~Parameter() {
	// TODO Auto-generated destructor stub
}

std::ostream& operator<<(std::ostream &s, const Parameter &p) {
  return s << "Parameter({High:<" <<
		  p.T60 << ", " <<
		  p.ED << ", " <<
		  p.C << ", " <<
		  p.TC << ", " <<
		  p.SC << ">}, " <<
		  "{Low: <" <<
		  p.D1 << ", " <<
		  p.DA << ", " <<
		  p.G1 << ", " <<
		  p.GC << ", " <<
		  p.G  << ">}).";

}

int Parameter::set_low_and_compute(double D1, double DA, double G1, double GC,
		double G) {

	/* Temporary values */

	double _T60, _ED, _C, _TC, _SC;

	_T60 = compute_T60(D1, G1, GC, G);
	if (_T60 < T60_MIN || _T60 > T60_MAX)
		return -1;

	_ED = compute_ED(D1, DA);
	if (_ED < ED_MIN || _ED > ED_MAX)
		return -1;

	_C = compute_C(G1, GC, G);

	if (_C < C_MIN || _C > C_MAX)
		return -1;

	_TC = compute_TC(D1, G1, DA);

	if (_TC < TC_MIN || _TC > TC_MAX)
		return -1;

	_SC = compute_SC(GC, SampleRate);

	if (_SC < SC_MIN || _SC > SC_MAX)
		return -1;

	T60 = _T60;
	ED = _ED;
	C = _C;
	TC = _TC;
	SC = _SC;

	return 0;

}

double Parameter::compute_T60(double _D1, double _G1, double _GC, double _G) {
	return -_D1/log(_G1)*(log(_G)+log(1-_GC)+log(GA)+3*log(10));
}

double Parameter::compute_ED(double _D1, double _DA) {
	return 2.078125/_D1/_DA;
}

double Parameter::compute_C(double _G1, double _GC, double _G) {
	double sum1 = 0;
	for (int k=1; k<7; k++){
		sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1, 2*powf(1.5, 1-k)));
	}
	return -10*log10(_G*_G*(1-_GC)/(1+_GC)*sum1);
}

double Parameter::compute_TC(double _D1, double _G1, double _DA) {
	double sum1 = 0, sum2 = 0, s1 = 0;
	for (int k=1; k<7; k++){
		s1 = powf(_GK(k), 2.0)/(1-powf(_GK(k), 2.0));
		sum2 += s1;
		sum1 += _DK(k)*s1/(1-powf(_GK(k),2.0));
	}

	return sum1/sum2 + _DA;
}

double Parameter::compute_SC(double _GC, double _SR) {
	double sum1 = 0, sum2 = 0;

	for (int k=0; k < int(_SR/2.0); k++){
		double lut = 1.0-2.0*_GC*cos(2.0*M_PI*k/_SR)+_GC*_GC;
		sum1 += k/lut;
		sum2 += 1/lut;
	}

	return sum1/sum2;
}

int Parameter::get_low_and_computer() {
	return 0;
}

int Parameter::set_high_and_compute(double T60, double ED, double C, double TC,
		double SC) {
	return 0;
}

double Parameter::TC_G1_D1_ED(double _D1, double _G1, double _ED, double t) {
	double sum1 = 0, sum2 = 0;

	for (int k=1; k<7; k++){
		sum1 += powf(1.5, -k+1)*_D1*
				powf(_G1, 2*powf(1.5, -k+1))/
				powf((-powf(_G1, 2*powf(1.5, -k+1))+1),2);
		sum2 += powf(_G1,2*powf(1.5,-k+1))/
				(-powf(_G1, 2*powf(1.5,-k+1))+1);
	}

	return sum1/sum2 + 20.78125*t/(_ED*_D1);
}

double Parameter::T60_D1_G1_C(double _D1, double _G1, double _GC, double _C) {
	double sum1 = 0;

	for (int k=1; k<7; k++){
		sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1,2*powf(1.5,1-k)));
	}

	return _D1*log(0.001/GA*exp(_C*log(10)/20)/
			(sqrt((1+_GC)/((1-_GC)*sum1))*(1-_GC)))/log(_G1);
}

double Parameter::D1_G1_T60_C(double _G1, double _GC, double _T60, double _C) {
	double sum1 = 0;

	for (int k=1; k<7; k++){
		sum1 += powf(_G1,(2*powf(1.5,1-k)))/(1-powf(_G1,2*powf(1.5,1-k)));
	}

	return _T60*log(_G1)/log(0.001/GA*exp(_C*log(10)/20)/(sqrt((1+_GC)/((1-_GC)*sum1))*(1-_GC)));
}

double Parameter::G_G1_C(double _G1, double _GC, double _C) {
	double sum1 = 0;
	for (int k=1; k<7; k++){
		sum1 += powf(_G1, 2*powf(1.5, 1-k))/(1-powf(_G1, 2*powf(1.5, 1-k)));
	}

	return sqrt(exp(-_C/10*log(10))/(1-_GC)*(1+_GC)/sum1);
}

double Parameter::DA_D1_ED(double _D1, double _ED, double t) {
	return 20.78125*t/_D1/_ED;
}


double Parameter::GC_SC(double _SC, double SR) {
	double p0 = 1.0,
			p1 = 1.0,
			p2 = 1.0,
			p3 = 1.0,
			p4 = 1.0;

	switch (int(SR)){
	case 8000:
		p0 = -7.545401867799958007e-15;
		p1 = 6.605438530163032344e-11;
		p2 = -2.339787304666883204e-07;
		p3 = -2.380996920823082359e-04;
		p4 = 1.003686549974644349e+00;
	break;
	case 11025:
		p0 = -2.090495384250846006e-15;
		p1 = 2.523110858472787412e-11;
		p2 = -1.232004017234833354e-07;
		p3 = -1.727859517219043628e-04;
		p4 = 1.003773048218970398e+00;
	break;
	case 16000:
		p0 = -4.707928872668820147e-16;
		p1 = 8.250344830381294683e-12;
		p2 = -5.847797160444276075e-08;
		p3 = -1.190647306137786752e-04;
		p4 = 1.003845141060948354e+00;
	break;
	case 21000:
		p0 = -1.585759862171924958e-16;
		p1 = 3.648228536317020713e-12;
		p2 = -3.394367375513113085e-08;
		p3 = -9.071942141117647239e-05;
		p4 = 1.003883355464454619e+00;
	break;
	case 22050:
		p0 = -1.304517046624577766e-16;
		p1 = 3.151372288676199373e-12;
		p2 = -3.078752354037690531e-08;
		p3 = -8.639997174340922025e-05;
		p4 = 1.003889194220798364e+00;
	break;
	case 32000:
		p0 = -2.939615050414309854e-17;
		p1 = 1.030821612883252373e-12;
		p2 = -1.461695000370224873e-08;
		p3 = -5.953732000720860076e-05;
		p4 = 1.003925597252757207e+00;
	break;
	case 37800:
		p0 = -1.509578131570550335e-17;
		p1 = 6.253538622497773237e-13;
		p2 = -1.047516110428303656e-08;
		p3 = -5.040267062802067523e-05;
		p4 = 1.003938012276286873e+00;
	break;
	case 44056:
		p0 = -8.179918755257568798e-18;
		p1 = 3.949666985075376660e-13;
		p2 = -7.711241596699441441e-09;
		p3 = -4.324592160336278773e-05;
		p4 = 1.003947751979878067e+00;
	break;
	case 44100:
		p0 = -8.147316012596786372e-18;
		p1 = 3.937855238616550567e-13;
		p2 = -7.695860699409773702e-09;
		p3 = -4.320277669373632828e-05;
		p4 = 1.003947810730627443e+00;
	break;
	case 47250:
		p0 = -6.182156102335834534e-18;
		p1 = 3.201545709180980074e-13;
		p2 = -6.703887839395126650e-09;
		p3 = -4.032277597168298445e-05;
		p4 = 1.003951733380953559e+00;
	break;
	case 48000:
		p0 = -5.804669671441499519e-18;
		p1 = 3.053789928757690291e-13;
		p2 = -6.496014866101249794e-09;
		p3 = -3.969277247952757479e-05;
		p4 = 1.003952591709977415e+00;
	break;
	case 50000:
		p0 = -4.930039877651698485e-18;
		p1 = 2.701762164752109105e-13;
		p2 = -5.986696651033372113e-09;
		p3 = -3.810515833478767255e-05;
		p4 = 1.003954755096216900e+00;
	break;
	case 50400:
		p0 = -4.775358854140600092e-18;
		p1 = 2.637937027466497354e-13;
		p2 = -5.892041055138477995e-09;
		p3 = -3.780275477003837485e-05;
		p4 = 1.003955167234265478e+00;
	break;
	case 88200:
		p0 = -5.090139768480031043e-19;
		p1 = 4.921424693356908116e-14;
		p2 = -1.923829523675108346e-09;
		p3 = -2.160214687014651151e-05;
		p4 = 1.003977276279834019e+00;
	break;
	case 96000:
		p0 = -3.626648836744412987e-19;
		p1 = 3.816597595693573777e-14;
		p2 = -1.623898000952206408e-09;
		p3 = -1.984703070895103917e-05;
		p4 = 1.003979674995812754e+00;
	break;
	case 176400:
		p0 = -3.180720104052761169e-20;
		p1 = 6.151207927460715513e-15;
		p2 = -4.809399117661735116e-10;
		p3 = -1.080127073758081756e-05;
		p4 = 1.003992048707243789e+00;
	break;
	case 192000:
		p0 = -2.266250752087159745e-20;
		p1 = 4.770337959813220317e-15;
		p2 = -4.059609201723455372e-10;
		p3 = -9.923682425709201677e-06;
		p4 = 1.003993250132436454e+00;
	break;
	case 352800:
		p0 = -1.987755090799793513e-21;
		p1 = 7.688647436953224411e-16;
		p2 = -1.202327621003883834e-10;
		p3 = -5.400685657755166908e-06;
		p4 = 1.003999444875642189e+00;
	break;
	case 2822400:
		p0 = -4.852498403591686415e-25;
		p1 = 1.501626536939688520e-18;
		p2 = -1.878606337281021417e-12;
		p3 = -6.750912735944337505e-07;
		p4 = 1.004005921979573834e+00;
	break;
	case 5644800:
		p0 = -3.032792721276874507e-26;
		p1 = 1.877027577372803603e-19;
		p2 = -4.696510359527020341e-13;
		p3 = -3.375458367765195566e-07;
		p4 = 1.004006384824945863e+00;
	break;
	}

	return p0*pow(_SC,4) +
			p1*pow(_SC,3)+
			p2*pow(_SC,2)+
			p3*pow(_SC,1)+
			p4;
}

int Parameter::get_high_and_compute() {
}

double Parameter::errorfunc(double _G1, double _GC, const matrix<double, 4, 1> &Target) {
	double kT60 = Target(0),
			kED = Target(1),
			kC  = Target(2),
			kTC = Target(3);

	double G = G_G1_C(_G1, _GC, kC);
	double D1 = D1_G1_T60_C(_G1, _GC, kT60, kC);
	double DA = DA_D1_ED(D1, kED, 0.1);
	double T60est = compute_T60(D1, _G1, _GC, G);
	double EDest = compute_ED(D1, DA);
	double Cest = compute_C(_G1, _GC, G);
	double TCest = compute_TC(D1, _G1, DA);

	matrix <double, 4, 1> est_vec, orig_vec, err_vec;

	float w0 = 1.0, w1=1.0, w2=1.0, w3=1.0; // Weights for errors.

	est_vec(0) = w0*_N(T60est, T60_MIN, T60_MAX);
	est_vec(1) = w1*_N(EDest, ED_MIN, ED_MAX);
	est_vec(2) = w2*_N(Cest, C_MIN, C_MAX);
	est_vec(3) = w3*_N(TCest, TC_MIN, TC_MAX);

	orig_vec(0) = w0*_N(Target(0), T60_MIN, T60_MAX);
	orig_vec(1) = w1*_N(Target(1), ED_MIN, ED_MAX);
	orig_vec(2) = w2*_N(Target(2), C_MIN, C_MAX);
	orig_vec(3) = w3*_N(Target(3), TC_MIN, TC_MAX);

	err_vec = dlib::abs(est_vec-orig_vec);

	return mean(pow(err_vec,2)) + variance(pow(err_vec,2));

}

matrix<double, 5, 1> Parameter::inverse_mapping(double _T60, double _ED,
		double _C, double _TC, double _SC, double SampleRate, double _G1) {
    /*
     * _G1 is a starting point for the search of G1
     */

	matrix<double, 5, 1> result;
	for (int i=0; i<5; i++) result(i) = 0.0;

	double GC = GC_SC(_SC, SampleRate);
	result(3) = GC;

	matrix<double,1,1> starting_point;
	matrix<double,4,1> Target;

	Target(0) = _T60;
	Target(1) = _ED;
	Target(2) = _C;
	Target(3) = _TC;

	starting_point = _G1; // Starting point is current G1, for faster convergence

	// Lambda function for objective function to minimize.

	auto func = [&] (const matrix<double, 0, 1> &x) { float g1 = x(0); return errorfunc(g1, GC, Target); };

    find_min_box_constrained(bfgs_search_strategy(),
                             objective_delta_stop_strategy(1e-9),
							 func, derivative(func), starting_point, G1_MIN, G1_MAX);

    double G1 = starting_point(0);

    result(2) = G1;

    double G = G_G1_C(G1,GC,_C);


    result(4) = G;

    double D1 = D1_G1_T60_C(G1, GC, _T60, _C);


    result(0) = D1;

    double DA = DA_D1_ED(D1, _ED);
    result(1) = DA;

	return result;
}