channelestimation.cpp

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

/******************************************************************************\
 * Technische Universitaet Darmstadt, Institut fuer Nachrichtentechnik
 * Copyright (c) 2001
 *
 * Author(s):
 *	Volker Fischer
 *
 * Description:
 *	Channel estimation and equalization
 *
 ******************************************************************************
 *
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 2 of the License, or (at your option) any later 
 * version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT 
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 
 * details.
 *
 * You should have received a copy of the GNU General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 
 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
\******************************************************************************/

#include "ChannelEstimation.h"


/* Implementation *************************************************************/
void CChannelEstimation::ProcessDataInternal(CParameter& ReceiverParam)
{
	int		i, j, k;
	int		iModSymNum;
	_REAL	rCurSNREst;
	_REAL	rOffsPDSEst;

	/* Check if symbol ID index has changed by the synchronization unit. If it
	   has changed, reinit this module */
	if ((*pvecInputData).GetExData().bSymbolIDHasChanged == TRUE)
	{
// FIXME: we loose one OFDM symbol by this call -> slower DRM signal acquisition
		SetInitFlag();
		return;
	}

	/* Move data in history-buffer (from iLenHistBuff - 1 towards 0) */
	for (j = 0; j < iLenHistBuff - 1; j++)
	{
		for (i = 0; i < iNumCarrier; i++)
			matcHistory[j][i] = matcHistory[j + 1][i];
	}

	/* Write new symbol in memory */
	for (i = 0; i < iNumCarrier; i++)
		matcHistory[iLenHistBuff - 1][i] = (*pvecInputData)[i];


	/* Time interpolation *****************************************************/
	/* Get symbol-counter for next symbol. Use the count from the frame 
	   synchronization (in OFDM.cpp). Call estimation routine */
	const _REAL rSNRAftTiInt = 
		pTimeInt->Estimate(pvecInputData, veccPilots, 
						   ReceiverParam.matiMapTab[(*pvecInputData).
						   GetExData().iSymbolID],
						   ReceiverParam.matcPilotCells[(*pvecInputData).
						   GetExData().iSymbolID],
						   /* The channel estimation is based on the pilots so
						      it needs the SNR on the pilots. Do a correction */
						   rSNREstimate * rSNRTotToPilCorrFact);

	/* Debar initialization of channel estimation in time direction */
	if (iInitCnt > 0)
	{
		iInitCnt--;

		/* Do not put out data in initialization phase */
		iOutputBlockSize = 0;

		/* Do not continue */
		return;
	}
	else
		iOutputBlockSize = iNumCarrier; 

	/* Define DC carrier for robustness mode D because there is no pilot */
	if (iDCPos != 0)
		veccPilots[iDCPos] = (CReal) 0.0;


	/* -------------------------------------------------------------------------
	   Use time-interpolated channel estimate for timing synchronization 
	   tracking */
	TimeSyncTrack.Process(ReceiverParam, veccPilots, 
		(*pvecInputData).GetExData().iCurTimeCorr, rLenPDSEst /* out */,
		rOffsPDSEst /* out */);

	/* Store current delay in history */
	vecrDelayHist.AddEnd(rLenPDSEst);


	/* Frequency-interploation ************************************************/
	switch (TypeIntFreq)
	{
	case FLINEAR:
		/**********************************************************************\
		 * Linear interpolation												  *
		\**********************************************************************/
		/* Set first pilot position */
		veccChanEst[0] = veccPilots[0];

		for (j = 0, k = 1; j < iNumCarrier - iScatPilFreqInt;
			j += iScatPilFreqInt, k++)
		{
			/* Set values at second time pilot position in cluster */
			veccChanEst[j + iScatPilFreqInt] = veccPilots[k];

			/* Interpolation cluster */
			for (i = 1; i < iScatPilFreqInt; i++)
			{
				/* E.g.: c(x) = (c_4 - c_0) / 4 * x + c_0 */
				veccChanEst[j + i] =
					(veccChanEst[j + iScatPilFreqInt] - veccChanEst[j]) /
					(_REAL) (iScatPilFreqInt) * (_REAL) i + veccChanEst[j];
			}
		}
		break;

	case FDFTFILTER:
		/**********************************************************************\
		 * DFT based algorithm												  *
		\**********************************************************************/
		/* ---------------------------------------------------------------------
		   Put all pilots at the beginning of the vector. The "real" length of
		   the vector "pcFFTWInput" is longer than the No of pilots, but we 
		   calculate the FFT only over "iNumCarrier / iScatPilFreqInt + 1"
		   values (this is the number of pilot positions) */
		/* Weighting pilots with window */
		veccPilots *= vecrDFTWindow;

		/* Transform in time-domain */
		veccPilots = Ifft(veccPilots, FftPlanShort);

		/* Set values outside a defined bound to zero, zero padding (noise
		   filtering). Copy second half of spectrum at the end of the new vector
		   length and zero out samples between the two parts of the spectrum */
		veccIntPil.Merge(
			/* First part of spectrum */
			veccPilots(1, iStartZeroPadding), 
			/* Zero padding in the middle, length: Total length minus length of
			   the two parts at the beginning and end */
			CComplexVector(Zeros(iLongLenFreq - 2 * iStartZeroPadding), 
			Zeros(iLongLenFreq - 2 * iStartZeroPadding)), 
			/* Set the second part of the actual spectrum at the end of the new
			   vector */
			veccPilots(iNumIntpFreqPil - iStartZeroPadding + 1, 
			iNumIntpFreqPil));

		/* Transform back in frequency-domain */
		veccIntPil = Fft(veccIntPil, FftPlanLong);

		/* Remove weighting with DFT window by inverse multiplication */
		veccChanEst = veccIntPil(1, iNumCarrier) * vecrDFTwindowInv;
		break;

	case FWIENER:
		/**********************************************************************\
		 * Wiener filter													   *
		\**********************************************************************/
		/* Wiener filter update --------------------------------------------- */
		/* Do not update filter in case of simulation */
		if (ReceiverParam.eSimType == CParameter::ST_NONE)
		{
			/* Update Wiener filter each OFDM symbol. Use current estimates */
			UpdateWienerFiltCoef(rSNRAftTiInt, rLenPDSEst / iNumCarrier,
				rOffsPDSEst / iNumCarrier);
		}


		/* Actual Wiener interpolation (FIR filtering) ---------------------- */
		/* FIR filter of the pilots with filter taps. We need to filter the
		   pilot positions as well to improve the SNR estimation (which 
		   follows this procedure) */
		for (j = 0; j < iNumCarrier; j++)
		{
// TODO: Do only calculate channel estimation for data cells, not for pilot
// cells (exeption: if we want to use SNR estimation based on pilots, we also
// need Wiener on these cells!)
			/* Convolution */
			veccChanEst[j] = _COMPLEX((_REAL) 0.0, (_REAL) 0.0);

			for (i = 0; i < iLengthWiener; i++)
			{
				veccChanEst[j] +=
					matcFiltFreq[j][i] * veccPilots[veciPilOffTab[j] + i];
			}
		}
		break;
	}


	/* Equalize the output vector ------------------------------------------- */
	/* Calculate squared magnitude of channel estimation */
	vecrSqMagChanEst = SqMag(veccChanEst);

	/* Write to output vector. Take oldest symbol of history for output. Also,
	   ship the channel state at a certain cell */
	for (i = 0; i < iNumCarrier; i++)
	{
		(*pvecOutputData)[i].cSig = matcHistory[0][i] / veccChanEst[i];

#ifdef USE_MAX_LOG_MAP
		/* In case of MAP we need the squared magnitude for the calculation of
		   the metric */
		(*pvecOutputData)[i].rChan = vecrSqMagChanEst[i];
#else
		/* In case of hard-desicions, we need the magnitude of the channel for
		   the calculation of the metric */
		(*pvecOutputData)[i].rChan = sqrt(vecrSqMagChanEst[i]);
#endif
	}


	/* -------------------------------------------------------------------------
	   Calculate symbol ID of the current output block and set parameter */
	(*pvecOutputData).GetExData().iSymbolID = 
		(*pvecInputData).GetExData().iSymbolID - iLenHistBuff + 1;


	/* SNR estimation ------------------------------------------------------- */
	/* Modified symbol ID, check range {0, ..., iNumSymPerFrame} */
	iModSymNum = (*pvecOutputData).GetExData().iSymbolID;

	while (iModSymNum < 0)
		iModSymNum += iNumSymPerFrame;

	/* Two different types of SNR estimation are available */
	switch (TypeSNREst)
	{
	case SNR_PIL:
		/* Use estimated channel and compare it to the received pilots. This
		   estimation works only if the channel estimation was successful */
		for (i = 0; i < iNumCarrier; i++)
		{
			/* Identify pilot positions. Use MODIFIED "iSymbolID" (See lines
			   above) */
			if (_IsScatPil(ReceiverParam.matiMapTab[iModSymNum][i]))
			{
				/* We assume that the channel estimation in "veccChanEst" is
				   noise free (e.g., the wiener interpolation does noise
				   reduction). Thus, we have an estimate of the received signal
				   power \hat{r} = s * \hat{h}_{wiener} */
				const _COMPLEX cModChanEst = veccChanEst[i] *
					ReceiverParam.matcPilotCells[iModSymNum][i];


				/* Calculate and average noise and signal estimates --------- */
				/* The noise estimation is difference between the noise reduced
				   signal and the noisy received signal
				   \tilde{n} = \hat{r} - r */
				IIR1(rNoiseEst, SqMag(matcHistory[0][i] - cModChanEst),
					rLamSNREstFast);

				/* The received signal power estimation is just \hat{r} */
				IIR1(rSignalEst, SqMag(cModChanEst), rLamSNREstFast);

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

				/* Average the SNR with a two sided recursion. Apply correction
				   factor, too */
				IIR1TwoSided(rSNREstimate, rCurSNREst / rSNRTotToPilCorrFact,
					rLamSNREstFast,	rLamSNREstSlow);
			}
		}
		break;

	case SNR_FAC:
		/* SNR estimation based on FAC cells and hard decisions */
		/* SNR estimation with initialization */
		if (iSNREstInitCnt > 0)
		{
			for (i = 0; i < iNumCarrier; i++)
			{
				/* Only use the last frame of the initialization phase for
				   initial SNR estimation to debar initialization phase of
				   synchronization and channel estimation units */
				if (iSNREstInitCnt < iNumSymPerFrame * iNumCarrier)
				{
					const int iCurCellFlag =
						ReceiverParam.matiMapTab[iModSymNum][i];

					/* Initial signal estimate. Use channel estimation from all
					   data and pilot cells. Apply averaging */
					if ((_IsData(iCurCellFlag)) || (_IsPilot(iCurCellFlag)))
					{
						/* Signal estimation */
						rSignalEst += vecrSqMagChanEst[i];
						iSNREstIniSigAvCnt++;
					}

					/* Noise estimation from all data cells from tentative
					   decisions */
					if (_IsFAC(iCurCellFlag)) /* FAC cell */
					{
						rNoiseEst += vecrSqMagChanEst[i] *
							SqMag(MinDist4QAM((*pvecOutputData)[i].cSig));

						iSNREstIniNoiseAvCnt++;
					}
				}
			}

			iSNREstInitCnt--;
		}
		else
		{
			/* Only right after initialization phase apply initial SNR
			   value */
			if (bSNRInitPhase == TRUE)
			{
				/* Normalize average */
				rSignalEst /= iSNREstIniSigAvCnt;
				rNoiseEst /= iSNREstIniNoiseAvCnt;

				bSNRInitPhase = FALSE;
			}

			for (i = 0; i < iNumCarrier; i++)
			{
				/* Only use FAC cells for this SNR estimation method */
				if (_IsFAC(ReceiverParam.matiMapTab[iModSymNum][i]))
				{
					/* Get tentative decision for this FAC QAM symbol. FAC is
					   always 4-QAM. Calculate all distances to the four
					   possible constellation points of a 4-QAM and use the
					   squared result of the returned distance vector */
					const CReal rCurErrPow =
						SqMag(MinDist4QAM((*pvecOutputData)[i].cSig));

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