// Josiah Boning
//
// frac.c
// Code to calculate the fractal dimension using various methods

#include "project.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/**
 * Calculate the fractal dimension of an audio file using the Variation method
 *
 * See Bigarelle and Iost, 1999
 * Currently rather inaccurate... need to implement cubic splines!
 */
double variation(double* datas, int size) {
	double d, tau;
	tau = 10.0/44100;
	d = 1.0/44100;
	//printf("d: %f\n", d);
	double sum = 0.0;
	double min,max;
	int x, i;
	for (x = 0; x < size; x++) {
		min = inf;
		max = -inf;
		for (i = x-1; i <= x+1; i++) {
			if (datas[i] < min) {
				min = datas[i];
			}
			if (datas[i] > max) {
				max = datas[i];
			}
		}
		sum += (max - min) * d;
	}
	//printf("sum: %f\n", sum);
	//printf("avg: %f\n", sum/size);
	//printf("log: %f\n", log(sum/size));
	//printf("den: %f\n", log(tau));
	//printf("frc: %f\n", log(sum/size) / log(tau));
	return 2.0 - (log(sum/size) / log(tau));
}

//double anam(double* datas, int size) {
//	double tau, dx, sum=0, subsum=0, alpha=1;
//	int i, j, l, k=1;
//	dx = 1.0/44100;
//	tau = k * dx;
//	for (i=k; i <= size - k; i++) {
//		subsum = 0;
//		for (j=0; j <= k; j++) {
//			for (l=0; l <= k; l++) {
//				subsum += pow(fabs(datas[i+j]-datas[i-l]), alpha);
//			}
//		}
//		sum += pow(subsum, 1.0/alpha);
//	}
//	double avg = pow(k+1,-2.0/alpha)/((double)size-2*k) * sum;
//	printf("avg: %f\n", avg);
//	printf("log: %f\n", log(avg));
//	printf("den: %f\n", log(tau));
//	printf("frc: %f\n", log(avg) / log(tau));
//	return 2.0 - (log(avg) / log(tau));
//}

/**
 * Calculate the fractal dimension of an audio file using the ANAM method
 *
 * See Bigarelle and Iost, 1999
 * Currently rather inaccurate... need to implement cubic splines!
 */
double anam(double* datas, int size) {
	double alpha, a, b, sum, tau=10;
	int x, t1, t2;
	a = tau; b = size-tau;
	sum = 0.0;
	for (x=a; x < b; x++) {
		double subsum = 0.0;
		for (t1 = 0; t1 <= tau; t1++) {
			for (t2 = 0; t2 <= tau; t2++) {
				subsum += fabs(datas[x+t1] - datas[x-t2]);
			}
		}
		//sum += fabs(datas[x] - datas[x-1]) + fabs(datas[x+1] - datas[x]) + fabs(datas[x+1] - datas[x-1]);
		sum += 1.0 / pow(tau,2) * subsum;
		//printf("%f, %f\n", subsum / pow(tau,2), sum / (x - a));
	}
	double avg = 1.0 / (b - a) * sum;
	//printf("sum: %f\n", sum);
	//printf("avg: %f\n", avg);
	//printf("log: %f\n", log(avg));
	//printf("den: %f\n", log(tau/44100.0));
	//printf("frc: %f\n", log(avg)/log(tau/44100.0));
	return 2.0 - (log(avg) / log(tau/44100.0));
}