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

added c++ class Parameter

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp/chourdakisreiss2016cpp/Parameter.cpp	Wed Dec 28 17:22:49 2016 +0000
@@ -0,0 +1,630 @@
+/*
+ * 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;
+}