timewiener.cpp

Dream.exe soft source (Visual C++)

	else
		iNumTapsSigEst = NUM_TAPS_USED4SIGMA_EST;

	/* Init vector for estimation of the correlation function in time direction
	   (IIR average) */
	veccTiCorrEst.Init(iNumTapsSigEst, (CReal) 0.0);

	/* Init time constant for IIR filter for averaging correlation estimation.
	   Consider averaging over frequency axis, too. Pilots in frequency
	   direction are "iScatPilTimeInt * iScatPilFreqInt" apart */
	const int iNumPilOneOFDMSym = iNumIntpFreqPil / iScatPilTimeInt;
	rLamTiCorrAv = IIR1Lam(TICONST_TI_CORREL_EST * iNumPilOneOFDMSym,
		(CReal) SOUNDCRD_SAMPLE_RATE / ReceiverParam.iSymbolBlockSize);

	/* Init update counter for Wiener filter update. We immediatly use the
	   filtered result although right at the beginning there is no averaging.
	   But sine the estimation usually starts with higher values and goes down
	   to the correct one, this should not be critical */
	iUpCntWienFilt = iNumSymPerFrame;

	/* Init averaging of SNR values */
	rAvSNR = (_REAL) 0.0;
	iAvSNRCnt = 0;


	/* Allocate memory for filter phases (Matrix) */
	matrFiltTime.Init(iNumFiltPhasTi, iLengthWiener);

	/* Length of the timing correction history buffer */
	iLenTiCorrHist = iLengthWiener * iNumFiltPhasTi;

	/* Init timing correction history with zeros */
	vecTiCorrHist.Init(iLenTiCorrHist, 0);


	/* Calculate optimal filter --------------------------------------------- */
	_REAL rSNR;

	/* Distinguish between simulation and regular receiver. When we run a
	   simulation, the parameters are taken from simulation init */
	if (ReceiverParam.eSimType == CParameter::ST_NONE)
	{
		/* Init SNR value */
		rSNR = pow(10, INIT_VALUE_SNR_WIEN_TIME_DB / 10);

		/* Init sigma with a large value. This make the acquisition more
		   robust in case of a large sample frequency offset. But we get more
		   aliasing in the time domain and this could make the timing unit
		   perform worse. Therefore, this is only a trade-off */
		rSigma = rSigmaMax;
	}
	else
	{
		/* Get SNR on the pilot positions */
		rSNR = pow(10, ReceiverParam.GetSysSNRdBPilPos() / 10);
	
		/* Sigma from channel profiles */
		switch (ReceiverParam.iDRMChannelNum)
		{
		case 1:
		case 2:
			rSigma = LOW_BOUND_SIGMA;
			break;

		case 4:
			rSigma = 1.0 / 2;
			break;

		case 3:
			rSigma = 1.2 / 2;
			break;

		case 5:
			rSigma = 2.0 / 2;
			break;

		case 8:
			rSigma = ReceiverParam.iChan8Doppler;
			break;

		default: /* Including channel number 6 */
			rSigma = rSigmaMax / 2;
			break;
		}

		/* Reset flag to inhibit parameter adaptation */
		bTracking = FALSE;
	}

	/* Calculate initialization wiener filter taps and init MMSE */
	rMMSE = UpdateFilterCoef(rSNR, rSigma);

	/* Return delay of channel equalization */
	return iLenHistBuff;
}

void CTimeWiener::GenFiltPhaseTable(CMatrix<int>& matiMapTab,
									const int iNumCarrier,
									const int iNumSymPerFrame,
									const int iScatPilTimeInt)
{
	/* Init matrix */
	matiFiltPhaseTable.Init(iNumCarrier, iNumSymPerFrame);

	/* Get the index of first symbol in a super-frame on where the first cell
	   (carrier-index = 0) is a pilot. This is needed for determining the
	   correct filter phase for the convolution */
	int iFirstSymWithPilot = 0;
	while (!_IsScatPil(matiMapTab[iFirstSymWithPilot][0]))
		iFirstSymWithPilot++;

	for (int i = 0; i < iNumCarrier; i++)
	{
		for (int j = 0; j < iNumSymPerFrame; j++)
		{
			/* Calculate filter phases for Wiener filter for each OFDM cell in
			   a DRM frame */
			matiFiltPhaseTable[i][j] = (iScatPilTimeInt -
				(iNumSymPerFrame - j + iFirstSymWithPilot + i) %
				iScatPilTimeInt) % iScatPilTimeInt;
		}
	}
}

