/*
 * WaveAnalyzer.cpp
 *
 *  Created on: 02.04.2011
 *      Author: sven
 */
#ifndef WAVEANALYZER_H_
	#include "WaveAnalyzer.h"
#endif

#ifndef _GLIBCXX_FSTREAM
	#include <fstream>
#endif

#ifndef _GLIBCXX_IOSTREAM
	#include <iostream>
#endif

#ifndef _GLIBCXX_CMATH
	#include <cmath>
#endif

#ifndef _GLIBCXX_MAP
	#include <map>
#endif

#ifndef _GLIBCXX_NUMERIC_LIMITS
	#include <limits>
#endif

WaveAnalyzer::WaveAnalyzer() {
	// TODO Auto-generated constructor stub

}

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

long double WaveAnalyzer::calcSNR(std::string& fullQuantFilename, std::string& lowQuantFilename)
{
	long double result = 0.0;
	double sumXnSquare = 0.0;
	double sumDiffSquare = 0.0;
	std::ifstream fullQuantFile;
	std::ifstream lowQuantFile;
	fullQuantFile.open(fullQuantFilename.c_str(), std::ios::in|std::ios::binary);
	lowQuantFile.open(lowQuantFilename.c_str(), std::ios::in|std::ios::binary);
	int dataLenOffset = 0x28;
	fullQuantFile.seekg(dataLenOffset,std::ios_base::cur);
	int numSamples = 0;
	fullQuantFile.read((char*) &numSamples, sizeof(numSamples));
	lowQuantFile.seekg(dataLenOffset+sizeof(numSamples),std::ios_base::cur);
	double fullVal = 0.0;
	double lowVal = 0.0;
	short data = 0;
	int diff = 0;
	for(int i=0; i < numSamples; ) {
		fullQuantFile.read((char*) &data, sizeof(data));
		fullVal = (double)(data);
		//std::cout << "full: " << data << " "<< fullVal << std::endl;
		lowQuantFile.read((char*) &data, sizeof(data));
		lowVal = (double)(data);
		//std::cout << "low: " << data << " "<< lowVal << std::endl;
		sumXnSquare += (fullVal * fullVal);
		//std::cout << "Sum Xn Square: " << sumXnSquare << std::endl;
		diff = lowVal - fullVal;
		sumDiffSquare += (diff * diff);
		//std::cout << "Sum Diff Square: " << sumDiffSquare << std::endl;
		i += sizeof(data);
	}
	if (sumDiffSquare != 0.0) {
		result = 10.0 * std::log10(sumXnSquare / sumDiffSquare);
	} else {
		result = std::numeric_limits<long double>::infinity();
	}
	fullQuantFile.close();
	lowQuantFile.close();
	return result;
}

long double WaveAnalyzer::calcArithAvg(std::string& wavefilename)
{
	long double result = 0.0;
	int sum = 0;
	std::ifstream waveFile;
	waveFile.open(wavefilename.c_str(), std::ios::in|std::ios::binary);
	int dataLenOffset = 0x28;
	waveFile.seekg(dataLenOffset,std::ios_base::cur);
	int numSamples = 0;
	waveFile.read((char*) &numSamples, sizeof(numSamples));
	short data = 0;
	for(int i=0 ; i < numSamples; ) {
		waveFile.read((char*) &data, sizeof(data));
		sum += data;
		i+=sizeof(data);
	}
	result = ((long double) sum)/((long double)numSamples);
	waveFile.close();
	return result;
}

long double WaveAnalyzer::calcVariance(std::string& wavefilename)
{
	long double variance = 0.0;
	std::map<short,statVal*> statMap;
	std::ifstream waveFile;
	waveFile.open(wavefilename.c_str(), std::ios::in|std::ios::binary);
	UInt32 dataLenOffset = 0x28;
	waveFile.seekg(dataLenOffset,std::ios_base::cur);
	UInt32 numBytes = 0;
	waveFile.read((char*) &numBytes, sizeof(numBytes));
	UInt32 numSamples = numBytes / sizeof(short);
	for(UInt32 i=0 ; i < numSamples; i++) {
		short data = 0;
		waveFile.read((char*) &data, sizeof(data));
		std::map<short,statVal*>::iterator findIt = statMap.find(data);
		if (findIt == statMap.end()) {
			statVal* vPtr = new statVal;
			vPtr->count = 1UL;
			vPtr->p = 0.0;
			statMap.insert(std::make_pair(data,vPtr));
		} else {
			findIt->second->count++;
		}
	}
	waveFile.close();
	long double expectation = 0.0;
	std::map<short,statVal*>::iterator it;
	long double pSum = 0.0;
	for(it = statMap.begin() ; it != statMap.end(); it++)
	{
		long double x = (long double) it->first;
		it->second->p = (long double) it->second->count / (long double) numSamples;
		long double actExp = (x * it->second->p);
		expectation += actExp;
		pSum += it->second->p;
	}
	std::cout << "expectation: " << expectation  << " pSum: " << pSum << std::endl;
	for(it = statMap.begin() ; it != statMap.end(); it++)
	{
		long double x = (long double) it->first;
		long double diff = x - expectation;
		long double diffQuad = (diff*diff);
		variance += diffQuad * it->second->p;
	}
	for (it = statMap.begin() ; it != statMap.end();) {
		delete it->second;
		it++;
	}
	statMap.clear();
	return variance;
}

long double WaveAnalyzer::calcStdDeviation(std::string& wavefilename)
{
	return std::sqrt(calcVariance(wavefilename));
}