
#include "my_signal.h"
#include "my_point.h"
#include <iostream>
#include <limits>


Signal::Signal(){
	sourceFile="";
	maxY=0;
	maxX=0;
	minY=0;
	minX=0;
	mean=0;
	variance=0;
	dispersion=0;
	stdDeviation=0;
	avgDeviation=0;
	skewness=0;
	kurtosis=0;
	minXdist=1;
	minYdist=1;
	points.clear();
}

Signal::Signal(const char * fname){
	sourceFile=fname;
	ifstream in(fname);
	if (!in.is_open()) return;
	while (!in.eof()){
		Point tmp;
		in>>tmp;
		points.push_back(tmp);
	}
	calcStats();
}

Signal::Signal(const Signal & sig){
	sourceFile=sig.sourceFile;
	points=sig.points;
	calcStats();
}

Signal::Signal(vector<Point> pts){
	points=pts;
	sourceFile="";
	calcStats();
}

Signal::~Signal(){
}

// Operators
Point & Signal::operator[](ul index) {
	return points[index];
}

ofstream & operator<<(ofstream & out, const Signal &sig){
	for (ul i=0; i<sig.points.size()-1; i++){
		out<<sig.points[i].x<<" "<<sig.points[i].y<<endl;
	}
	out<<sig.points[sig.points.size()-1].x<<" "<<sig.points[sig.points.size()-1].y;

	return out;
}

ostream & operator<<(ostream & out, const Signal &sig){
	for (ul i=0; i<sig.points.size()-1; i++){
		out<<sig.points[i].x<<" "<<sig.points[i].y<<endl;
	}
	out<<sig.points[sig.points.size()-1].x<<" "<<sig.points[sig.points.size()-1].y;

	return out;
}

ifstream & operator>>(ifstream & in, Signal &sig){
	if (!in.is_open()) return in;
	while (!in.eof()){
		Point tmp;
		in>>tmp;
		sig.points.push_back(tmp);
	}
	sig.calcStats();
	return in;
}

Signal& Signal::operator+=(const Point& pt){
	points.push_back(pt);
	calcStats();
	return *this;
}

/*
 * To have meaning adding 2 Signals we need to ensure that:
 *	- They both have the same step on x-values
 *	- The max of minX1 n minX2 exists in the signal
 *	  with the min of minX1 n minX2
 */
Signal Signal::operator+(const Signal& sig){
	// Empty sig to return
	Signal ret;

	// Check the signals if big enought...to avoid segmentation fault
	if (points.size()<2 || sig.points.size()<2) return ret;

	// Check the step     .... WRONG DOUBLE COMPARSION??? 0.2 != 0.2 ?????????
	/*
	 * The 2 steps may differ for 3.123332e-14... even if there are "equal"
	 * this happens because of the precision used on the memory for floats
	 * and doubles. We will accept only the first MAX_DEC_DIGITS of the number
	 * while comparing and "cut" the others by assigning the value
	 * pow(10.0,MAX_DEC_DIGITS)*stepdiff to an int.
	 */
	double stepdiff=(points[1].x-points[0].x)-(sig.points[1].x-sig.points[0].x);
	long long int accepted_diff=(long long int)(pow(10.0,MAX_DEC_DIGITS)*stepdiff);
	
	if (accepted_diff!=0) return ret;


	// Find the signal with the min minX...
	Signal min_minXsig;
	Signal max_minXsig;
	if (minX < sig.minX){
		min_minXsig=*this;
		max_minXsig=sig;
	}
	else{
		min_minXsig=sig;
		max_minXsig=*this;
	}

	// Scan min_minXsig for minX of max_minXsig...
	ul max_start=0;
	bool found=false;
	
	for (ul i=0; i<min_minXsig.size(); i++){
		// I HATE DOUBLE COMPS
	
		if ((long long int)(min_minXsig[i].x*pow(10.0,MAX_DEC_DIGITS))== (long long int)(max_minXsig.getMinX()*pow(10.0,MAX_DEC_DIGITS))) {
			found=true;
			max_start=i;
			break;
		}
	}
	
	// Check if ok
	if (!found) return ret;

	// Populate points of the new signal
	
	ul minSigPos=0;
	ul maxSigPos=0;


	for (ul i=0; ; i++){
	
		if (min_minXsig.size()<=minSigPos && max_minXsig.size()<=maxSigPos){
			break;
		}else if (min_minXsig.size()<=minSigPos) {
			// min_minXsig has ended
			ret+=max_minXsig[maxSigPos];
			maxSigPos++;
		}
		else if (max_minXsig.size()<=maxSigPos){
			// max_minXsig has ended
			ret+=min_minXsig[minSigPos];
			minSigPos++;
		}else if (i<max_start){
			// Min sig only in this area
			ret+=min_minXsig[minSigPos];
			minSigPos++;
		}
		else {
			// Real adding is done here(common x areas)...
			Point p;
			p.x=min_minXsig[minSigPos].x;
			p.y=min_minXsig[minSigPos].y+max_minXsig[maxSigPos].y;
			// Fix ZERO PROBLEM... (0.00000000000000001263 or -1.9436e -17)
			if ((long long int)(p.y*pow(10.0,MAX_DEC_DIGITS))==0) p.y=0.0;
			maxSigPos++;
			minSigPos++;
			ret+=p;
		}
	}
	
	return ret;
}

