timesynctrack.cpp

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

			Parameter.rResampleOffset -= CONTR_SAMP_OFF_INT_FTI * rSamOffset;
#endif
		}
	}


	/* Delay spread length estimation --------------------------------------- */
	/* Estimate the noise energy using the minimum statistic. We assume that
	   the noise variance is equal on all samples of the impulse response.
	   Therefore we subtract the variance on each sample. The total estimated
	   signal energy is the total energy minus the noise energy */
	/* Calculate total energy */
	const CReal rTotEgy = Sum(vecrAvPoDeSpRot);

	/* Sort the values of the PDS to get the smallest values */
	CRealVector rSortAvPoDeSpRot(Sort(vecrAvPoDeSpRot));

	/* Average the result of smallest values and overestimate result */
	const CReal rSigmaNoise =
		Sum(rSortAvPoDeSpRot(1, NUM_SAM_IR_FOR_MIN_STAT - 1)) /
		NUM_SAM_IR_FOR_MIN_STAT * OVER_EST_FACT_MIN_STAT;

	/* Calculate signal energy by subtracting the noise energy from total
	   energy (energy cannot by negative -> bound at zero) */
	const CReal rSigEnergyBound =
		Max(rTotEgy - rSigmaNoise * iNumIntpFreqPil, (CReal) 0.0);

	/* From left to the right -> search for end of PDS */
	rEstPDSEnd = (CReal) (iNumIntpFreqPil - 1);
	rCurEnergy = (CReal) 0.0;
	bPDSResultFound = FALSE;
	for (i = 0; i < iNumIntpFreqPil; i++)
	{
		if (bPDSResultFound == FALSE)
		{
			if (rCurEnergy > rSigEnergyBound)
			{
				/* Delay index */
				rEstPDSEnd = (CReal) i;

				bPDSResultFound = TRUE;
			}

			/* Accumulate signal energy, subtract noise on each sample */
			rCurEnergy += vecrAvPoDeSpRot[i] - rSigmaNoise;
		}
	}

	/* From right to the left -> search for beginning of PDS */
	rEstPDSBegin = (CReal) 0.0;
	rCurEnergy = (CReal) 0.0;
	bPDSResultFound = FALSE;
	for (i = iNumIntpFreqPil - 1; i >= 0; i--)
	{
		if (bPDSResultFound == FALSE)
		{
			if (rCurEnergy > rSigEnergyBound)
			{
				/* Delay index */
				rEstPDSBegin = (CReal) i;

				bPDSResultFound = TRUE;
			}

			/* Accumulate signal energy, subtract noise on each sample */
			rCurEnergy += vecrAvPoDeSpRot[i] - rSigmaNoise;
		}
	}

	/* If the signal energy is too low it can happen that the estimated
	   beginning of the impulse response is before the end -> correct */
	if (rEstPDSBegin > rEstPDSEnd)
	{
		/* Set beginning and end to their maximum (minimum) value */
		rEstPDSBegin = (CReal) 0.0;
		rEstPDSEnd = (CReal) (iNumIntpFreqPil - 1);
	}

	/* Correct estimates of begin and end of PDS by the rotation */
	const CReal rPDSLenCorrection = iNumIntpFreqPil - iStPoRot + 1;
	rEstPDSBegin -= rPDSLenCorrection;
	rEstPDSEnd -= rPDSLenCorrection;

	/* Set return parameters */
	rLenPDS = rEstPDSEnd - rEstPDSBegin;
	rOffsPDS = rEstPDSBegin;
}

