/*
 * DFTProcessor.cpp
 *
 *  Created on: 29.04.2011
 *      Author: sven
 */
#ifndef DFTPROCESSOR_H_
	#include "DFTProcessor.h"
#endif

#ifndef _GLIBCXX_IOSTREAM
	#include <iostream>
#endif

#ifndef _GLIBCXX_FSTREAM
	#include <fstream>
#endif

#include <cmath>
#include <string.h>
#include <stdlib.h>
#include <iomanip>

const double PI = 3.14159265358979323846264338327950288419716939;

const char* ORIG_FILENAME = "orig.wav";
const char* IDFT_FILENAME = "idft.wav";

DFTProcessor::DFTProcessor(DFTAppParameters& params)
:mParams(params)
,mNumDftSamples(mParams.getNumSamples() / 2)
,mInputData(mParams.getNumSamples())
,mRe(mNumDftSamples)
,mIm(mNumDftSamples)
,mIdft(mParams.getNumSamples())
{
	std::cout << "DFT processing " << mParams.getNumSamples() <<
			" samples of" << mParams.getInputFilename() << std::endl;
	std::ifstream inFile;
	std::ofstream origFile;
	std::ofstream idftFile;
	inFile.open(mParams.getInputFilename().c_str(), std::ios::in | std::ios::binary);
	origFile.open(ORIG_FILENAME, std::ios::out | std::ios::binary);
	idftFile.open(IDFT_FILENAME, std::ios::out | std::ios::binary);
	int dataStartOffset = 0x2C;
	char header[dataStartOffset];
	inFile.read(&header[0],dataStartOffset);
	memcpy(&mSampleRate,((&header[0])+0x18),4);
	origFile.write((char*) &header[0],dataStartOffset);
	idftFile.write((char*) &header[0],dataStartOffset);
	short data;
	std::cout << "Original file:" << std::endl;
	for (unsigned i=0;i<mParams.getNumSamples();i++)
	{
		inFile.read((char*) &data,sizeof(short));
		mInputData.at(i) = (double) data;
		origFile.write((char*) &data, sizeof(short));
		if (i < 5) {
			if (i != 0) {
				std::cout << std::endl;
			}
			std::cout << mInputData.at(i);
		}
	}
	inFile.close();
	origFile.close();

	double nd,id,Nd,xn;
	Nd = (double) mParams.getNumSamples();
	for (unsigned i = 0; i < mNumDftSamples; i++) {
		id = (double) i;
		mRe.at(i) = 0.0;
		mIm.at(i) = 0.0;
		for (unsigned n=0; n < mParams.getNumSamples(); n++) {
			nd = (double) n;
			xn = mInputData.at(n);
			mRe.at(i) += xn * std::cos(2.0 * PI * nd * id / Nd );
			mIm.at(i) += xn * std::sin(2.0 * PI * nd * id / Nd );
		}
		mIm.at(i) = -mIm.at(i);
	}
	std::cout << std::endl << "dft done" << std::endl;

	// idft
	double rex,imx,fd;
	std::cout << "idft file:" << std::endl;
	for (unsigned n=0 ; n<mParams.getNumSamples() ; n++)
	{
		nd = (double) n;
		double reSum = 0.0;
		double imSum = 0.0;
		for (unsigned f=0; f < mNumDftSamples; f++)
		{
			fd = (double) f;
			if (f == 0) {
				rex = mRe.at(0) / Nd;
			} else if (f == mNumDftSamples) {
				rex = mRe.at(mNumDftSamples) / Nd;
			} else {
				rex = mRe.at(f) / (Nd / 2.0);
			}
			imx = - mIm.at(f) / (Nd / 2.0);
			reSum += rex * std::cos(2.0 * PI * nd * fd / Nd);
			imSum += imx * std::sin(2.0 * PI * nd * fd / Nd);
		}
		mIdft.at(n) = reSum + imSum;
	}
	for (unsigned n=0 ; n<mParams.getNumSamples(); n++) {
		data = (short) round(mIdft.at(n));
		idftFile.write((char*) &data, sizeof(short));
		if (n < 5) {
			if (n != 0) {
				std::cout << std::endl;
			}
			double dout = mIdft.at(n);
			std::cout << std::setprecision(16) <<  dout;
		}
	}
	idftFile.close();
	std::cout << std::endl << "idft done" << std::endl;

	// amplitude spectrum
	double A1, A2;
	unsigned maxTimeDomFreq = mSampleRate / 2;
	unsigned freqIntervalWidth = maxTimeDomFreq / mNumDftSamples;
	for (unsigned f = 0 ; f < mNumDftSamples ; f++) {
		A1 = mIm.at(f);
		A2 = mRe.at(f);
		SpectralData actSpecData;
		actSpecData.mAmplitude = std::sqrt( (A1*A1) + (A2*A2) );
		actSpecData.mPhase = std::atan(A2/A1);
		actSpecData.mFreq = f * freqIntervalWidth;
		mSpectrum.push_back(actSpecData);
	}
	writeSpectrumToCSV();
	startGnuPlot();
}

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

void DFTProcessor::writeSpectrumToCSV()
{
	std::ofstream csvFile;
	csvFile.open("spectrum.csv", std::ios::out);
	std::vector<SpectralData>::const_iterator specIter;
	for(specIter = mSpectrum.begin() ; specIter != mSpectrum.end() ; specIter++) {
		const SpectralData& actSpec = *specIter;
		if (specIter != mSpectrum.begin()) {
			csvFile << std::endl;
		}
		csvFile << actSpec.mFreq << ";";
		csvFile << std::setprecision (12) << std::fixed << actSpec.mAmplitude;
		csvFile << ";";
		csvFile << std::setprecision (12) << std::fixed << actSpec.mPhase;
	}
	csvFile.close();
}

void DFTProcessor::startGnuPlot()
{
	int i = system("gnuplot spectrum.plt -p");
	if (i != 0) {
		std::cerr << "Error starting gnuplot: " << i << std::endl;
	}
}