Signal Signal::operator*(const Signal& sig){

	// Empty sig to return
	Signal ret;

	// Check the signals if big enought...to avoid segmentation fault
	if (points.size()<2 || sig.points.size()<2) return ret;

	double stepdiff=(points[1].x-points[0].x)-(sig.points[1].x-sig.points[0].x);
	long long int accepted_diff=(long long int)(pow(10.0,MAX_DEC_DIGITS)*stepdiff);
	
	if (accepted_diff!=0) return ret;
	
	// Find the signal with the min minX...
	Signal min_minXsig;
	Signal max_minXsig;
	if (minX < sig.minX){
		min_minXsig=*this;
		max_minXsig=sig;
	}
	else{
		min_minXsig=sig;
		max_minXsig=*this;
	}
	
	// Scan min_minXsig for minX of max_minXsig...
	bool found=false;
	
	for (ul i=0; i<min_minXsig.size(); i++){
		if ((long long int)(min_minXsig[i].x*pow(10.0,MAX_DEC_DIGITS))== (long long  int)(max_minXsig.getMinX()*pow(10.0,MAX_DEC_DIGITS))) {
			found=true;
			break;
		}
	}
	
	// Check if ok
	if (!found) return ret;

	// Populate points of the new signal
	
	ul minSigPos=0;
	ul maxSigPos=0;


	for (ul i=0; ; i++){
		// Check FIRST if something ENDED
		if (min_minXsig.size()<=minSigPos && max_minXsig.size()<=maxSigPos){
			break;
		}else if (min_minXsig.size()<=minSigPos) {
			// min_minXsig has ended
			ret+=Point(max_minXsig[maxSigPos].x,0);
			maxSigPos++;
		}else if ((long long int)(min_minXsig[i].x*pow(10.0,MAX_DEC_DIGITS)) < (long long  int)(max_minXsig.getMinX()*pow(10.0,MAX_DEC_DIGITS)) || max_minXsig.size()<=maxSigPos){
			// Min sig only in this area
			ret+=Point(min_minXsig[minSigPos].x,0);
			minSigPos++;
		}
		else {
			// Real Multiplication is done here(common x areas)...
			Point p;
			p.x=min_minXsig[minSigPos].x;
			p.y=min_minXsig[minSigPos].y*max_minXsig[maxSigPos].y;
			maxSigPos++;
			minSigPos++;
			ret+=p;
		}
	}
	
	return ret;
	//Sorry For Copy/Paste :-))
}