void CTimeSyncTrack::Init(CParameter& Parameter, int iNewSymbDelay)
{
	iNumCarrier = Parameter.iNumCarrier;
	iScatPilFreqInt = Parameter.iScatPilFreqInt;
	iNumIntpFreqPil = Parameter.iNumIntpFreqPil;
	iDFTSize = Parameter.iFFTSizeN;

	/* Timing correction history */
	iSymDelay = iNewSymbDelay;
	vecTiCorrHist.Init(iSymDelay, 0);

	/* History for new measurements (corrections) */
	veciNewMeasHist.Init(iSymDelay - 1, 0);

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

	/* Vector for averaged power delay spread estimation */
	vecrAvPoDeSp.Init(iNumIntpFreqPil, (CReal) 0.0);

	/* Lambda for IIR filter for averaging the PDS */
	rLamAvPDS = IIR1Lam(TICONST_PDS_EST_TISYNC, (CReal) SOUNDCRD_SAMPLE_RATE /
		Parameter.iSymbolBlockSize);

	/* Vector for rotated result */
	vecrAvPoDeSpRot.Init(iNumIntpFreqPil);

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

	/* Get the hamming window taps. The window is to reduce the leakage effect
	   of a DFT transformation */
	vecrHammingWindow.Init(iNumIntpFreqPil);
	vecrHammingWindow = Hamming(iNumIntpFreqPil);

	/* Weights for peak bound calculation, in Eq. (19), special values for
	   robustness mode D! */
	if (Parameter.GetWaveMode() == RM_ROBUSTNESS_MODE_D)
	{
		rConst1 = pow(10, (_REAL) -TETA1_DIST_FROM_MAX_DB_RMD / 10);
		rConst2 = pow(10, (_REAL) TETA2_DIST_FROM_MIN_DB_RMD / 10);
	}
	else
	{
		rConst1 = pow(10, (_REAL) -TETA1_DIST_FROM_MAX_DB / 10);
		rConst2 = pow(10, (_REAL) TETA2_DIST_FROM_MIN_DB / 10);
	}

	/* Define start point for rotation of detection vector for acausal taps.
	   Per definition is this point somewhere in the region after the
	   actual guard-interval window */
	if ((int) rGuardSizeFFT > iNumIntpFreqPil)
		iStPoRot = iNumIntpFreqPil;
	else
	{
		/* "+ 1" because of "Matlab indices" used in ".Merge()" function */
		iStPoRot = (int) (rGuardSizeFFT +
			Ceil((iNumIntpFreqPil - rGuardSizeFFT) / 2) + 1);
	}

	/* Init fractional part of timing correction to zero and fractional part
	   of controlling */
	rFracPartTiCor = (CReal) 0.0;
	rFracPartContr = (CReal) 0.0;

	/* Inits for the time synchronization tracking type */
	SetTiSyncTracType(TypeTiSyncTrac);

	/* Init begin and end of PDS estimation with zero and the length of guard-
	   interval respectively */
	rEstPDSBegin = (CReal) 0.0;
	rEstPDSEnd = rGuardSizeFFT;

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


	/* Inits for sample rate offset estimation ------------------------------ */
	/* Calculate number of symbols for a given time span as defined for the
	   length of the sample rate offset estimation history size */
	iLenCorrectionHist = (int) ((_REAL) SOUNDCRD_SAMPLE_RATE *
		HIST_LEN_SAM_OFF_EST_TI_CORR / Parameter.iSymbolBlockSize);

	/* Init count for acquisition */
	iResOffAcqCntMax = (int) ((_REAL) SOUNDCRD_SAMPLE_RATE *
		SAM_OFF_EST_TI_CORR_ACQ_LEN / Parameter.iSymbolBlockSize);

	/* Init sample rate offset estimation acquisition count */
	iResOffsetAcquCnt = iResOffAcqCntMax;

	veciSRTiCorrHist.Init(iLenCorrectionHist, 0); /* Init with zeros */
	iIntegTiCorrections = 0;

	/* Symbol block size converted in domain of estimated PDS */
	rSymBloSiIRDomain =
		(CReal) Parameter.iSymbolBlockSize * iNumCarrier / iDFTSize;

	/* Init variable for storing the old difference of maximum position */
	iOldNonZeroDiff = 0;
}

CReal CTimeSyncTrack::GetSamOffHz(int iDiff, int iLen)
{
	/* Calculate actual sample rate offset in Hertz */
	const CReal rCurSampOffsNorm = (CReal) iDiff / iLen / rSymBloSiIRDomain;

	return (CReal) SOUNDCRD_SAMPLE_RATE * ((CReal) 1.0 -
		(CReal) 1.0 / ((CReal) 1.0 + rCurSampOffsNorm));
}

void CTimeSyncTrack::SetTiSyncTracType(ETypeTiSyncTrac eNewTy)
{
	TypeTiSyncTrac = eNewTy;

	switch (TypeTiSyncTrac)
	{
	case TSFIRSTPEAK:
		/* Define target position for first path. Should be close to zero but
		   not exactely zero because even small estimation errors would lead to
		   ISI. The target timing position must be at least 2 samples away from
		   the guard-interval border */
		iTargetTimingPos = (int) (rGuardSizeFFT / TARGET_TI_POS_FRAC_GUARD_INT);
		if (iTargetTimingPos < 2)
			iTargetTimingPos = 2;
		break;

	case TSENERGY:
		/* No target timing position needed */
		iTargetTimingPos = 0;
		break;
	}
}

