channelestimation.cpp

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

					IIR1(rSignalEst, vecrSqMagChanEst[i],
						rLamSNREstFast);
				}
			}

			/* Calculate final result (signal to noise ratio) */
			rCurSNREst = CalAndBoundSNR(rSignalEst, rNoiseEst);

			/* Consider correction factor for average signal energy */
			rSNREstimate = rCurSNREst * rSNRFACSigCorrFact;
		}
		break;
	}


	/* WMER on MSC cells estimation ----------------------------------------- */
	for (i = 0; i < iNumCarrier; i++)
	{
		/* Use MSC cells for this SNR estimation method */
		if (_IsMSC(ReceiverParam.matiMapTab[iModSymNum][i]))
		{
			CReal rCurErrPow;

			/* Get tentative decision for this MSC QAM symbol and calculate
			   squared distance as a measure for the noise. MSC can be 16 or
			   64 QAM */
			switch (ReceiverParam.eMSCCodingScheme)
			{
			case CParameter::CS_3_SM:
			case CParameter::CS_3_HMSYM:
			case CParameter::CS_3_HMMIX:
				rCurErrPow = SqMag(MinDist64QAM((*pvecOutputData)[i].cSig));
				break;

			case CParameter::CS_2_SM:
				rCurErrPow = SqMag(MinDist16QAM((*pvecOutputData)[i].cSig));
				break;
			}

			/* Use decision together with channel estimate to get
			   estimates for signal and noise (for each carrier) */
			IIR1(vecrNoiseEstMSC[i], rCurErrPow * vecrSqMagChanEst[i],
				rLamMSCSNREst);

			IIR1(vecrSigEstMSC[i], vecrSqMagChanEst[i],
				rLamMSCSNREst);

			/* Calculate MER on MSC cells */
			IIR1(rNoiseEstMSCMER, rCurErrPow, rLamSNREstFast);
		}
	}


	/* Interferer consideration --------------------------------------------- */
	if (bInterfConsid == TRUE)
	{
		for (i = 0; i < iNumCarrier; i++)
		{
			/* Weight the channel estimates with the SNR estimate of the current
			   carrier to consider the higher noise variance caused by
			   interferers */
			(*pvecOutputData)[i].rChan *=
				CalAndBoundSNR(vecrSigEstMSC[i], vecrNoiseEstMSC[i]);
		}
	}
}

