viterbidecoder.cpp

Dream.exe soft source (Visual C++)

/******************************************************************************\
 * Technische Universitaet Darmstadt, Institut fuer Nachrichtentechnik
 * Copyright (c) 2001
 *
 * Author(s):
 *	Volker Fischer, Alexander Kurpiers
 *
 * Description:
 *
 *
 ******************************************************************************
 *
 * 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 "ViterbiDecoder.h"


/* Implementation *************************************************************/
_REAL CViterbiDecoder::Decode(CVector<CDistance>& vecNewDistance,
							  CVector<_DECISION>& vecOutputBits)
{
	int				i;
	int				iDistCnt;
	int				iCurDecState;
	_VITMETRTYPE*	pCurTrelMetric;
	_VITMETRTYPE*	pOldTrelMetric;

#ifdef USE_SIMD
	/* -------------------------------------------------------------------------
	   Since the metric is 8-bit fixed-point type, we need to scale the input
	   metrics to avoid overflows */
	/* Calculate average value of input metrics */
	_REAL rAverage = (_REAL) 0.0;
	for (i = 0; i < iNumOutBitsWithMemory; i++)
	{
		rAverage += vecNewDistance[i].rTow0;
		rAverage += vecNewDistance[i].rTow1;
	}

	/* Scale input metrics */
	const _REAL rAmp = rAverage / (2 * iNumOutBitsWithMemory) / 10;
	for (i = 0; i < iNumOutBitsWithMemory; i++)
	{
		vecNewDistance[i].rTow0 /= rAmp;
		vecNewDistance[i].rTow1 /= rAmp;
	}
#endif

	/* Init pointers for old and new trellis state */
	pCurTrelMetric = vecTrelMetric1;
	pOldTrelMetric = vecTrelMetric2;

	/* Reset all metrics in the trellis. We initialize all states exept of
	   the zero state with a high metric, because we KNOW that the state "0"
	   is the transmitted state */
	pOldTrelMetric[0] = (_VITMETRTYPE) 0;
	for (i = 1; i < MC_NUM_STATES; i++)
		pOldTrelMetric[i] = MC_METRIC_INIT_VALUE;

	/* Reset counter for puncturing and distance (from metric) */
	iDistCnt = 0;


	/* Main loop over all bits ---------------------------------------------- */
	for (i = 0; i < iNumOutBitsWithMemory; i++)
	{
		/* Calculate all possible metrics ----------------------------------- */
		/* There are only a small set of possible puncturing patterns used for
		   DRM: 0001, 0101, 0011, 0111, 1111. These need different numbers of
		   input bits (increment of "iDistCnt" is dependent on pattern!). To
		   optimize the calculation of the metrics, a "subset" of bits are first
		   calculated which are used to get the final result. In this case,
		   redundancy can be avoided.
		   Note, that not all possible bit-combinations are used in the coder,
		   only a subset of numbers: 0, 2, 4, 6, 9, 11, 13, 15 (compare numbers
		   in the BUTTERFLY() calls) */

		/* Get first position in input vector (is needed for all cases) */
		const int iPos0 = iDistCnt;
		iDistCnt++;

		if (veciTablePuncPat[i] == PP_TYPE_0001)
		{
			/* Pattern 0001 */
			METRICSET(i)[ 0] = vecNewDistance[iPos0].rTow0;
			METRICSET(i)[ 2] = vecNewDistance[iPos0].rTow0;
			METRICSET(i)[ 4] = vecNewDistance[iPos0].rTow0;
			METRICSET(i)[ 6] = vecNewDistance[iPos0].rTow0;
			METRICSET(i)[ 9] = vecNewDistance[iPos0].rTow1;
			METRICSET(i)[11] = vecNewDistance[iPos0].rTow1;
			METRICSET(i)[13] = vecNewDistance[iPos0].rTow1;
			METRICSET(i)[15] = vecNewDistance[iPos0].rTow1;
		}
		else
		{
			/* The following patterns need one more bit */
			const int iPos1 = iDistCnt;
			iDistCnt++;

			/* Calculate "subsets" of bit-combinations. "rIRxx00" means that
			   the fist two bits are used, others are x-ed. "IR" stands for
			   "intermediate result" */
			const _REAL rIRxx00 =
				vecNewDistance[iPos1].rTow0 + vecNewDistance[iPos0].rTow0;
			const _REAL rIRxx10 =
				vecNewDistance[iPos1].rTow1 + vecNewDistance[iPos0].rTow0;
			const _REAL rIRxx01 =
				vecNewDistance[iPos1].rTow0 + vecNewDistance[iPos0].rTow1;
			const _REAL rIRxx11 =
				vecNewDistance[iPos1].rTow1 + vecNewDistance[iPos0].rTow1;

			if (veciTablePuncPat[i] == PP_TYPE_0101)
			{
				/* Pattern 0101 */
				METRICSET(i)[ 0] = rIRxx00;
				METRICSET(i)[ 2] = rIRxx00;
				METRICSET(i)[ 4] = rIRxx10;
				METRICSET(i)[ 6] = rIRxx10;
				METRICSET(i)[ 9] = rIRxx01;
				METRICSET(i)[11] = rIRxx01;
				METRICSET(i)[13] = rIRxx11;
				METRICSET(i)[15] = rIRxx11;
			}
			else if (veciTablePuncPat[i] == PP_TYPE_0011)
			{
				/* Pattern 0011 */
				METRICSET(i)[ 0] = rIRxx00;
				METRICSET(i)[ 2] = rIRxx10;
				METRICSET(i)[ 4] = rIRxx00;
				METRICSET(i)[ 6] = rIRxx10;
				METRICSET(i)[ 9] = rIRxx01;
				METRICSET(i)[11] = rIRxx11;
				METRICSET(i)[13] = rIRxx01;
				METRICSET(i)[15] = rIRxx11;
			}
			else
			{
				/* The following patterns need one more bit */
				const int iPos2 = iDistCnt;
				iDistCnt++;

				if (veciTablePuncPat[i] == PP_TYPE_0111)
				{
					/* Pattern 0111 */
					METRICSET(i)[ 0] = vecNewDistance[iPos2].rTow0 + rIRxx00;
					METRICSET(i)[ 2] = vecNewDistance[iPos2].rTow0 + rIRxx10;
					METRICSET(i)[ 4] = vecNewDistance[iPos2].rTow1 + rIRxx00;
					METRICSET(i)[ 6] = vecNewDistance[iPos2].rTow1 + rIRxx10;
					METRICSET(i)[ 9] = vecNewDistance[iPos2].rTow0 + rIRxx01;
					METRICSET(i)[11] = vecNewDistance[iPos2].rTow0 + rIRxx11;
					METRICSET(i)[13] = vecNewDistance[iPos2].rTow1 + rIRxx01;
					METRICSET(i)[15] = vecNewDistance[iPos2].rTow1 + rIRxx11;
				}
				else
				{
					/* Pattern 1111 */
					/* This pattern needs all four bits */
					const int iPos3 = iDistCnt;
					iDistCnt++;

					/* Calculate "subsets" of bit-combinations. "rIRxx00" means
					   that the last two bits are used, others are x-ed.
					   "IR" stands for "intermediate result" */
					const _REAL rIR00xx = vecNewDistance[iPos3].rTow0 +
						vecNewDistance[iPos2].rTow0;
					const _REAL rIR10xx = vecNewDistance[iPos3].rTow1 +
						vecNewDistance[iPos2].rTow0;
					const _REAL rIR01xx = vecNewDistance[iPos3].rTow0 +
						vecNewDistance[iPos2].rTow1;
					const _REAL rIR11xx = vecNewDistance[iPos3].rTow1 +
						vecNewDistance[iPos2].rTow1;

					METRICSET(i)[ 0] = rIR00xx + rIRxx00; /* 0 */
					METRICSET(i)[ 2] = rIR00xx + rIRxx10; /* 2 */
					METRICSET(i)[ 4] = rIR01xx + rIRxx00; /* 4 */
					METRICSET(i)[ 6] = rIR01xx + rIRxx10; /* 6 */
					METRICSET(i)[ 9] = rIR10xx + rIRxx01; /* 9 */
					METRICSET(i)[11] = rIR10xx + rIRxx11; /* 11 */
					METRICSET(i)[13] = rIR11xx + rIRxx01; /* 13 */
					METRICSET(i)[15] = rIR11xx + rIRxx11; /* 15 */
				}
			}
		}


		/* Update trellis --------------------------------------------------- */
#ifdef USE_SIMD
		/* Use the butterfly unroll for reordering the metrics for SIMD
		   trellis */
#define BUTTERFLY(cur, next, prev0, prev1, met0, met1) \
		{ \
			/* At this point we convert from float to char! No overflow-check
			   is done here */ \
			chMet1[prev0] = (_VITMETRTYPE) METRICSET(i)[met0]; \
			chMet2[prev0] = (_VITMETRTYPE) METRICSET(i)[met1]; \
		}
#else
		/* c++ version of trellis update */
#define BUTTERFLY(cur, next, prev0, prev1, met0, met1) \
		{ \
			/* First state in this set ------------------------------------ */ \
			/* Calculate metrics from the two previous states, use the old
			   metric from the previous states plus the "transition-metric" */ \
			const _VITMETRTYPE rFiStAccMetricPrev0 = \
				pOldTrelMetric[prev0] + METRICSET(i)[met0]; \
			const _VITMETRTYPE  rFiStAccMetricPrev1 = \
				pOldTrelMetric[prev1] + METRICSET(i)[met1]; \
			\
			/* Take path with smallest metric */ \
			if (rFiStAccMetricPrev0 < rFiStAccMetricPrev1) \
			{ \
				/* Save minimum metric for this state and store decision */ \
				pCurTrelMetric[cur] = rFiStAccMetricPrev0; \
				matdecDecisions[i][cur] = 0; \
			} \
			else \
			{ \
				/* Save minimum metric for this state and store decision */ \
				pCurTrelMetric[cur] = rFiStAccMetricPrev1; \
				matdecDecisions[i][cur] = 1; \
			} \
			\
			/* Second state in this set ----------------------------------- */ \
			/* The only difference is that we swapped the matric sets */ \
			const _VITMETRTYPE rSecStAccMetricPrev0 = \
				pOldTrelMetric[prev0] + METRICSET(i)[met1]; \
			const _VITMETRTYPE  rSecStAccMetricPrev1 = \
				pOldTrelMetric[prev1] + METRICSET(i)[met0]; \
			\
			/* Take path with smallest metric */ \
			if (rSecStAccMetricPrev0 < rSecStAccMetricPrev1) \
			{ \
				/* Save minimum metric for this state and store decision */ \
				pCurTrelMetric[next] = rSecStAccMetricPrev0; \
				matdecDecisions[i][next] = 0; \
			} \
			else \
			{ \
				/* Save minimum metric for this state and store decision */ \
				pCurTrelMetric[next] = rSecStAccMetricPrev1; \
				matdecDecisions[i][next] = 1; \
			} \
		}
#endif

		/* Unroll butterflys to avoid loop overhead. For c++ version, the
		   actual calculation of the trellis update is done here, for MMX
		   version, only the reordering of the new metrics is done here */
		BUTTERFLY( 0,  1,  0, 32,  0, 15)
		BUTTERFLY( 2,  3,  1, 33,  6,  9)
		BUTTERFLY( 4,  5,  2, 34, 11,  4)
		BUTTERFLY( 6,  7,  3, 35, 13,  2)
		BUTTERFLY( 8,  9,  4, 36, 11,  4)
		BUTTERFLY(10, 11,  5, 37, 13,  2)
		BUTTERFLY(12, 13,  6, 38,  0, 15)
		BUTTERFLY(14, 15,  7, 39,  6,  9)
		BUTTERFLY(16, 17,  8, 40,  4, 11)
		BUTTERFLY(18, 19,  9, 41,  2, 13)
		BUTTERFLY(20, 21, 10, 42, 15,  0)
		BUTTERFLY(22, 23, 11, 43,  9,  6)
		BUTTERFLY(24, 25, 12, 44, 15,  0)
		BUTTERFLY(26, 27, 13, 45,  9,  6)
		BUTTERFLY(28, 29, 14, 46,  4, 11)
		BUTTERFLY(30, 31, 15, 47,  2, 13)
		BUTTERFLY(32, 33, 16, 48,  9,  6)
		BUTTERFLY(34, 35, 17, 49, 15,  0)
		BUTTERFLY(36, 37, 18, 50,  2, 13)
		BUTTERFLY(38, 39, 19, 51,  4, 11)
		BUTTERFLY(40, 41, 20, 52,  2, 13)
		BUTTERFLY(42, 43, 21, 53,  4, 11)
		BUTTERFLY(44, 45, 22, 54,  9,  6)
		BUTTERFLY(46, 47, 23, 55, 15,  0)
		BUTTERFLY(48, 49, 24, 56, 13,  2)
		BUTTERFLY(50, 51, 25, 57, 11,  4)
		BUTTERFLY(52, 53, 26, 58,  6,  9)
		BUTTERFLY(54, 55, 27, 59,  0, 15)
		BUTTERFLY(56, 57, 28, 60,  6,  9)
		BUTTERFLY(58, 59, 29, 61,  0, 15)
		BUTTERFLY(60, 61, 30, 62, 13,  2)
		BUTTERFLY(62, 63, 31, 63, 11,  4)


#undef BUTTERFLY

#ifdef USE_SIMD
		/* Do actual trellis update in separate file (assembler implementation) */
#ifdef USE_MMX
		TrellisUpdateMMX(
#endif
#ifdef USE_SSE2
		TrellisUpdateSSE2(
#endif
			&matdecDecisions[i][0], pCurTrelMetric, pOldTrelMetric,
			chMet1, chMet2);
#endif

#ifdef USE_MAX_LOG_MAP
		/* Store accumulated metrics for backward Viterbi */
		for (int j = 0; j < MC_NUM_STATES; j++)
			matrAlpha[i][j] = pOldTrelMetric[j];
#endif

		/* Swap trellis data pointers (old -> new, new -> old) */
		_VITMETRTYPE* pTMPTrelMetric = pCurTrelMetric;
		pCurTrelMetric = pOldTrelMetric;
		pOldTrelMetric = pTMPTrelMetric;
	}


#ifdef USE_MAX_LOG_MAP
	/* MAX-LOG MAP implementation ------------------------------------------- */