freqsyncacq.cpp

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

#ifdef _DEBUG_
/* Stores all important parameters for last shot */
FILE* pFile1 = fopen("test/freqacq.dat", "w");
int iPeakCnt = 0;
for (i = 1; i < iSearchWinSize; i++) {
	_REAL rPM, rFinPM; /* rPM: peak marker, rFinPM: final peak marker */
	if (iPeakCnt < iNumDetPeaks) {
		if (i == veciPeakIndex[iPeakCnt]) {
			rPM = vecrPSDPilCor[i];
			if (vecbFlagVec[iPeakCnt] == 1) rFinPM = vecrPSDPilCor[i]; else rFinPM = 0;
			iPeakCnt++;
		} else {rPM = 0; rFinPM = 0;}
	} else {rPM = 0; rFinPM = 0;}
	fprintf(pFile1, "%e %e %e\n", vecrPSDPilCor[i], rPM, rFinPM);
}
fclose(pFile1);
// close all;load freqacq.dat;load freqacqFilt.dat;subplot(211),semilogy(freqacq(:,1));hold;plot(freqacq(:,2),'*y');plot(freqacq(:,3),'*k');subplot(212),semilogy(freqacqFilt)
// close all;load freqacq.dat;semilogy(freqacq(:,1));hold;plot(freqacq(:,2),'*y');plot(freqacq(:,3),'*k');
#endif


						/* -----------------------------------------------------
						   An acquisition frequency offest estimation was
						   found */
						/* Calculate frequency offset and set global parameter
						   for offset */
						ReceiverParam.rFreqOffsetAcqui =
							(_REAL) iMaxIndex / iFrAcFFTSize;

						/* Reset acquisition flag */
						bAquisition = FALSE;


						/* Send out the data stored for FFT calculation ----- */
						/* This does not work for bandpass filter. TODO: make
						   this possible for bandpass filter, too */
						if (bUseRecFilter == FALSE)
						{
							iOutputBlockSize = iHistBufSize;

							/* Frequency offset correction */
							const _REAL rNormCurFreqOffsFst = (_REAL) 2.0 * crPi *
								(ReceiverParam.rFreqOffsetAcqui - rInternIFNorm);

							for (i = 0; i < iHistBufSize; i++)
							{
								/* Multiply with exp(j omega t) */
								(*pvecOutputData)[i] = vecrFFTHistory[i] *
									_COMPLEX(Cos(i * rNormCurFreqOffsFst),
									Sin(-i * rNormCurFreqOffsFst));
							}

							/* Init "exp-step" for regular frequency shift which
							   is used in tracking mode to get contiuous mixing
							   signal */
							cCurExp =
								_COMPLEX(Cos(iHistBufSize * rNormCurFreqOffsFst),
								Sin(-iHistBufSize * rNormCurFreqOffsFst));
						}
					}
				}
			}
		}
	}
	else
	{
		/* If synchronized DRM input stream is used, overwrite the detected
		   frequency offest estimate by the desired frequency, because we know
		   this value */
		if (bSyncInput == TRUE)
		{
			ReceiverParam.rFreqOffsetAcqui =
				(_REAL) ReceiverParam.iIndexDCFreq / ReceiverParam.iFFTSizeN;
		}

		/* Use the same block size as input block size */
		iOutputBlockSize = iInputBlockSize;


		/* Frequency offset correction -------------------------------------- */
		/* Total frequency offset from acquisition and tracking (we calculate
		   the normalized frequency offset) */
		const _REAL rNormCurFreqOffset =
			(_REAL) 2.0 * crPi * (ReceiverParam.rFreqOffsetAcqui +
			ReceiverParam.rFreqOffsetTrack - rInternIFNorm);

		/* New rotation vector for exp() calculation */
		const _COMPLEX cExpStep =
			_COMPLEX(Cos(rNormCurFreqOffset), Sin(rNormCurFreqOffset));

		/* Input data is real, make complex and compensate for frequency
		   offset */
		for (i = 0; i < iOutputBlockSize; i++)
		{
			(*pvecOutputData)[i] = (*pvecInputData)[i] * Conj(cCurExp);

			/* Rotate exp-pointer on step further by complex multiplication with
			   precalculated rotation vector cExpStep. This saves us from
			   calling sin() and cos() functions all the time (iterative
			   calculation of these functions) */
			cCurExp *= cExpStep;
		}


		/* Bandpass filter -------------------------------------------------- */
		if (bUseRecFilter == TRUE)
			BPFilter.Process(*pvecOutputData);
	}
}