_REAL CTimeWiener::UpdateFilterCoef(const _REAL rNewSNR, const _REAL rNewSigma)
{
	/* Calculate MMSE for wiener filtering for all phases and average */
	_REAL rMMSE = (_REAL) 0.0;

	/* One filter for all possible filter phases */
	for (int j = 0; j < iNumFiltPhasTi; j++)
	{
		/* We have to define the dependency between the difference between the
		   current pilot to the observed symbol in the history buffer and the
		   indizes of the FiltTime array. Definition:
		   Largest distance = index zero, index increases to smaller
		   distances */
		const int iCurrDiffPhase = -(iLenHistBuff - j - 1);

		/* Calculate filter phase and average MMSE */
		rMMSE += TimeOptimalFilter(matrFiltTime[j], iScatPilTimeInt,
			iCurrDiffPhase,	rNewSNR, rNewSigma, rTs, iLengthWiener);
	}

	/* Normalize averaged MMSE */
	rMMSE /= iNumFiltPhasTi;

	return rMMSE;
}

CReal CTimeWiener::TimeOptimalFilter(CRealVector& vecrTaps, const int iTimeInt,
									 const int iDiff, const CReal rNewSNR,
									 const CReal rNewSigma, const CReal rTs,
									 const int iLength)
{
	int i;

	CRealVector vecrRpp(iLength);
	CRealVector vecrRhp(iLength);

	/* Factor for the argument of the exponetial function to generate the
	   correlation function */
	const CReal rFactorArgExp = 
		(CReal) -2.0 * crPi * crPi * rTs * rTs * rNewSigma * rNewSigma;

	/* Doppler-spectrum for short-wave channel is Gaussian
	   (Calculation of R_hp!) */
	for (i = 0; i < iLength; i++)
	{
		const int iCurPos = i * iTimeInt + iDiff;

		vecrRhp[i] = exp(rFactorArgExp * iCurPos * iCurPos);
	}

	/* Doppler-spectrum for short-wave channel is Gaussian
	   (Calculation of R_pp!) */
	for (i = 0; i < iLength; i++)
	{
		const int iCurPos = i * iTimeInt;

		vecrRpp[i] = exp(rFactorArgExp * iCurPos * iCurPos);
	}

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

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

	/* Return MMSE for the current wiener filter */
	return (CReal) 1.0 - Sum(vecrRhp * vecrTaps);
}

CReal CTimeWiener::ModLinRegr(const CComplexVector& veccCorrEst)
{
	/* Modified linear regression to estimate the "sigma" of the Gaussian
	   correlation function */
	/* Get vector length */
	const int iVecLen = Size(veccCorrEst);

	/* Init vectors and variables */
	CReal		rSigmaRet;
	CRealVector Tau(iVecLen);
	CRealVector Z(iVecLen);
	CRealVector W(iVecLen);
	CRealVector Wmrem(iVecLen);
	CReal		Wm, Zm;
	CReal		A1;

	/* Generate the tau vector */
	for (int i = 0; i < iVecLen; i++)
		Tau[i] = (CReal) (i * iScatPilTimeInt);

	/* Linearize acf equation:  y = a * exp(-b * x ^ 2)
	   z = ln(y); w = x ^ 2
	   -> z = a0 + a1 * w */
	Z = Log(Abs(veccCorrEst));

	W = Tau * Tau;

	Wm = Mean(W);
	Zm = Mean(Z);

	Wmrem = W - Wm; /* Remove mean of W */

	A1 = Sum(Wmrem * (Z - Zm)) / Sum(Wmrem * Wmrem);

	/* Final sigma calculation from estimation and assumed Gaussian model */
	rSigmaRet = (CReal) 0.5 / crPi * Sqrt((CReal) -2.0 * A1) / rTs;

	/* Bound estimated sigma value */
	if (rSigmaRet > rSigmaMax)
		rSigmaRet = rSigmaMax;
	if (rSigmaRet < LOW_BOUND_SIGMA)
		rSigmaRet = LOW_BOUND_SIGMA;

	return rSigmaRet;
}