channelestimation.cpp

数字音频广播中的信道估计算法在计算机上的算法实现

}

void CChannelEstimation::UpdateWienerFiltCoef(CReal rNewSNR, CReal rRatPDSLen,
											  CReal rRatPDSOffs)
{
	int	j, i;
	int	iDiff;
	int	iCurPil;

	/* Calculate all possible wiener filters */
	for (j = 0; j < iNumWienerFilt; j++)
	{
		matcWienerFilter[j] = FreqOptimalFilter(iScatPilFreqInt, j, rNewSNR,
			rRatPDSLen, rRatPDSOffs, iLengthWiener);
	}


#if 0
/* Save filter coefficients */
static FILE* pFile = fopen("test/wienerfreq.dat", "w");
for (j = 0; j < iNumWienerFilt; j++)
{
	for (i = 0; i < iLengthWiener; i++)
	{
		fprintf(pFile, "%e %e\n", Real(matcWienerFilter[j][i]),
			Imag(matcWienerFilter[j][i]));
	}
}
fflush(pFile);
#endif


	/* Set matrix with filter taps, one filter for each carrier */
	for (j = 0; j < iNumCarrier; j++)
	{
		/* We define the current pilot position as the last pilot which the
		   index "j" has passed */
		iCurPil = (int) (j / iScatPilFreqInt);

		/* Consider special cases at the edges of the DRM spectrum */
		if (iCurPil < iPilOffset)
		{
			/* Special case: left edge */
			veciPilOffTab[j] = 0;
		}
		else if (iCurPil - iPilOffset > iNumIntpFreqPil - iLengthWiener)
		{
			/* Special case: right edge */
			veciPilOffTab[j] = iNumIntpFreqPil - iLengthWiener;
		}
		else
		{
			/* In the middle */
			veciPilOffTab[j] = iCurPil - iPilOffset;
		}

		/* Special case for robustness mode D, since the DC carrier is not used
		   as a pilot and therefore we use the same method for the edges of the
		   spectrum also in the middle of robustness mode D */
		if (iDCPos != 0)
		{
			if ((iDCPos - iCurPil < iLengthWiener) && (iDCPos - iCurPil > 0))
			{
				/* Left side of DC carrier */
				veciPilOffTab[j] = iDCPos - iLengthWiener;
			}

			if ((iCurPil - iDCPos < iLengthWiener) && (iCurPil - iDCPos > 0))
			{
				/* Right side of DC carrier */
				veciPilOffTab[j] = iDCPos + 1;
			}
		}

		/* Difference between the position of the first pilot (for filtering)
		   and the position of the observed carrier */
		iDiff = j - veciPilOffTab[j] * iScatPilFreqInt;

		/* Copy correct filter in matrix */
		for (i = 0; i < iLengthWiener; i++)
			matcFiltFreq[j][i] = matcWienerFilter[iDiff][i];
	}
}

_REAL CChannelEstimation::CalAndBoundSNR(const _REAL rSignalEst,
										 const _REAL rNoiseEst)
{
	_REAL rReturn;

	/* "rNoiseEst" must not be zero */
	if (rNoiseEst != (_REAL) 0.0)
		rReturn = rSignalEst / rNoiseEst;
	else
		rReturn = (_REAL) 1.0;

	/* Bound the SNR at 0 dB */
	if (rReturn < (_REAL) 1.0)
		rReturn = (_REAL) 1.0;

	return rReturn;
}

_BOOLEAN CChannelEstimation::GetSNREstdB(_REAL& rSNREstRes) const
{
	const _REAL rNomBWSNR = rSNREstimate * rSNRSysToNomBWCorrFact;

	/* Bound the SNR at 0 dB */
	if ((rNomBWSNR > (_REAL) 1.0) && (bSNRInitPhase == FALSE))
		rSNREstRes = (_REAL) 10.0 * log10(rNomBWSNR);
	else
		rSNREstRes = (_REAL) 0.0;

	if (bSNRInitPhase == TRUE)
		return FALSE;
	else
		return TRUE;
}

