timewiener.cpp

Dream.exe soft source (Visual C++)

/******************************************************************************\
 * Technische Universitaet Darmstadt, Institut fuer Nachrichtentechnik
 * Copyright (c) 2001
 *
 * Author(s):
 *	Volker Fischer
 *
 * Description:
 *	Wiener filter in time direction for channel estimation
 *
 ******************************************************************************
 *
 * 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 "TimeWiener.h"


/* Implementation *************************************************************/
_REAL CTimeWiener::Estimate(CVectorEx<_COMPLEX>* pvecInputData,
						    CComplexVector& veccOutputData,
						    CVector<int>& veciMapTab,
						    CVector<_COMPLEX>& veccPilotCells, _REAL rSNR)
{
	int			j, i;
	int			iPiHiIdx;
	int			iTimeDiffNew;
	_COMPLEX	cNewPilot;

	/* Timing correction history -------------------------------------------- */
	/* Shift old vaules and add a "0" at the beginning of the vector */
	vecTiCorrHist.AddBegin(0);

	/* Add new one to all history values except of the current one */
	for (i = 1; i < iLenTiCorrHist; i++)
		vecTiCorrHist[i] += (*pvecInputData).GetExData().iCurTimeCorr;


	/* Update histories for channel estimates at the pilot positions -------- */
	for (i = 0, iPiHiIdx = 0; i < iNumCarrier; i += iScatPilFreqInt, iPiHiIdx++)
	{
		/* Identify and calculate transfer function at the pilot positions */
		if (_IsScatPil(veciMapTab[i]))
		{
			/* Save channel estimates at the pilot positions for each carrier
			   Move old estimates and put new value. Use reversed order to
			   prepare vector for convolution */
			for (j = iLengthWiener - 1; j > 0; j--)
			{
				matcChanAtPilPos[j][iPiHiIdx] =
					matcChanAtPilPos[j - 1][iPiHiIdx];
			}

			/* Add new channel estimate: h = r / s, h: transfer function of the
			   channel, r: received signal, s: transmitted signal */
			matcChanAtPilPos[0][iPiHiIdx] =
				(*pvecInputData)[i] / veccPilotCells[i];


			/* Estimation of the channel correlation function --------------- */
			/* We calcuate the estimation for one symbol first and average this
			   result */
			for (j = 0; j < iNumTapsSigEst; j++)
			{
				/* Correct pilot information for phase rotation */
				iTimeDiffNew = vecTiCorrHist[iScatPilTimeInt * j];
				cNewPilot = Rotate(matcChanAtPilPos[j][iPiHiIdx], i,
					iTimeDiffNew);

				/* Use IIR filtering for averaging */
				IIR1(veccTiCorrEst[j],
					Conj(matcChanAtPilPos[0][iPiHiIdx]) * cNewPilot,
					rLamTiCorrAv);
			}
		}


		/* Wiener interpolation and filtering ------------------------------- */
		/* This check is for robustness mode D since "iScatPilFreqInt" is "1"
		   in this case it would include the DC carrier in the for-loop */
		if (!_IsDC(veciMapTab[i]))
		{
			/* Read current filter phase from table */
			const int iCurrFiltPhase = matiFiltPhaseTable[iPiHiIdx]
				[(*pvecInputData).GetExData().iSymbolID];

			/* Convolution with one phase of the optimal filter */
			/* Init sum */
			_COMPLEX cCurChanEst = _COMPLEX((_REAL) 0.0, (_REAL) 0.0);
			for (j = 0; j < iLengthWiener; j++)
			{
				/* We need to correct pilots due to timing corrections */
				/* Calculate timing difference */
				iTimeDiffNew =
					vecTiCorrHist[j * iScatPilTimeInt + iCurrFiltPhase] -
					vecTiCorrHist[iLenHistBuff - 1];

				/* Correct pilot information for phase rotation */
				cNewPilot =
					Rotate(matcChanAtPilPos[j][iPiHiIdx], i, iTimeDiffNew);

				/* Actual convolution with filter phase */
				cCurChanEst += cNewPilot * matrFiltTime[iCurrFiltPhase][j];
			}

			/* Copy channel estimation from current symbol in output buffer */
			veccOutputData[iPiHiIdx] = cCurChanEst;
		}
	}


	/* Update sigma estimation ---------------------------------------------- */
	if (bTracking == TRUE)
	{
		/* Update filter coefficients once in one DRM frame */
		if (iUpCntWienFilt > 0)
		{
			iUpCntWienFilt--;

			/* Average estimated SNR values */
			rAvSNR += rSNR;
			iAvSNRCnt++;
		}
		else
		{
			/* Actual estimation of sigma */
			rSigma = ModLinRegr(veccTiCorrEst);

			/* Use overestimated sigma for filter update */
			const _REAL rSigOverEst = rSigma * SIGMA_OVERESTIMATION_FACT;

			/* Update the wiener filter, use averaged SNR */
			if (rSigOverEst < rSigmaMax)
				rMMSE = UpdateFilterCoef(rAvSNR / iAvSNRCnt, rSigOverEst);
			else
				rMMSE = UpdateFilterCoef(rAvSNR / iAvSNRCnt, rSigmaMax);

			/* If no SNR improvent is achieved by the optimal filter, use
			   SNR estimation for MMSE */
			_REAL rNewSNR = (_REAL) 1.0 / rMMSE;
			if (rNewSNR < rSNR)
				rMMSE = (_REAL) 1.0 / rSNR;

			/* Reset counter and sum (for SNR) */
			iUpCntWienFilt = iNumSymPerFrame;
			iAvSNRCnt = 0;
			rAvSNR = (_REAL) 0.0;
		}
	}

	/* Return the SNR improvement by Wiener interpolation in time direction */
	return (_REAL) 1.0 / rMMSE;
}