/* Finite Impulse Response, FIR (Peperasmenh kroustikh apokrish)
 *
 * y(n) = S(k=-M/2+1,k=M/2) x(k)h(n-k) ... M=length, x(M)???
 * y(n) = S(k=0,k=N) x(k)h(n-k) ... N=x.size + h.size -1
 */
Signal Signal::convolution(const Signal& h){
	Signal ret;
	
	// Check the signals if big enought...to avoid segmentation fault
	if (points.size()<2 || h.points.size()<2) return ret;

	double stepdiff=(points[1].x-points[0].x)-(h.points[1].x-h.points[0].x);
	long long int accepted_diff = (long long int)(pow(10.0,MAX_DEC_DIGITS)*stepdiff);
	

	if (accepted_diff!=0) return ret;
	
	double step=points[1].x-points[0].x;

	ul N = points.size()+h.points.size()-1;
	
	for (ul n=0; n<N; n++){
		Point np(n*step,0);
		
		
		for (ul k=0; k<points.size(); k++){
			
			double h_y=0;
			double x_y=0;
			
			if ( n-k > 0 && n-k <h.points.size()) h_y=h.points[n-k].y; 
			x_y=points[k].y;
			 
			np.y += x_y*h_y;
		}
		ret+=np;
	}

	
	return ret;
}

void Signal::fft(Signal &real, Signal &imag, Signal &power){
	// Fix signal size
	ul size = sizeToPowOfTow(points.size());

	// Init 2 double *
	double real_data[size];
	for (ul i=0; i<size; i++) real_data[i]=points[i].y;
	
	double imag_data[size];
	for (ul i=0; i<size; i++) imag_data[i]=0.0;
	
	// Real Call to fft
	fft::fft(real_data, imag_data, size, 1);
		
	// Create Signals
	for (ul i=0; i<size; i++){
		real+=Point(points[i].x,real_data[i]);
		imag+=Point(points[i].x,imag_data[i]);
		power+=Point(points[i].x,sqrt( pow(real_data[i],2)+pow(imag_data[i],2 ) ) );
	}
	
}




// Public methods
string Signal::toString() const{
	string temp="";
	for (ul i=0; i<points.size(); i++){
		char line[50];
		sprintf((char *)line,"%f %f\n",points[i].x,points[i].y);
		temp+=line;
	}
	
	return temp;
}
void Signal::save(const char * fname){
	
	sourceFile=fname;
	ofstream out(fname);
	out<<*this;
	out.close();
}

void Signal::load(const char * fname){
	//Clear Vector , if any
	points.clear();

	ifstream in(fname,ifstream::in);
	if (!in.is_open()) return;

	sourceFile=fname;
	
	while (!in.eof()){
		Point tmp;
		in>>tmp;
		points.push_back(tmp);
	}
	
	calcStats();
}

bool Signal::isEmpty(){
	return points.size()==0;
}

void Signal::moveHor(double pts){
	for (ul i=0; i<points.size(); i++)
		points[i].x+=pts;
	calcStats();
}

void Signal::moveVer(double pts){
	for (ul i=0; i<points.size(); i++)
		points[i].y+=pts;
	calcStats();
}

ul Signal::size(){
	return points.size();
}

void Signal::addNoise(double max){
	srand(time(NULL));
	for (ul i=0; i<points.size(); i++){
		double r = (double)rand() / ((double)(RAND_MAX) +(double)1); // r -> [0,0.1)
		// decide a sign
		int up_down=1;
		if (r<0.5) up_down=-1;
		// make percent
		double per=r*max;  // per -> [0,max)
		points[i].y+=(per*(maxY-minY)/100) * up_down;
	}
	calcStats();
}

string Signal::getFileName(){
	return sourceFile;
}