_REAL CChannelEstimation::GetMSCMEREstdB()
{
	/* Calculate final result (signal to noise ratio) and consider correction
	   factor for average signal energy */
	const _REAL rCurMSCMEREst = rSNRSysToNomBWCorrFact *
		CalAndBoundSNR(AV_DATA_CELLS_POWER, rNoiseEstMSCMER);

	/* Bound the MCS MER at 0 dB */
	if (rCurMSCMEREst > (_REAL) 1.0)
		return (_REAL) 10.0 * log10(rCurMSCMEREst);
	else
		return (_REAL) 0.0;
}

_REAL CChannelEstimation::GetMSCWMEREstdB()
{
	_REAL rAvNoiseEstMSC = (_REAL) 0.0;
	_REAL rAvSigEstMSC = (_REAL) 0.0;

	/* Lock resources */
	Lock(); /* Lock start */
	{
		/* Average results from all carriers */
		for (int i = 0; i < iNumCarrier; i++)
		{
			rAvNoiseEstMSC += vecrNoiseEstMSC[i];
			rAvSigEstMSC += vecrSigEstMSC[i];
		}
	}
	Unlock(); /* Lock end */

	/* Calculate final result (signal to noise ratio) and consider correction
	   factor for average signal energy */
	const _REAL rCurMSCWMEREst = rSNRSysToNomBWCorrFact *
		CalAndBoundSNR(rAvSigEstMSC, rAvNoiseEstMSC);

	/* Bound the MCS MER at 0 dB */
	if (rCurMSCWMEREst > (_REAL) 1.0)
		return (_REAL) 10.0 * log10(rCurMSCWMEREst);
	else
		return (_REAL) 0.0;
}

_BOOLEAN CChannelEstimation::GetSigma(_REAL& rSigma)
{
	/* Doppler estimation is only implemented in the Wiener time interpolation
	   module */
	if (TypeIntTime == TWIENER)
	{
		rSigma = TimeWiener.GetSigma();
		return TRUE;
	}
	else
	{
		rSigma = (_REAL) 0.0;
		return FALSE;
	}
}

_REAL CChannelEstimation::GetDelay() const
{
	/* Delay in ms */
	return rLenPDSEst * iFFTSizeN / 
		(SOUNDCRD_SAMPLE_RATE * iNumIntpFreqPil * iScatPilFreqInt) * 1000;
}

_REAL CChannelEstimation::GetMinDelay()
{
	/* Lock because of vector "vecrDelayHist" access */
	Lock();

	/* Return minimum delay in history */
	_REAL rMinDelay = vecrDelayHist[0];
	for (int i = 0; i < iLenDelayHist; i++)
	{
		if (rMinDelay > vecrDelayHist[i])
			rMinDelay = vecrDelayHist[i];
	}

	Unlock();

	/* Return delay in ms */
	return rMinDelay * iFFTSizeN / 
		(SOUNDCRD_SAMPLE_RATE * iNumIntpFreqPil * iScatPilFreqInt) * 1000;
}