void CChannelEstimation::InitInternal(CParameter& ReceiverParam)
{
	/* Get parameters from global struct */
	iScatPilTimeInt = ReceiverParam.iScatPilTimeInt;
	iScatPilFreqInt = ReceiverParam.iScatPilFreqInt;
	iNumIntpFreqPil = ReceiverParam.iNumIntpFreqPil;
	iNumCarrier = ReceiverParam.iNumCarrier;
	iFFTSizeN = ReceiverParam.iFFTSizeN;
	iNumSymPerFrame = ReceiverParam.iNumSymPerFrame;

	/* Length of guard-interval with respect to FFT-size! */
	rGuardSizeFFT = (_REAL) iNumCarrier *
		ReceiverParam.RatioTgTu.iEnum / ReceiverParam.RatioTgTu.iDenom;

	/* If robustness mode D is active, get DC position. This position cannot
	   be "0" since in mode D no 5 kHz mode is defined (see DRM-standard). 
	   Therefore we can also use this variable to get information whether
	   mode D is active or not (by simply write: "if (iDCPos != 0)") */
	if (ReceiverParam.GetWaveMode() == RM_ROBUSTNESS_MODE_D)
	{
		/* Identify DC carrier position */
		for (int i = 0; i < iNumCarrier; i++)
		{
			if (_IsDC(ReceiverParam.matiMapTab[0][i]))
				iDCPos = i;
		}
	}
	else
		iDCPos = 0;

	/* FFT must be longer than "iNumCarrier" because of zero padding effect (
	   not in robustness mode D! -> "iLongLenFreq = iNumCarrier") */
	iLongLenFreq = iNumCarrier + iScatPilFreqInt - 1;

	/* Init vector for received data at pilot positions */
	veccPilots.Init(iNumIntpFreqPil);

	/* Init vector for interpolated pilots */
	veccIntPil.Init(iLongLenFreq);

	/* Init plans for FFT (faster processing of Fft and Ifft commands) */
	FftPlanShort.Init(iNumIntpFreqPil);
	FftPlanLong.Init(iLongLenFreq);

	/* Choose time interpolation method and set pointer to correcponding 
	   object */
	switch (TypeIntTime)
	{
	case TLINEAR:
		pTimeInt = &TimeLinear;
		break;

	case TWIENER:
		pTimeInt = &TimeWiener;
		break;
	}

	/* Init time interpolation interface and set delay for interpolation */
	iLenHistBuff = pTimeInt->Init(ReceiverParam);

	/* Init time synchronization tracking unit */
	TimeSyncTrack.Init(ReceiverParam, iLenHistBuff);

	/* Set channel estimation delay in global struct. This is needed for 
	   simulation */
	ReceiverParam.iChanEstDelay = iLenHistBuff;


	/* Init window for DFT operation for frequency interpolation ------------ */
	/* Init memory */
	vecrDFTWindow.Init(iNumIntpFreqPil);
	vecrDFTwindowInv.Init(iNumCarrier);

	/* Set window coefficients */
	switch (eDFTWindowingMethod)
	{
	case DFT_WIN_RECT:
		vecrDFTWindow = Ones(iNumIntpFreqPil);
		vecrDFTwindowInv = Ones(iNumCarrier);
		break;

	case DFT_WIN_HAMM:
		vecrDFTWindow = Hamming(iNumIntpFreqPil);
		vecrDFTwindowInv = (CReal) 1.0 / Hamming(iNumCarrier);
		break;
	}


	/* Set start index for zero padding in time domain for DFT method */
	iStartZeroPadding = (int) rGuardSizeFFT;
	if (iStartZeroPadding > iNumIntpFreqPil)
		iStartZeroPadding = iNumIntpFreqPil;

	/* Allocate memory for channel estimation */
	veccChanEst.Init(iNumCarrier);
	vecrSqMagChanEst.Init(iNumCarrier);

	/* Allocate memory for history buffer (Matrix) and zero out */
	matcHistory.Init(iLenHistBuff, iNumCarrier,
		_COMPLEX((_REAL) 0.0, (_REAL) 0.0));

	/* After an initialization we do not put out data before the number symbols
	   of the channel estimation delay have been processed */
	iInitCnt = iLenHistBuff - 1;

	/* SNR correction factor for the different SNR estimation types. For the
	   FAC method, the average signal power has to be considered. For the pilot
	   based method, only the SNR on the pilots are evaluated. Therefore, to get
	   the total SNR, a correction has to be applied */
	rSNRFACSigCorrFact = ReceiverParam.rAvPowPerSymbol / CReal(iNumCarrier);
	rSNRTotToPilCorrFact = ReceiverParam.rAvScatPilPow *
		(_REAL) iNumCarrier / ReceiverParam.rAvPowPerSymbol;

	/* Correction factor for transforming the estimated system SNR in the SNR
	   where the noise bandwidth is according to the nominal DRM bandwidth */
	rSNRSysToNomBWCorrFact = ReceiverParam.GetSysToNomBWCorrFact();

	/* Inits for SNR estimation (noise and signal averages) */
	rSignalEst = (_REAL) 0.0;
	rNoiseEst = (_REAL) 0.0;
	rNoiseEstMSCMER = (_REAL) 0.0;
	rSNREstimate = (_REAL) pow(10, INIT_VALUE_SNR_ESTIM_DB / 10);
	vecrNoiseEstMSC.Init(iNumCarrier, (_REAL) 0.0);
	vecrSigEstMSC.Init(iNumCarrier, (_REAL) 0.0);

	/* For SNR estimation initialization */
	iSNREstIniSigAvCnt = 0;
	iSNREstIniNoiseAvCnt = 0;

	/* We only have an initialization phase for SNR estimation method based on
	   the tentative decisions of FAC cells */
	if (TypeSNREst == SNR_FAC)
		bSNRInitPhase = TRUE;
	else
		bSNRInitPhase = FALSE;

	/* 5 DRM frames to start initial SNR estimation after initialization phase
	   of other units */
	iSNREstInitCnt = 5 * iNumSymPerFrame;

	/* Lambda for IIR filter */
	rLamSNREstFast = IIR1Lam(TICONST_SNREST_FAST, (CReal) SOUNDCRD_SAMPLE_RATE /
		ReceiverParam.iSymbolBlockSize);
	rLamSNREstSlow = IIR1Lam(TICONST_SNREST_SLOW, (CReal) SOUNDCRD_SAMPLE_RATE /
		ReceiverParam.iSymbolBlockSize);
	rLamMSCSNREst = IIR1Lam(TICONST_SNREST_MSC, (CReal) SOUNDCRD_SAMPLE_RATE /
		ReceiverParam.iSymbolBlockSize);

	/* Init delay spread length estimation (index) */
	rLenPDSEst = (_REAL) 0.0;

	/* Init history for delay values */
	/* Duration of OFDM symbol */
	const _REAL rTs = (CReal) (ReceiverParam.iFFTSizeN +
		ReceiverParam.iGuardSize) / SOUNDCRD_SAMPLE_RATE;

	iLenDelayHist = (int) (LEN_HIST_DELAY_LOG_FILE_S / rTs);
	vecrDelayHist.Init(iLenDelayHist, (CReal) 0.0);


	/* Inits for Wiener interpolation in frequency direction ---------------- */
	/* Length of wiener filter */
	switch (ReceiverParam.GetWaveMode())
	{
	case RM_ROBUSTNESS_MODE_A:
		iLengthWiener = LEN_WIENER_FILT_FREQ_RMA;
		break;

	case RM_ROBUSTNESS_MODE_B:
		iLengthWiener = LEN_WIENER_FILT_FREQ_RMB;
		break;

	case RM_ROBUSTNESS_MODE_C:
		iLengthWiener = LEN_WIENER_FILT_FREQ_RMC;
		break;

	case RM_ROBUSTNESS_MODE_D:
		iLengthWiener = LEN_WIENER_FILT_FREQ_RMD;
		break;
	}


	/* Inits for wiener filter ---------------------------------------------- */
	/* In frequency direction we can use pilots from both sides for 
	   interpolation */
	iPilOffset = iLengthWiener / 2;

	/* Allocate memory */
	matcFiltFreq.Init(iNumCarrier, iLengthWiener);

	/* Pilot offset table */
	veciPilOffTab.Init(iNumCarrier);

	/* Number of different wiener filters */
	iNumWienerFilt = (iLengthWiener - 1) * iScatPilFreqInt + 1;

	/* Allocate temporary matlib vector for filter coefficients */
	matcWienerFilter.Init(iNumWienerFilt, iLengthWiener);

	/* Distinguish between simulation and regular receiver. When we run a
	   simulation, the parameters are taken from simulation init */
	if (ReceiverParam.eSimType == CParameter::ST_NONE)
	{
		/* Initial Wiener filter. Use initial SNR definition and assume that the
		   PDS ranges from the beginning of the guard-intervall to the end */
		UpdateWienerFiltCoef(pow(10, INIT_VALUE_SNR_WIEN_FREQ_DB / 10),
			(_REAL) ReceiverParam.RatioTgTu.iEnum /
			ReceiverParam.RatioTgTu.iDenom, (CReal) 0.0);
	}
	else
	{
		/* Get simulation SNR on the pilot positions and set PDS to entire
		   guard-interval length */
		UpdateWienerFiltCoef(pow(10, ReceiverParam.GetSysSNRdBPilPos() / 10),
			(_REAL) ReceiverParam.RatioTgTu.iEnum /
			ReceiverParam.RatioTgTu.iDenom, (CReal) 0.0);
	}


	/* Define block-sizes for input and output */
	iInputBlockSize = iNumCarrier;
	iMaxOutputBlockSize = iNumCarrier; 
}