void CFreqSyncAcq::InitInternal(CParameter& ReceiverParam)
{
	/* Needed for calculating offset in Hertz in case of synchronized input
	   (for simulation) */
	iFFTSize = ReceiverParam.iFFTSizeN;

	/* We using parameters from robustness mode B as pattern for the desired
	   frequency pilot positions */
	veciTableFreqPilots[0] =
		iTableFreqPilRobModB[0][0] * NUM_BLOCKS_4_FREQ_ACQU;
	veciTableFreqPilots[1] =
		iTableFreqPilRobModB[1][0] * NUM_BLOCKS_4_FREQ_ACQU;
	veciTableFreqPilots[2] =
		iTableFreqPilRobModB[2][0] * NUM_BLOCKS_4_FREQ_ACQU;

	/* Size of FFT */
	iFrAcFFTSize = RMB_FFT_SIZE_N * NUM_BLOCKS_4_FREQ_ACQU;


	/* -------------------------------------------------------------------------
	   Set start- and endpoint of search window for DC carrier after the
	   correlation with the known pilot structure */
	/* Normalize the desired position and window size which are in Hertz */
	const _REAL rNormDesPos = rCenterFreq / SOUNDCRD_SAMPLE_RATE * 2;
	const _REAL rNormHalfWinSize = rWinSize / SOUNDCRD_SAMPLE_RATE;

	/* Length of the half of the spectrum of real input signal (the other half
	   is the same because of the real input signal). We have to consider the
	   Nyquist frequency ("iFrAcFFTSize" is always even!) */
	iHalfBuffer = iFrAcFFTSize / 2 + 1;

	/* Search window is smaller than haft-buffer size because of correlation
	   with pilot positions */
	iSearchWinSize = iHalfBuffer - veciTableFreqPilots[2];

	/* Calculate actual indices of start and end of search window */
	iStartDCSearch =
		(int) Floor((rNormDesPos - rNormHalfWinSize) * iHalfBuffer);
	iEndDCSearch = (int) Ceil((rNormDesPos + rNormHalfWinSize) * iHalfBuffer);

	/* Check range. If out of range -> correct */
	if (!((iStartDCSearch > 0) && (iStartDCSearch < iSearchWinSize)))
		iStartDCSearch = 0;

	if (!((iEndDCSearch > 0) && (iEndDCSearch < iSearchWinSize)))
		iEndDCSearch = iSearchWinSize;

	/* Set bound for ratio between filtered signal to signal. Use a lower bound
	   if the search window is smaller */
	if (((_REAL) iEndDCSearch - iStartDCSearch) / iHalfBuffer < (_REAL) 0.042)
		rPeakBoundFiltToSig = PEAK_BOUND_FILT2SIGNAL_0_042;
	else
		rPeakBoundFiltToSig = PEAK_BOUND_FILT2SIGNAL_1;


	/* Init vectors and FFT-plan -------------------------------------------- */
	/* Allocate memory for FFT-histories and init with zeros */
	iHistBufSize = iFrAcFFTSize * NUM_BLOCKS_USED_FOR_AV;
	vecrFFTHistory.Init(iHistBufSize, (_REAL) 0.0);
	vecrFFTInput.Init(iFrAcFFTSize);
	vecrSqMagFFTOut.Init(iHalfBuffer);

	/* Allocate memory for PSD after pilot correlation */
	vecrPSDPilCor.Init(iHalfBuffer);

	/* Init vectors for filtering in frequency direction */
	vecrFiltResLR.Init(iHalfBuffer);
	vecrFiltResRL.Init(iHalfBuffer);
	vecrFiltRes.Init(iHalfBuffer);

	/* Index memory for detected peaks (assume worst case with the size) */
	veciPeakIndex.Init(iHalfBuffer);

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

	/* Init Hamming window */
	vecrHammingWin.Init(iFrAcFFTSize);
	vecrHammingWin = Hamming(iFrAcFFTSize);

	/* Init moving average class for SqMag FFT results */
	vvrPSDMovAv.InitVec(NUM_FFT_RES_AV_BLOCKS, iHalfBuffer);


	/* Frequency correction */
	/* Start with phase null (arbitrary) */
	cCurExp = (_REAL) 1.0;
	rInternIFNorm = (_REAL) ReceiverParam.iIndexDCFreq / iFFTSize;


	/* Init bandpass filter object */
	BPFilter.Init(ReceiverParam.iSymbolBlockSize, VIRTUAL_INTERMED_FREQ,
		ReceiverParam.GetSpectrumOccup(), CDRMBandpassFilt::FT_RECEIVER);


	/* Define block-sizes for input (The output block size is set inside
	   the processing routine, therefore only a maximum block size is set
	   here) */
	iInputBlockSize = ReceiverParam.iSymbolBlockSize;

	/* We have to consider that the next module can take up to two symbols per
	   step. This can be satisfied be multiplying with "3". We also want to ship
	   the whole FFT buffer after finishing the frequency acquisition so that
	   these samples can be reused for synchronization and do not get lost */
	iMaxOutputBlockSize = 3 * ReceiverParam.iSymbolBlockSize + iHistBufSize;
}

void CFreqSyncAcq::SetSearchWindow(_REAL rNewCenterFreq, _REAL rNewWinSize)
{
	/* Set internal parameters */
	rCenterFreq = rNewCenterFreq;
	rWinSize = rNewWinSize;

	/* Set flag to initialize the module to the new parameters */
	SetInitFlag();
}

void CFreqSyncAcq::StartAcquisition()
{
	/* Set flag so that the actual acquisition routine is entered */
	bAquisition = TRUE;

	/* Reset (or init) counters */
	iAquisitionCounter = NUM_BLOCKS_4_FREQ_ACQU;
	iAverageCounter = NUM_FFT_RES_AV_BLOCKS;

	/* Reset FFT-history */
	vecrFFTHistory.Reset((_REAL) 0.0);
}