timesynctrack.cpp

数字音频广播中的同步算法在pc上的实现 主要运用了导频、最大似然估计等算法结合

/******************************************************************************\
 * Technische Universitaet Darmstadt, Institut fuer Nachrichtentechnik
 * Copyright (c) 2001
 *
 * Author(s):
 *	Volker Fischer
 *
 * Description:
 *	Time synchronization tracking using information of scattered pilots
 *
 *	Algorithm proposed by Baoguo Yang in "Timing Recovery for OFDM
 *	Transmission", IEEE November 2000
 *
 ******************************************************************************
 *
 * 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 "TimeSyncTrack.h"


/* Implementation *************************************************************/
void CTimeSyncTrack::Process(CParameter& Parameter,
							 CComplexVector& veccChanEst, int iNewTiCorr,
							 _REAL& rLenPDS, _REAL& rOffsPDS)
{
	int			i, j;
	int			iIntShiftVal;
	int			iFirstPathDelay;
	CReal		rPeakBound;
	CReal		rPropGain;
	CReal		rCurEnergy;
	CReal		rWinEnergy;
	CReal		rMaxWinEnergy;
	_BOOLEAN	bDelayFound;
	_BOOLEAN	bPDSResultFound;

	/* Rotate the averaged PDP to follow the time shifts -------------------- */
	/* Update timing correction history (shift register) */
	vecTiCorrHist.AddEnd(iNewTiCorr);

	/* Calculate the actual shift of the timing correction. Since we do the
	   timing correction at the sound card sample rate (48 kHz) and the
	   estimated impulse response has a different sample rate (since the
	   spectrum is only one little part of the sound card frequency range)
	   we have to correct the timing correction by a certain bandwidth factor */
	const CReal rActShiftTiCor = rFracPartTiCor -
		(_REAL) vecTiCorrHist[0] * iNumCarrier / iDFTSize;

	/* Integer part of shift */
	const int iIntPartTiCorr = (int) Round(rActShiftTiCor);

	/* Extract the fractional part since we can only correct integer timing
	   shifts */
	rFracPartTiCor = rActShiftTiCor - iIntPartTiCorr;

	/* Shift the values in the vector storing the averaged impulse response. We
	   have to consider two cases for shifting (left and right shift) */
	if (rActShiftTiCor < 0)
		iIntShiftVal = iIntPartTiCorr + iNumIntpFreqPil;
	else
		iIntShiftVal = iIntPartTiCorr;

	/* If new correction is out of range, do not apply rotation */
	if ((iIntShiftVal > 0) && (iIntShiftVal < iNumIntpFreqPil))
	{
		/* Actual rotation of vector */
		vecrAvPoDeSp.Merge(vecrAvPoDeSp(iIntShiftVal + 1, iNumIntpFreqPil),
			vecrAvPoDeSp(1, iIntShiftVal));
	}


	/* New estimate for impulse response ------------------------------------ */
	/* Apply hamming window, Eq (15) */
	veccPilots = veccChanEst * vecrHammingWindow;

	/* Transform in time-domain to get an estimate for the delay power profile,
	   Eq (15) */
	veccPilots = Ifft(veccPilots, FftPlan);

	/* Average result, Eq (16) (Should be a moving average function, for
	   simplicity we have chosen an IIR filter here) */
	IIR1(vecrAvPoDeSp, SqMag(veccPilots), rLamAvPDS);

	/* Rotate the averaged result vector to put the earlier peaks
	   (which can also detected in a certain amount) at the beginning of
	   the vector */
	vecrAvPoDeSpRot.Merge(vecrAvPoDeSp(iStPoRot, iNumIntpFreqPil),
		vecrAvPoDeSp(1, iStPoRot - 1));


	/* Different timing algorithms ------------------------------------------ */
	switch (TypeTiSyncTrac)
	{
	case TSFIRSTPEAK:
		/* Detect first peak algorithm proposed by Baoguo Yang */
		/* Lower and higher bound */
		rBoundHigher = Max(vecrAvPoDeSpRot) * rConst1;
		rBoundLower = Min(vecrAvPoDeSpRot) * rConst2;

		/* Calculate the peak bound, Eq (19) */
		rPeakBound = Max(rBoundHigher, rBoundLower);

		/* Get final estimate, Eq (18) */
		bDelayFound = FALSE; /* Init flag */
		for (i = 0; i < iNumIntpFreqPil - 1; i++)
		{
			/* We are only interested in the first peak */
			if (bDelayFound == FALSE)
			{
				if ((vecrAvPoDeSpRot[i] > vecrAvPoDeSpRot[i + 1]) &&
					(vecrAvPoDeSpRot[i] > rPeakBound))
				{
					/* The first peak was found, store index */
					iFirstPathDelay = i;

					/* Set the flag */
					bDelayFound = TRUE;
				}
			}
		}
		break;

	case TSENERGY:
		/* Determin position of window with maximum energy in guard-interval.
		   A window with the size of the guard-interval is moved over the entire
		   profile and the energy inside this window is calculated. The window
		   position which maximises this energy is taken as the new timing
		   position */
		rMaxWinEnergy = (CReal) 0.0;
		iFirstPathDelay = 0;
		for (i = 0; i < iNumIntpFreqPil - 1 - rGuardSizeFFT; i++)
		{
			rWinEnergy = (CReal) 0.0;

			/* Energy IN the guard-interval */
			for (j = 0; j < rGuardSizeFFT; j++)
				rWinEnergy += vecrAvPoDeSpRot[i + j];

			/* Get maximum */
			if (rWinEnergy > rMaxWinEnergy)
			{
				rMaxWinEnergy = rWinEnergy;
				iFirstPathDelay = i;
			}
		}

		/* We always have a valid measurement, set flag */
		bDelayFound = TRUE;
		break;
	}


	/* Only apply timing correction if search was successful and tracking is
	   activated */
	if ((bDelayFound == TRUE) && (bTiSyncTracking == TRUE))
	{
		/* Consider the rotation introduced for earlier peaks in path delay.
		   Since the "iStPoRot" is the position of the beginning of the block
		   at the end for cutting out, "iNumIntpFreqPil" must be substracted.
		   (Actually, a part of the following line should be look like this:
		   "iStPoRot - 1 - iNumIntpFreqPil + 1" but the "- 1 + 1" compensate
		   each other) */
		iFirstPathDelay += iStPoRot - iNumIntpFreqPil - iTargetTimingPos - 1;


		/* Correct timing offset -------------------------------------------- */
		/* Final offset is target position in comparision to the estimated first
		   path delay. Since we have a delay from the channel estimation, the
		   previous correction is subtracted "- vecrNewMeasHist[0]". If the
		   "real" correction arrives after the delay, this correction is
		   compensated. The length of the history buffer (vecrNewMeasHist) must
		   be equal to the delay of the channel estimation.
		   The corrections must be quantized to the upsampled output sample
		   rate ("* iDFTSize / iNumCarrier") */
		const CReal rTiOffset = (CReal) -iFirstPathDelay *
			iDFTSize / iNumCarrier - veciNewMeasHist[0];

		/* Different controlling parameters for different types of tracking */
		switch (TypeTiSyncTrac)
		{
		case TSFIRSTPEAK:
			/* Adapt the linear control parameter to the region, where the peak
			   was found. The region left of the desired timing position is
			   critical, because we immediately get ISI if a peak appers here.
			   Therefore we apply fast correction here. At the other positions,
			   we smooth the controlling to improve the immunity against false
			   peaks */
			if (rTiOffset > 0)
				rPropGain = CONT_PROP_BEFORE_GUARD_INT;
			else
				rPropGain = CONT_PROP_IN_GUARD_INT;
			break;

		case TSENERGY:
			rPropGain = CONT_PROP_ENERGY_METHOD;
			break;
		}

		/* In case of sample rate offset acquisition phase, use faster timing
		   corrections */
		if (bSamRaOffsAcqu == TRUE)
			rPropGain *= 2;

		/* Apply proportional control and fix result to sample grid */
		const CReal rCurCorrValue = rTiOffset * rPropGain + rFracPartContr;
		const int iContrTiOffs = (int) Fix(rCurCorrValue);

		/* Calculate new fractional part of controlling */
		rFracPartContr = rCurCorrValue - iContrTiOffs;

		/* Manage correction history */
		veciNewMeasHist.AddEnd(0);
		for (i = 0; i < iSymDelay - 1; i++)
			veciNewMeasHist[i] += iContrTiOffs;

		/* Apply correction */
		Parameter.iTimingOffsTrack = -iContrTiOffs;
	}


	/* Sample rate offset estimation ---------------------------------------- */
	/* This sample rate offset estimation is based on the movement of the
	   estimated PDS with time. The movement per symbol (or a number of symbols)
	   is proportional to the sample rate offset. It has to be considered the
	   PDS shiftings of the timing correction unit ("rActShiftTiCor" can be used
	   for that). The procedere is to detect the maximum peak in the PDS and use
	   this as a reference, assuming tha delay of this peak is not changing. The
	   problem is when another peak get higher than this due to fading. In this
	   case we must adjust the history to this new peak (the new reference) */
	int		iMaxInd;
	CReal	rMax;

	/* Find index of maximum peak in PDS estimation. This is our reference
	   for this estimation method */
	Max(rMax, iMaxInd, vecrAvPoDeSpRot);

	/* Integration of timing corrections
	   FIXME: Check for overrun of "iIntegTiCorrections" variable! */
	iIntegTiCorrections += (long int) iIntPartTiCorr;

	/* We need to consider the timing corrections done by the timing unit. What
	   we want to estimate is only the real movement of the detected maximum
	   peak */
	const int iCurRes = iIntegTiCorrections + iMaxInd;
	veciSRTiCorrHist.AddEnd(iCurRes);

	/* We assumed that the detected peak is always the same peak in the actual
	   PDS. But due to fading the maximum can change to a different peak. In
	   this case the estimation would be wrong. We try to detect the detection
	   of a different peak by defining a maximum sample rate change. The sample
	   rate offset is very likely to be very constant since usually crystal
	   oscialltors are used. Thus, if a larger change of sample rate offset
	   happens, we assume that the maximum peak has changed */
	const int iNewDiff = veciSRTiCorrHist[iLenCorrectionHist - 2] - iCurRes;

	/* If change is larger than 2, it is most likely that a new peak was chosen
	   by the maximum function. Also, if the sign has changed of the difference
	   (and it is not zero), we also say that a new peak was selected */
	if ((abs(iNewDiff) > 2) ||
		((Sign(iOldNonZeroDiff) != Sign(iNewDiff)) && (iNewDiff != 0)))
	{
		/* Correct the complete history to the new reference peak. Reference
		   peak was already added, therefore do not use last element */
		for (i = 0; i < iLenCorrectionHist - 1; i++)
			veciSRTiCorrHist[i] -= iNewDiff;
	}

	/* Store old difference if it is not zero */
	if (iNewDiff != 0)
		iOldNonZeroDiff = iNewDiff;


	/* Check, if we are in acquisition phase */
	if (iResOffsetAcquCnt > 0)
	{
		/* Acquisition phase */
		iResOffsetAcquCnt--;
	}
	else
	{
		/* Apply the result from acquisition only once */
		if (bSamRaOffsAcqu == TRUE)
		{
			/* End of acquisition phase */
			bSamRaOffsAcqu = FALSE;

			/* Set sample rate offset to initial estimate. We consider the
			   initialization phase of channel estimation by "iSymDelay" */
			CReal rInitSamOffset = GetSamOffHz(iCurRes - veciSRTiCorrHist[
				iLenCorrectionHist - (iResOffAcqCntMax - iSymDelay)],
				iResOffAcqCntMax - iSymDelay - 1);

#ifndef USE_SAMOFFS_TRACK_FRE_PIL
			/* Apply initial sample rate offset estimation */
			Parameter.rResampleOffset -= rInitSamOffset;
#endif

			/* Reset estimation history (init with zeros) since the sample
			   rate offset was changed */
			veciSRTiCorrHist.Init(iLenCorrectionHist, 0);
			iIntegTiCorrections = 0;
		}
		else
		{
			/* Tracking phase */
			/* Get actual sample rate offset in Hertz */
			const CReal rSamOffset = GetSamOffHz(iCurRes - veciSRTiCorrHist[0],
				iLenCorrectionHist - 1);

#ifndef USE_SAMOFFS_TRACK_FRE_PIL
			/* Apply result from sample rate offset estimation */