void CChannelEstimation::GetTransferFunction(CVector<_REAL>& vecrData,
											 CVector<_REAL>& vecrGrpDly,
											 CVector<_REAL>& vecrScale)
{
	/* Init output vectors */
	vecrData.Init(iNumCarrier, (_REAL) 0.0);
	vecrGrpDly.Init(iNumCarrier, (_REAL) 0.0);
	vecrScale.Init(iNumCarrier, (_REAL) 0.0);

	/* Do copying of data only if vector is of non-zero length which means that
	   the module was already initialized */
	if (iNumCarrier != 0)
	{
		_REAL rDiffPhase, rOldPhase;

		/* Lock resources */
		Lock();

		/* Init constants for normalization */
		const _REAL rTu = (CReal) iFFTSizeN / SOUNDCRD_SAMPLE_RATE;
		const _REAL rNormData = (_REAL) _MAXSHORT * _MAXSHORT;

		/* Copy data in output vector and set scale 
		   (carrier index as x-scale) */
		for (int i = 0; i < iNumCarrier; i++)
		{
			/* Transfer function */
			const _REAL rNormSqMagChanEst = SqMag(veccChanEst[i]) / rNormData;

			if (rNormSqMagChanEst > 0)
				vecrData[i] = (_REAL) 10.0 * Log10(rNormSqMagChanEst);
			else
				vecrData[i] = RET_VAL_LOG_0;

			/* Group delay */
			if (i == 0)
			{
				/* At position 0 we cannot calculate a derivation -> use
				   the same value as position 0 */
				rDiffPhase = Angle(veccChanEst[1]) - Angle(veccChanEst[0]);
			}
			else
				rDiffPhase = Angle(veccChanEst[i]) - rOldPhase;

			/* Take care of wrap around of angle() function */
			if (rDiffPhase > WRAP_AROUND_BOUND_GRP_DLY)
				rDiffPhase -= 2.0 * crPi;
			if (rDiffPhase < -WRAP_AROUND_BOUND_GRP_DLY)
				rDiffPhase += 2.0 * crPi;

			/* Apply normalization */
			vecrGrpDly[i] = rDiffPhase * rTu * 1000.0 /* ms */;

			/* Store old phase */
			rOldPhase = Angle(veccChanEst[i]);

			/* Scale (carrier index) */
			vecrScale[i] = i;
		}

		/* Release resources */
		Unlock();
	}
}

void CChannelEstimation::GetSNRProfile(CVector<_REAL>& vecrData,
									   CVector<_REAL>& vecrScale)
{
	int i;

	/* Init output vectors */
	vecrData.Init(iNumCarrier, (_REAL) 0.0);
	vecrScale.Init(iNumCarrier, (_REAL) 0.0);

	/* Do copying of data only if vector is of non-zero length which means that
	   the module was already initialized */
	if (iNumCarrier != 0)
	{
		/* Lock resources */
		Lock();

		/* We want to suppress the carriers on which no SNR measurement can be
		   performed (DC carrier, frequency pilots carriers) */
		int iNewNumCar = 0;
		for (i = 0; i < iNumCarrier; i++)
		{
			if (vecrSigEstMSC[i] != (_REAL) 0.0)
				iNewNumCar++;
		}

		/* Init output vectors for new size */
		vecrData.Init(iNewNumCar, (_REAL) 0.0);
		vecrScale.Init(iNewNumCar, (_REAL) 0.0);

		/* Copy data in output vector and set scale 
		   (carrier index as x-scale) */
		int iCurOutIndx = 0;
		for (i = 0; i < iNumCarrier; i++)
		{
			/* Suppress carriers where no SNR measurement is possible */
			if (vecrSigEstMSC[i] != (_REAL) 0.0)
			{
				/* Calculate final result (signal to noise ratio). Consider
				   correction factor for average signal energy */
				const _REAL rNomBWSNR =
					CalAndBoundSNR(vecrSigEstMSC[i], vecrNoiseEstMSC[i]) *
					rSNRFACSigCorrFact * rSNRSysToNomBWCorrFact;

				/* Bound the SNR at 0 dB */
				if ((rNomBWSNR > (_REAL) 1.0) && (bSNRInitPhase == FALSE))
					vecrData[iCurOutIndx] = (_REAL) 10.0 * Log10(rNomBWSNR);
				else
					vecrData[iCurOutIndx] = (_REAL) 0.0;

				/* Scale (carrier index) */
				vecrScale[iCurOutIndx] = i;

				iCurOutIndx++;
			}
		}

		/* Release resources */
		Unlock();
	}
}

void CChannelEstimation::GetAvPoDeSp(CVector<_REAL>& vecrData,
									 CVector<_REAL>& vecrScale,
									 _REAL& rLowerBound, _REAL& rHigherBound,
									 _REAL& rStartGuard, _REAL& rEndGuard,
									 _REAL& rPDSBegin, _REAL& rPDSEnd)
{
	/* Lock resources */
	Lock();

	TimeSyncTrack.GetAvPoDeSp(vecrData, vecrScale, rLowerBound,
		rHigherBound, rStartGuard, rEndGuard, rPDSBegin, rPDSEnd);

	/* Release resources */
	Unlock();
}