Signal Signal::getSampledBy(int step){
	Signal ret;
	for (ul i=0; i<points.size(); i+=step)
		ret+=points[i];
	return ret;
}

Signal Signal::getQuantedWith(int levels){
	Signal ret;
	double qStep=(maxY-minY) / (levels-1);

	for (ul i=0; i<points.size(); i++){
		Point p;
		p.x=points[i].x;
		
		double dist1=fabs( fmod(points[i].y-minY,qStep) );	//<- going to "down"
		double dist2=fabs( fmod((maxY-points[i].y),qStep) );	//<- away from "up"
		
		if (dist2<dist1) p.y=points[i].y+dist2;
		else  p.y=points[i].y-dist1;
		
		ret+=p;
	}
	return ret;
}

double Signal::getMean(){return mean;}
double Signal::getMaxY(){return maxY;}
double Signal::getMaxX(){return maxX;}
double Signal::getMinY(){return minY;}
double Signal::getMinX(){return minX;}
double Signal::getVariance(){return variance;}
double Signal::getDispersion(){return dispersion;}
double Signal::getStdDeviation(){return stdDeviation;}
double Signal::getAvgDeviation(){return avgDeviation;}
double Signal::getSkewness(){return skewness;}
double Signal::getKurtosis(){return kurtosis;}
double Signal::getMinXdist(){return minXdist;}
double Signal::getMinYdist(){return minYdist;}

// Private
void Signal::calcStats(){
	// Update All Statistics
	calcMinMaxMean();
	calcVariance();
	calcDispersion();
	calcStdDeviation();
	calcAvgDeviation();
	calcSkewness();
	calcKurtosis();
}

void Signal::calcMinMaxMean(){
	
	// init
	minX=points[0].x;
	minY=points[0].y;
	maxX=points[0].x;
	maxY=points[0].y;
	
	minXdist=10000.0;
	minYdist=10000.0;
	
	double sum=points[0].y;
	
	//Scan Vector
	for(ul i=1; i<points.size(); i++){
		// Min
		if (points[i].x < minX) minX=points[i].x;
		if (points[i].y < minY) minY=points[i].y;
		// Max
		if (points[i].x > maxX) maxX=points[i].x;
		if (points[i].y > maxY) maxY=points[i].y;
		//Sum
		sum+=points[i].y;
		// Min distance
		double xdist = fabs(points[i].x - points[i-1].x);
		double ydist = fabs(points[i].y - points[i-1].y);
		if (xdist<minXdist && xdist>0) minXdist=xdist;
		if (ydist<minYdist && ydist>0) minYdist=ydist;
	}
	
	mean=sum/(double)points.size();
}

void Signal::calcVariance(){
	double sum=0;
	for (ul i=0; i<points.size(); i++){
		sum+=pow(points[i].y-mean,2);
	}
	variance=sum/(double)(points.size()-1);
}

void Signal::calcDispersion(){}

void Signal::calcStdDeviation(){
	stdDeviation=sqrt(variance);
}

void Signal::calcAvgDeviation(){
	double sum=0;
	for (ul i=0; i<points.size(); i++){
		sum+=fabs(points[i].y-mean);
	}
	avgDeviation=sum/(double)points.size();
}

void Signal::calcSkewness(){
	double sum=0;
	for (ul i=0; i<points.size(); i++){
		sum+=pow((points[i].y-mean)/stdDeviation,3);
	}
	skewness=sum/(double)points.size();
}

void Signal::calcKurtosis(){
	double sum=0;
	for (ul i=0; i<points.size(); i++){
		sum+=pow((points[i].y-mean)/stdDeviation,4);
	}
	kurtosis=(sum/(double)points.size())-3;
}


ul Signal::sizeToPowOfTow(ul size){
	ul old_size=0;
	
	for (int i=0; i<1000; i++){
		if (pow((double)2,i)>size) break;
		old_size=(ul) pow((double)2,i);
	}
	return old_size;
}
