基于 SP++ Matlab A 计权的 C++ 实现

/*
	int :data 输入数据
		 responseType true ：slow 1.0s    false：fast 0.125s

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

	if (responseType == true)
		duration = 1.0f;
	else
		duration = 0.125;

	int N1 = ceil(duration*Fs);
	N1 = NextPow2(N1);
	//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(ceil(datalength - N1,N1)+1);
	Vector<float>temp(N1);
	//float* temp = new float[N];
	for (int i = 0; i < dBA.size(); i++)
	{
		for (int j = 0; j < N1; j++) {
			temp[j] = data[i*N1 + 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;
}

double estimateLevel(Vector<float>data_in, float Fs, int A)
{
	Vector<double> X;
	Vector<double>data(data_in.size());
	for (int i = 0; i < data_in.size(); i++)
	{
		data[i] = (double)data_in[i];
	}
	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<double> f = linspace((double)0, (double)(X.size() - 1), X.size());
	f.operator*=(Fs);
	f.operator/=(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.
	double totalEnergy = sum(pow(X,(double)2)) / X.size();
	double meanEnergy = totalEnergy / ((1 / Fs)*data.size());
	double dBA = 10 * log10(meanEnergy) + A;
	return dBA;
}


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

	for (int i = 0; i < data.size(); i++)
	{
		if (data[i] == 0.f)
			data[i] = 1e-17;
	}

	data = splab::pow(data,(double)2);
	Vector<double> num = splab::pow(data, (double)4);
	num.operator*=(c1);
	Vector<double>den1 = splab::pow(data.operator+=(c2), (double)2);
	Vector<double>den2 = data.operator+=(c3);
	Vector<double>den3 = data.operator+=(c4);
	Vector<double>den4 = splab::pow(data.operator+=(c5), (double)2);
	Vector<double>den = den1*den2*den3*den4;
	Vector<double> A = num /den;
	//for (int i = 0; i < A.size(); i++)
	//{
	//	cout << A[i] << endl;
	//}
	return A;
}