int CTimeWiener::Init(CParameter& ReceiverParam)
{
	/* Init base class, must be at the beginning of this init! */
	CPilotModiClass::InitRot(ReceiverParam);

	/* Set local parameters */
	iNumCarrier = ReceiverParam.iNumCarrier;
	iScatPilTimeInt = ReceiverParam.iScatPilTimeInt;
	iScatPilFreqInt = ReceiverParam.iScatPilFreqInt;
	iNumSymPerFrame = ReceiverParam.iNumSymPerFrame;
	const int iNumIntpFreqPil = ReceiverParam.iNumIntpFreqPil;

	/* Generate filter phase table for Wiener filter */
	GenFiltPhaseTable(ReceiverParam.matiMapTab, iNumCarrier, iNumSymPerFrame,
		iScatPilTimeInt);

	/* Init length of filter and maximum value of sigma (doppler) */
	switch (ReceiverParam.GetWaveMode())
	{
	case RM_ROBUSTNESS_MODE_A:
		iLengthWiener = LEN_WIENER_FILT_TIME_RMA;
		rSigmaMax = MAX_SIGMA_RMA;
		break;

	case RM_ROBUSTNESS_MODE_B:
		iLengthWiener = LEN_WIENER_FILT_TIME_RMB;
		rSigmaMax = MAX_SIGMA_RMB;
		break;

	case RM_ROBUSTNESS_MODE_C:
		iLengthWiener = LEN_WIENER_FILT_TIME_RMC;
		rSigmaMax = MAX_SIGMA_RMC;
		break;

	case RM_ROBUSTNESS_MODE_D:
		iLengthWiener = LEN_WIENER_FILT_TIME_RMD;
		rSigmaMax = MAX_SIGMA_RMD;
		break;
	}

	/* Set delay of this channel estimation type. The longer the delay is, the
	   more "acausal" pilots can be used for interpolation. We use the same
	   amount of causal and acausal filter taps here. Make sure that we get
	   R_hp's which have the most energy collected:
	   L = Np * TiPi - TiPi + 1 is the total number of cells which span our
	   interpolation when we set a pilot on the left-most and right-most
	   Ceil(L / 2) is the middle of the range, now we only have to consider
	   half of the TiPi which is Floor(TiPi / 2) */
	const int iSymDelyChanEst = (int) Ceil((CReal) (
		iLengthWiener * iScatPilTimeInt - iScatPilTimeInt + 1) / 2) +
		(int) Floor((CReal) iScatPilTimeInt / 2) - 1;

	/* Set number of phases for wiener filter */
	iNumFiltPhasTi = iScatPilTimeInt;

	/* Set length of history-buffer */
	iLenHistBuff = iSymDelyChanEst + 1;

	/* Duration of useful part plus-guard interval */
	rTs = (_REAL) ReceiverParam.iSymbolBlockSize / SOUNDCRD_SAMPLE_RATE;

	/* Total number of interpolated pilots in frequency direction. We have to
	   consider the last pilot at the end ("+ 1") */
	const int iTotNumPiFreqDir = iNumCarrier / iScatPilFreqInt + 1;

	/* Allocate memory for Channel at pilot positions (matrix) and init with
	   ones */
	matcChanAtPilPos.Init(iLengthWiener, iTotNumPiFreqDir,
		_COMPLEX((_REAL) 1.0, (_REAL) 0.0));

	/* Set number of taps for sigma estimation */
	if (iLengthWiener < NUM_TAPS_USED4SIGMA_EST)
		iNumTapsSigEst = iLengthWiener;