#pragma once
#include <stdlib.h>
#include "xyz.h"
#include <cmath>
#include <valarray>
#include <complex>

using namespace std;
class Fourier;
typedef complex<double> Complex;
typedef valarray<Complex> CArray;
const double PI = 3.141592653589793238460;

class window
{
private:
	const int size;
	xyz* elements;
public:
	CArray* fftx;
	CArray* ffty;
	CArray* fftz;
	int curNumOfElements;
	double meanx, meany, meanz;
	double correlationXY, correlationXZ, correlationYZ;
	window(int size) :size(size), curNumOfElements(0)
	{
		elements = (xyz*)malloc(size * sizeof(xyz));
		fftx = new CArray(size);
		ffty = new CArray(size);
		fftz = new CArray(size);
	}
	~window()
	{
		free(elements);
		delete(fftx);
		delete(ffty);
		delete(fftz);
	}
	bool add(xyz &element)
	{
		if (isFull()) return false;
		else
		{
			elements[curNumOfElements] = element;
			curNumOfElements++;
			return true;
		}
	}
	bool isFull()
	{
		return curNumOfElements >= size;
	}
	bool isOverthreshold()
	{
		return curNumOfElements >= size / 2;
	}
	void calculate()
	{
		calculateMean();
		calculateCorrelation();
		calculateFFT();
	}
	void calculateFFT() {
		for (int i = 0; i < curNumOfElements; i++) {
			fftx[i] = elements[i].x;
			ffty[i] = elements[i].y;
			fftz[i] = elements[i].z;
		}
		fft_axis(*fftx);
		fft_axis(*ffty);
		fft_axis(*fftz);
	}

	void fft_axis(CArray& a) {
		const size_t N = a.size();
		if (N <= 1) return;

		// divide
		CArray even = a[slice(0, N / 2, 2)];
		CArray  odd = a[slice(1, N / 2, 2)];

		// conquer
		fft_axis(even);
		fft_axis(odd);

		// combine
		for (size_t k = 0; k < N / 2; ++k)
		{
			Complex t = polar(1.0, -2 * PI * k / N) * odd[k];
			a[k] = even[k] + t;
			a[k + N / 2] = even[k] - t;
		}
	}

	void calculateMean()
	{
		int sumx = 0;
		int sumy = 0;
		int sumz = 0;
		for (int i = 0; i < curNumOfElements; i++)
		{
			sumx += elements[i].x;
			sumy += elements[i].y;
			sumz += elements[i].z;
		}
		meanx = sumx / curNumOfElements;
		meany = sumy / curNumOfElements;
		meanz = sumz / curNumOfElements;
	}
	void calculateCorrelation()
	{
		double tempXY = 0;
		double tempXZ = 0;
		double tempYZ = 0;
		double tempVarx = 0;
		double tempVary = 0;
		double tempVarz = 0;
		for (int i = 0; i < curNumOfElements; i++)
		{
			tempXY += (elements[i].x - meanx) * (elements[i].y - meany);
			tempXZ += (elements[i].x - meanx) * (elements[i].z - meanz);
			tempYZ += (elements[i].y - meany) * (elements[i].z - meanz);
			tempVarx += (elements[i].x - meanx) * (elements[i].x - meanx);
			tempVary += (elements[i].y - meany) * (elements[i].y - meany);
			tempVarz += (elements[i].z - meanz) * (elements[i].z - meanz);
		}
		double covXY = tempXY / curNumOfElements;
		double covXZ = tempXZ / curNumOfElements;
		double covYZ = tempYZ / curNumOfElements;
		double varX = tempVarx / curNumOfElements;
		double varY = tempVary / curNumOfElements;
		double varZ = tempVarz / curNumOfElements;
		correlationXY = covXY / sqrt((varX*varY));
		correlationXZ = covXZ / sqrt((varX*varZ));
		correlationYZ = covXZ / sqrt((varZ*varY));
		cout << "correlationXY is:" << correlationXY << endl;
	}

	xyz getElement(int i) const {
		return elements[i];
	}

};