/*
	int :data 輸入資料
		 responseType 0 ：slow 1.0s    1：fast 0.125s

*/
Vector<float> A_weight(float * data,int datalength,int responseType,float Fs) {
	//% Note: A quite location will be ~55 dBA.
	int A = 55;
	float duration;

	if (responseType == 0)
		duration = 1.0f;
	else if (responseType == 1)
		duration = 0.125;
	else
		return 0.f;

	int N = ceil(duration*Fs);
	N = NextPow2(N);
	//Vector<int> windowStart(2 + (datalength - N) % N);
	//for (int i = 1; i < (datalength - N) % N; i++)
	//{
	//	windowStart[i] = 1 + i*N;
	//}
	//windowStart[0] = 1;
	//windowStart[windowStart.size() - 1] = datalength - N ;
	//Vector<float>dBA(windowStart.size(), 0.f);
	Vector<float>dBA(2 + (datalength - N) % N);
	Vector<float>temp(N);
	//float* temp = new float[N];
	for (int i = 0; i < dBA.size(); i++)
	{
		for (int j = 0; j < N; j++) {
			temp[j] = data[i*N + j];
		}
		dBA[i] = estimateLevel(temp, Fs, A);
	}

	return dBA;
}
//臨近的較大的2的整數次冪
unsigned int NextPow2(unsigned int v)
{
	v--;
	v |= v >> 1;
	v |= v >> 2;
	v |= v >> 4;
	v |= v >> 8;
	v |= v >> 16;
	v++;
	return v;
}

float estimateLevel(Vector<float>data, float Fs, int A)
{
	Vector<float> X;
	X= abs(fft(data));

	//% Add offset to prevent taking the log of zero.
	for (int i = 0; i < X.size(); i++)
	{
		if (X[i] == 0.f)
			X[i] = 1e-17;
	}

	//% Retain frequencies below Nyquist rate.
	Vector<float> f = linspace(0.f, (float)(X.size() - 1), X.size()).operator*=(Fs / X.size());
	X= splab::wkeep(X, X.size() / 2, 0);
	f = splab::wkeep(f, f.size() / 2, 0);
	//% Apply A - weighting filter.
	X = filterA(f)*X;
	//% Estimate dBA value using Parseval's relation.
	float totalEnergy = sum(sqrt(X)) / X.size();
	float meanEnergy = totalEnergy / ((1 / Fs)*data.size());
	float dBA = 10 * log10(meanEnergy) + A;
	return dBA;
}


Vector<float> filterA(Vector<float>data) {
	double c1 = 3.5041384e16;
	double c2 = sqrt(20.598997);
	double c3 = sqrt(107.65265);
	double c4 = sqrt(737.86223);
	double c5 = sqrt(12194.217);

	for (int i = 0; i < data.size(); i++)
	{
		if (data[i] == 0.f)
			data[i] = 1e-17;
	}
	data = sqrt(data);
	Vector<float> num = pow(data, 4.f);
	num.operator*=(c1);
	Vector<float>den = sqrt(data.operator+=(c2))* sqrt(data.operator+=(c3))*sqrt(data.operator+=(c4))*sqrt(data.operator+=(c5));
	Vector<float> A = num / den;
	return A;
}