CComplexVector CChannelEstimation::FreqOptimalFilter(int iFreqInt, int iDiff,
													 CReal rSNR,
													 CReal rRatPDSLen,
													 CReal rRatPDSOffs,
													 int iLength)
{
/* 
	We assume that the power delay spread is a rectangle function in the time
	domain (sinc-function in the frequency domain). Length and position of this
	window are adapted according to the current estimated PDS.
*/
	int				i;
	CRealVector		vecrRpp(iLength);
	CRealVector		vecrRhp(iLength);
	CRealVector		vecrH(iLength);
	CComplexVector	veccH(iLength);

	/* Calculation of R_hp, this is the SHIFTED correlation function */
	for (i = 0; i < iLength; i++)
	{
		const int iCurPos = i * iFreqInt - iDiff;

		vecrRhp[i] = Sinc((CReal) iCurPos * rRatPDSLen);
	}

	/* Calculation of R_pp */
	for (i = 0; i < iLength; i++)
	{
		const int iCurPos = i * iFreqInt;

		vecrRpp[i] = Sinc((CReal) iCurPos * rRatPDSLen);
	}

	/* Add SNR at first tap */
	vecrRpp[0] += (CReal) 1.0 / rSNR;

	/* Call levinson algorithm to solve matrix system for optimal solution */
	vecrH = Levinson(vecrRpp, vecrRhp);

	/* Correct the optimal filter coefficients. Shift the rectangular
	   function in the time domain to the correct position (determined by
	   the "rRatPDSOffs") by multiplying in the frequency domain
	   with exp(j w T) */
	for (i = 0; i < iLength; i++)
	{
		const int iCurPos = i * iFreqInt - iDiff;

		const CReal rArgExp =
			(CReal) crPi * iCurPos * (rRatPDSLen + rRatPDSOffs * 2);

		veccH[i] = vecrH[i] * CComplex(Cos(rArgExp), Sin(rArgExp));
	}

	return veccH;