void CTimeSyncTrack::GetAvPoDeSp(CVector<_REAL>& vecrData,
								 CVector<_REAL>& vecrScale,
								_REAL& rLowerBound, _REAL& rHigherBound,
								_REAL& rStartGuard, _REAL& rEndGuard,
								_REAL& rPDSBegin, _REAL& rPDSEnd)
{
	int		i;
	int		iHalfSpec;
	_REAL	rScaleIncr;
	_REAL	rScaleAbs;

	/* Init output vectors */
	vecrData.Init(iNumIntpFreqPil, (_REAL) 0.0);
	vecrScale.Init(iNumIntpFreqPil, (_REAL) 0.0);
	rHigherBound = (_REAL) 0.0;
	rLowerBound = (_REAL) 0.0;
	rStartGuard = (_REAL) 0.0;
	rEndGuard = (_REAL) 0.0;
	rPDSBegin = (_REAL) 0.0;
	rPDSEnd = (_REAL) 0.0;

	/* Do copying of data only if vector is of non-zero length which means that
	   the module was already initialized */
	if (iNumIntpFreqPil != 0)
	{
		/* With this setting we only define the position of the guard-interval
		   in the plot. With this setting we position it centered */
		iHalfSpec = (int) ((iNumIntpFreqPil - rGuardSizeFFT) / 2);

		/* Init scale (in "ms") */
		rScaleIncr = (_REAL) iDFTSize /
			(SOUNDCRD_SAMPLE_RATE * iNumIntpFreqPil) * 1000 / iScatPilFreqInt;

		/* Let the target timing position be the "0" time */
		rScaleAbs = -(iHalfSpec + iTargetTimingPos) * rScaleIncr;

		/* Copy first part of data in output vector */
		for (i = 0; i < iHalfSpec; i++)
		{
			const _REAL rCurPDSVal =
				vecrAvPoDeSp[iNumIntpFreqPil - iHalfSpec + i];

			if (rCurPDSVal > 0)
				vecrData[i] = (_REAL) 10.0 * log10(rCurPDSVal);
			else
				vecrData[i] = RET_VAL_LOG_0;

			/* Scale */
			vecrScale[i] = rScaleAbs;
			rScaleAbs += rScaleIncr;
		}

		/* Save scale point because this is the start point of guard-interval */
		rStartGuard = rScaleAbs;

		/* Copy second part of data in output vector */
		for (i = iHalfSpec; i < iNumIntpFreqPil; i++)
		{
			const _REAL rCurPDSVal = vecrAvPoDeSp[i - iHalfSpec];

			if (rCurPDSVal > 0)
				vecrData[i] = (_REAL) 10.0 * log10(rCurPDSVal);
			else
				vecrData[i] = RET_VAL_LOG_0;

			/* Scale */
			vecrScale[i] = rScaleAbs;
			rScaleAbs += rScaleIncr;
		}

		/* Return bounds */
		switch (TypeTiSyncTrac)
		{
		case TSFIRSTPEAK:
			if (rBoundHigher > 0)
				rHigherBound = (_REAL) 10.0 * log10(rBoundHigher);
			else
				rHigherBound = RET_VAL_LOG_0;

			if (rBoundLower > 0)
				rLowerBound = (_REAL) 10.0 * log10(rBoundLower);
			else
				rLowerBound = RET_VAL_LOG_0;
			break;

		case TSENERGY:
			/* No bounds needed for energy type, set both values to "defined
			   infinity value", so it does not show up in the plot */
			rHigherBound = RET_VAL_LOG_0;
			rLowerBound = RET_VAL_LOG_0;
			break;
		}

		/* End point of guard interval */
		rEndGuard = rScaleIncr * (rGuardSizeFFT - iTargetTimingPos);

		/* Estmiated begin and end of estimated PDS */
		rPDSBegin = rScaleIncr * (rEstPDSBegin - iTargetTimingPos);
		rPDSEnd = rScaleIncr * (rEstPDSEnd - iTargetTimingPos);
	}
}