timesync.cpp

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

/******************************************************************************\
 * Technische Universitaet Darmstadt, Institut fuer Nachrichtentechnik
 * Copyright (c) 2001
 *
 * Author(s):
 *	Volker Fischer
 *
 * Description:
 *	Time synchronization
 * This module can have different amounts of input data. If two
 * possible FFT-window positions are found, the next time no new block is
 * requested.
 *
 * Robustness-mode detection
 *
 ******************************************************************************
 *
 * 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 "TimeSync.h"


/* Implementation *************************************************************/
void CTimeSync::ProcessDataInternal(CParameter& ReceiverParam)
{
	int				i, j, k;
	int				iMaxIndex;
	int				iNewStIndCount; 
	int				iIntDiffToCenter;
	int				iCurPos;
	int				iDetectedRModeInd;
	int				iDecInpuSize;
	CReal			rMaxValue;
	CReal			rMaxValRMCorr;
	CReal			rSecHighPeak;
	CReal			rFreqOffsetEst;
	CComplexVector	cvecInpTmp;
	CRealVector		rResMode(NUM_ROBUSTNESS_MODES);

	/* Max number of detected peaks ("5" for safety reasons. Could be "2") */
	CVector<int>	iNewStartIndexField(5);

	/* Init start index (in case no timing could be detected or init phase) */
	int iStartIndex = iSymbolBlockSize;

	/* Write new block of data at the end of shift register */
	HistoryBuf.AddEnd(*pvecInputData, iInputBlockSize);

	/* In case the time domain frequency offset estimation method is activated,
	   the hilbert filtering of input signal must always be applied */
#ifndef USE_FRQOFFS_TRACK_GUARDCORR
	if ((bTimingAcqu == TRUE) || (bRobModAcqu == TRUE))
#endif
	{
		/* ---------------------------------------------------------------------
		   Data must be band-pass-filtered before applying the algorithms,
		   because we do not know which mode is active when we synchronize
		   the timing and we must assume the worst-case, therefore use only
		   from DC to 4.5 kHz. If the macro "USE_10_KHZ_HILBFILT" is defined,
		   a bandwith of approx. 10 kHz is used. In this case, we have better
		   performance with the most probable 10 kHz mode but may have worse
		   performance with the 4.5 or 5 kHz modes (for the acquisition) */

		/* The FIR filter intermediate buffer must be adjusted to the new
		   input block size since the size can be vary with time */
		cvecInpTmp.Init(iInputBlockSize);

		/* Copy CVector data in CMatlibVector */
		for (i = 0; i < iInputBlockSize; i++)
			cvecInpTmp[i] = (*pvecInputData)[i];

		/* Complex Hilbert filter. We use the copy constructor for storing
		   the result since the sizes of the output vector varies with time.
		   We decimate the signal with this function, too, because we only
		   analyze a spectrum bandwith of approx. 5 [10] kHz */
		CComplexVector cvecOutTmp(
			FirFiltDec(cvecB, cvecInpTmp, cvecZ, GRDCRR_DEC_FACT));

		/* Get size of new output vector */
		iDecInpuSize = Size(cvecOutTmp);

		/* Copy data from Matlib vector in regular vector for storing in
		   shift register
		   TODO: Make vector types compatible (or maybe only use matlib vectors
		   everywhere) */
		cvecOutTmpInterm.Init(iDecInpuSize);
		for (i = 0; i < iDecInpuSize; i++)
			cvecOutTmpInterm[i] = cvecOutTmp[i];

		/* Write new block of data at the end of shift register */
		HistoryBufCorr.AddEnd(cvecOutTmpInterm, iDecInpuSize);


#ifdef USE_FRQOFFS_TRACK_GUARDCORR
		/* ---------------------------------------------------------------------
		   Frequency offset tracking estimation method based on guard-interval
		   correlation.
		   Problem: The tracking frequency range is "GRDCRR_DEC_FACT"-times
		   smaller than the one in the frequency domain tracking algorithm.
		   Therefore, the acquisition unit must work more precisely
		   (see FreqSyncAcq.h) */

		/* Guard-interval correlation at ML estimated timing position */
		/* Calculate start points for correlation. Consider delay from
		   Hilbert-filter */
		const int iHalHilFilDelDec = NUM_TAPS_HILB_FILT / 2 / GRDCRR_DEC_FACT;
		const int iCorrPosFirst = iDecSymBS + iHalHilFilDelDec;
		const int iCorrPosSec =
			iDecSymBS + iHalHilFilDelDec + iLenUsefPart[iSelectedMode];

		/* Actual correlation over entire guard-interval */
		CComplex cGuardCorrFreqTrack = 0.0;
		for (i = 0; i < iLenGuardInt[iSelectedMode]; i++)
		{
			/* Use start point from ML timing estimation. The input stream is
			   automatically adjusted to have this point at "iDecSymBS" */
			cGuardCorrFreqTrack += HistoryBufCorr[i + iCorrPosFirst] *
				Conj(HistoryBufCorr[i + iCorrPosSec]);
		}

		/* Average vector, real and imaginary part separately */
		IIR1(cFreqOffAv, cGuardCorrFreqTrack, rLamFreqOff);

		/* Calculate argument */
		const CReal rFreqOffsetEst = Angle(cFreqOffAv);

		/* Correct measurement average for actually applied frequency
		   correction */
		cFreqOffAv *= CComplex(Cos(-rFreqOffsetEst), Sin(-rFreqOffsetEst));

		/* Integrate the result for controling the frequency offset, normalize
		   estimate */
		ReceiverParam.rFreqOffsetTrack -= rFreqOffsetEst * rNormConstFOE;
#endif


		if ((bTimingAcqu == TRUE) || (bRobModAcqu == TRUE))
		{
			/* Guard-interval correlation ----------------------------------- */
			/* Set position pointer back for this block */
			iTimeSyncPos -= iDecInpuSize;

			/* Init start-index count */
			iNewStIndCount = 0;

			/* We use the block in the middle of the buffer for observation */
			for (i = iDecSymBS + iDecSymBS - iDecInpuSize;
				i < iDecSymBS + iDecSymBS; i++)
			{
				/* Only every "iStepSizeGuardCorr"'th value is calculated for
				   efficiency reasons */
				if (i == iTimeSyncPos)
				{
					/* Do the following guard interval correlation for all
					   possible robustness modes (this is needed for robustness
					   mode detection) */
					for (j = 0; j < NUM_ROBUSTNESS_MODES; j++)
					{
						/* Guard-interval correlation ----------------------- */
						/* Speed optimized calculation of the guard-interval
						   correlation. We devide the total block, which has to
						   be computed, in parts of length "iStepSizeGuardCorr".
						   The results of these blocks are stored in a vector.
						   Now, only one new part has to be calculated and one
						   old one has to be subtracted from the global result.
						   Special care has to be taken since "iGuardSize" must
						   not be a multiple of "iStepSizeGuardCorr". Therefore
						   the "if"-condition */
						/* First subtract correlation values shifted out */
						cGuardCorr[j] -= 
							veccIntermCorrRes[j][iPosInIntermCResBuf[j]];
						rGuardPow[j] -= 
							vecrIntermPowRes[j][iPosInIntermCResBuf[j]];

						/* Calculate new block and add in memory */
						for (k = iLengthOverlap[j]; k < iLenGuardInt[j]; k++)
						{
							/* Actual correlation */
							iCurPos = iTimeSyncPos + k;
							cGuardCorrBlock[j] += HistoryBufCorr[iCurPos] * 
								Conj(HistoryBufCorr[iCurPos + iLenUsefPart[j]]);

							/* Energy calculation for ML solution */
							rGuardPowBlock[j] +=
								SqMag(HistoryBufCorr[iCurPos]) +
								SqMag(HistoryBufCorr[iCurPos + iLenUsefPart[j]]);

							/* If one complete block is ready -> store it. We
							   need to add "1" to the k, because otherwise
							   "iLengthOverlap" would satisfy the
							   "if"-condition */
							if (((k + 1) % iStepSizeGuardCorr) == 0)
							{
								veccIntermCorrRes[j][iPosInIntermCResBuf[j]] =
									cGuardCorrBlock[j];

								vecrIntermPowRes[j][iPosInIntermCResBuf[j]] =
									rGuardPowBlock[j];

								/* Add the new block to the global result */
								cGuardCorr[j] += cGuardCorrBlock[j];
								rGuardPow[j] += rGuardPowBlock[j];

								/* Reset block result */
								cGuardCorrBlock[j] = (CReal) 0.0;
								rGuardPowBlock[j] = (CReal) 0.0;

								/* Increase position pointer and test if wrap */
								iPosInIntermCResBuf[j]++;
								if (iPosInIntermCResBuf[j] == iLengthIntermCRes[j])
									iPosInIntermCResBuf[j] = 0;
							}
						}

						/* Save correlation results in shift register */
						for (k = 0; k < iRMCorrBufSize - 1; k++)
							vecrRMCorrBuffer[j][k] = vecrRMCorrBuffer[j][k + 1];

						/* ML solution */
						vecrRMCorrBuffer[j][iRMCorrBufSize - 1] =
							abs(cGuardCorr[j] + cGuardCorrBlock[j]) - 
							(rGuardPow[j] + rGuardPowBlock[j]) / 2;
					}

					/* Energy of guard intervall calculation and detection of
					   peak is only needed if timing aquisition is true */
					if (bTimingAcqu == TRUE)
					{
						/* Start timing detection not until initialization phase
						   is finished */
						if (iTiSyncInitCnt > 1)
						{
							/* Decrease counter */
							iTiSyncInitCnt--;
						}
						else
						{
							/* Average the correlation results */
							IIR1(vecCorrAvBuf[iCorrAvInd],
								vecrRMCorrBuffer[iSelectedMode][iRMCorrBufSize - 1],
								1 - rLambdaCoAv);


							/* Energy of guard-interval correlation calculation
							   (this is simply a moving average operation) */
							vecrGuardEnMovAv.Add(vecCorrAvBuf[iCorrAvInd]);


							/* Taking care of correlation average buffer ---- */
							/* We use a "cyclic buffer" structure. This index
							   defines the position in the buffer */
							iCorrAvInd++;
							if (iCorrAvInd == iMaxDetBufSize)
							{
								/* Adaptation of the lambda parameter for
								   guard-interval correlation averaging IIR
								   filter. With this adaptation we achieve
								   better averaging results. A lower bound is
								   defined for this parameter */
								if (rLambdaCoAv <= 0.1)
									rLambdaCoAv = 0.1;
								else
									rLambdaCoAv /= 2;

								iCorrAvInd = 0;
							}


							/* Detection buffer ----------------------------- */
							/* Update buffer for storing the moving average
							   results */
							pMaxDetBuffer.AddEnd(vecrGuardEnMovAv.GetAverage());

							/* Search for maximum */
							iMaxIndex = 0;
							rMaxValue = (CReal) -_MAXREAL; /* Init value */
							for (k = 0; k < iMaxDetBufSize; k++)
							{
								if (pMaxDetBuffer[k] > rMaxValue)
								{
									rMaxValue = pMaxDetBuffer[k];
									iMaxIndex = k;
								}
							}

							/* If maximum is in the middle of the interval, mark
							   position as the beginning of the FFT window */
							if (iMaxIndex == iCenterOfMaxDetBuf)
							{
								/* The optimal start position for the FFT-window
								   is the middle of the "MaxDetBuffer" */
								iNewStartIndexField[iNewStIndCount] = 
									iTimeSyncPos * GRDCRR_DEC_FACT -
									iSymbolBlockSize / 2 -
									/* Compensate for Hilbert-filter delay. The
									   delay is introduced in the downsampled
									   domain, therefore devide it by
									   "GRDCRR_DEC_FACT" */
									NUM_TAPS_HILB_FILT / 2 / GRDCRR_DEC_FACT;

								iNewStIndCount++;
							}
						}
					}

					/* Set position pointer to next step */
					iTimeSyncPos += iStepSizeGuardCorr;
				}
			}


			/* Robustness mode detection ------------------------------------ */
			if (bRobModAcqu == TRUE)
			{
				/* Start robustness mode detection not until the buffer is
				   filled */
				if (iRobModInitCnt > 1)
				{
					/* Decrease counter */
					iRobModInitCnt--;
				}
				else
				{
					/* Correlation of guard-interval correlation with prepared
					   cos-vector. Store highest peak */
					rMaxValRMCorr = (CReal) 0.0;
					for (j = 0; j < NUM_ROBUSTNESS_MODES; j++)
					{
						/* Correlation with symbol rate frequency (Correlations
						   must be normalized to be comparable!
						   ("/ iGuardSizeX")) */
						rResMode[j] =
							Abs(Sum(vecrRMCorrBuffer[j] * vecrCos[j])) /
							iLenGuardInt[j];

						/* Search for maximum */
						if (rResMode[j] > rMaxValRMCorr)
						{
							rMaxValRMCorr = rResMode[j];
							iDetectedRModeInd = j;
						}
					}

					/* Get second highest peak */
					rSecHighPeak = (CReal) 0.0;
					for (j = 0; j < NUM_ROBUSTNESS_MODES; j++)
					{
						if ((rResMode[j] > rSecHighPeak) &&
							(iDetectedRModeInd != j))
						{
							rSecHighPeak = rResMode[j];
						}
					}

					/* Find out if we have a reliable measure
					   (distance to next peak) */
					if ((rMaxValRMCorr / rSecHighPeak) > THRESHOLD_RELI_MEASURE)
					{
						/* Reset aquisition flag for robustness mode detection */
						bRobModAcqu = FALSE;

						/* Set wave mode */
						if (ReceiverParam.
							SetWaveMode(GetRModeFromInd(iDetectedRModeInd)) == TRUE)
						{
							/* Reset output cyclic-buffer because wave mode has
							   changed and the data written in the buffer is not
							   valid anymore */
							SetBufReset1();
						}
					}
				}
			}
		}
	}

	if (bTimingAcqu == TRUE)
	{
		/* Use all measured FFT-window start points for determining the "real"
		   one */	
		for (i = 0; i < iNewStIndCount; i++)
		{
			/* Check if new measurement is in range of predefined bound. This
			   bound shall eliminate outliers for the calculation of the
			   filtered result */
			if (((iNewStartIndexField[i] < (iCenterOfBuf + TIMING_BOUND_ABS)) &&
				(iNewStartIndexField[i] > (iCenterOfBuf - TIMING_BOUND_ABS))))
			{
				/* New measurement is in range -> use it for filtering */
				/* Low-pass filter detected start of frame */
				IIR1(rStartIndex, (CReal) iNewStartIndexField[i],
					LAMBDA_LOW_PASS_START);

				/* Reset counters for non-linear correction algorithm */
				iCorrCounter = 0;
				iAveCorr = 0;

				/* GUI message that timing is ok */
				PostWinMessage(MS_TIME_SYNC, 0);
				ReceiverParam.ReceptLog.SetSync(TRUE);

				/* Acquisition was successful, reset init flag (just in case it
				   was not reset by the non-linear correction unit */
				bInitTimingAcqu = FALSE;
			}
			else
			{
				/* Non-linear correction of the filter-output to ged rid of
				   large differences between current measurement and
				   filter-output */
				iCorrCounter++;

				/* Average the NUM_SYM_BEFORE_RESET measurements for reset
				   rStartIndex */
				iAveCorr += iNewStartIndexField[i];