FFT数据简单校正频率,相位,幅度

BOOL	CWaveBox::DataFftInfo(short *pData,int nCount, FFT_INFO &fft)
{
	Value2D		*pIn, *pOut;
	int			i, nPow, nIndex, nAllCount, nRealCount, nHi, nTotal;
	int			TopIndex[FFT_FI_CountMax*2];
	Value2D		v, v1, v2;
	double		ang1, ang2, ang, dIvt, dBase, dSumAm, dAm, sAmStd, cn, ang_dlt, f_dlt;
	double		*pFreqAm;
	int			x1, x2;
	double		y1, y2, a;
	BOOL		bLeft;

	pIn = (Value2D *)m_FftIn.GetPtr();
	pOut = (Value2D *)m_FftOut.GetPtr();
	pFreqAm = (double *)m_FftAm.GetPtr();

	nPow = (int)(log((double)nCount)/log(2.0));
	nPow = FORCE_RANGE(nPow, 10, 16);
	nAllCount = 1 << nPow;
	nRealCount = min(nCount, nAllCount);

	nHi = 0;
	nTotal = 0;
	memset(pIn, 0, sizeof(Value2D)*nAllCount);
	for(i=0; i<nRealCount; i++)
	{
		if(pData[i] > 0)
			nHi ++;
		pIn[i].Real = pData[i];
		if(pData[i])
			nTotal ++;
	}
	fft.dDutyR = 1.0*nHi / nTotal;
	NFourier::Fft_Double2(pIn, pOut, nPow);

	dIvt = 1.0/m_nFreq;
	dBase = NFourier::CalcuFreqAm(pOut, pFreqAm, nPow, dIvt);
	GetTopValuesBig(pFreqAm, nAllCount/2, TopIndex, FFT_FI_CountMax);

	for(i=0; i<FFT_FI_CountMax; i++)
	{
		nIndex = TopIndex[i];
		fft.Top[i].nIndex = nIndex;
		fft.Top[i].dFreq = dBase * nIndex;

		dAm = pFreqAm[nIndex];
		Level2Std((int)dAm, sAmStd);
		sAmStd = sAmStd * nAllCount / nRealCount;
		fft.Top[i].dAm = sAmStd;

		v1 = pOut[nIndex];
		ang1 = atan2(v1.Imgy, v1.Real);
		ang1 = ang1*360/2/PI;
		if(ang1<0) ang1+=360;
		fft.Top[i].dAng = ang1;
	}

	nIndex = TopIndex[0];
	if(nIndex == 0)
		nIndex = TopIndex[1];

	if(nIndex > 1)
	{
		x1 = nIndex;
		y1 = pFreqAm[x1];
		v1 = pOut[nIndex];

		if(pFreqAm[nIndex+1] >= pFreqAm[nIndex-1])
		{
			bLeft = FALSE;
			x2 = nIndex+1;
			v2 = pOut[x2];
			y2 = pFreqAm[x2];
			a = y2/(y1+y2);
			cn = x1+a;
		}
		else
		{
			bLeft = TRUE;
			x2 = nIndex-1;
			v2 = pOut[x2];
			y2 = pFreqAm[x2];
			a = y2/(y1+y2);
			cn = x1-a;
		}

		fft.dFreq = dBase * cn;
		f_dlt = dBase * a;
		dAm = y1 * a * PI / (sin(a*PI));
		fft.dAm = dAm;

		ang1 = atan2(v1.Imgy, v1.Real);
		ang2 = atan2(v2.Imgy, v2.Real);

		ang_dlt = f_dlt*PI*dIvt*(1<<nPow);
		ang_dlt = PI * a;
		ang_dlt = fabs(ang2-ang1) * a;
		if(bLeft)
			ang = PI/2 -ang1 + ang_dlt;
		else
			ang = PI/2 -ang1 - ang_dlt;

		ang = ang*360/2/PI;
		if(ang<0.0) ang+=360.0;
		fft.dAng = ang;

		NSys::Trace("ang1: %.1f ang2:%.1f ang_dlt:%.1f\n", ang1*180/PI, ang2*180/PI, ang_dlt*180/PI);
	}
	else
	{
		fft.dFreq = dBase * nIndex;
		fft.dAm = pFreqAm[nIndex];

		v1 = pOut[nIndex];
		ang1 = atan2(v1.Imgy, v1.Real);
		ang1 = ang1*360/2/PI;
		if(ang1<0.0) ang1+=360.0;
		fft.dAng = ang1;
	}

	dSumAm = 0.0;
	for(i=0; i<nAllCount/2; i++)
	{
		dSumAm += pFreqAm[i]*pFreqAm[i];
	}
	fft.dMainPst = fft.dAm*fft.dAm/dSumAm;
	NData::DbgArry10K(pFreqAm);

	Level2Std((int)fft.dAm, sAmStd);
	fft.dAm = sAmStd * nAllCount / nRealCount;

	CalcuDutyRatio(pData, nCount, fft);
	return TRUE;
}