log_map_decode_c.c

用Matlab仿真turbocode的算法

/******************************************************************************/
/*                                                                            */
/* log_map_decode_c.c                                                         */
/*                                                                            */
/******************************************************************************/
/* Author(s)    : Andreas Dittrich                                            */
/* Project start: 01.02.2001                                                  */
/* Completed    : 05.02.2001                                                  */
/* Last change  : 04.04.2001                                                  */
/* Time         : 10:15                                                       */
/******************************************************************************/
/* Changes      :                                                             */
/******************************************************************************/



/******************************************************************************/
/* [L_ud, L_cd] = log_MAP_decode_c( L_u, L_c, L_Anfangszustand, L_Endzustand, */
/*                                  State_Table )                             */
/*                                                                            */
/*    L_ud = Soft-out decoded Symbols  qu x N - Matrix                        */
/*    L_cd = Reliability coded (input-)Symbols  qc x N*itp - Matrix           */
/*                                                                            */
/*    L_u = A-Priori Log-Prob. about decoded Symbols  qu x N - Matrix         */
/*    L_c = Soft-in coded Symbols (Log-Prob)   qc x N*itp - Matrix            */
/*                                                                            */
/*    L_Anfangsz.: log-Wahrscheinlichkeiten f黵 den Anfangszustand   qu^L x 1 */
/*    L_Endz.: log-Wahrscheinlichkeiten f黵 den Endzustand       qu^L x 1     */
/*                                                                            */
/*    State_Table: (itp+3) x qu^K - Matrix                                    */
/*    State_Table(1,:) : zugeh鰎ige Eingangssymbole pro Moore-Zustand  1..qu  */
/*    State_Table(2:(itp+1),:) : Ausgangssymbole pro Moore-Zustand   1..qc    */
/*    State_Table(itp+2,:) : linker Knoten (Mealy-Zustand)  1..qu^L           */
/*    State_Table(itp+3,:) : rechter Knoten (Mealy-Zustand)   1..qu^L         */
/*                                                                            */
/*    qu = M鋍htigkeit des Eingangs-Symbolalphabets                           */
/*    qc = M鋍htigkeit des Ausgangs-Symbolalphabets                           */
/******************************************************************************/


#include <math.h>
#include <stdlib.h>
#include "mex.h"
#include "matrix.h"

#ifdef LINKED

#include "map_decoder.h"

#endif


#define INF 1.0E10
#define Nargs_rhs_str "5"
#define Nargs_rhs 5
#define Nargs_lhs_str "2"
#define Nargs_lhs 2
#define PROGNAME "log_MAP_decode_c"


void MAP_c( double L_ud[], double L_cd[],
            double L_u[], double L_c[], double L_Anfangszustand[], double L_Endzustand[],
            int ST[], const int qu, const int qc,
            const int qhL, const int qhK, const int itp, const int N );


/******************************************************************************/
/*  Gateway-Routine                                                           */
/******************************************************************************/

#ifndef LINKED

void mexFunction(
                 int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]
		 )
{
	double *L_ud, *L_cd;
	double *L_u, *L_c, *L_Anfangszustand, *L_Endzustand, *State_Table;
	int *State_Table_int;
	mxArray *State_Table_int_array_ptr;
  	int State_Table_int_dims[2];

	int qu, qc, qhK, qhL, N;
	int itp;
	int m, n;
	int i;

   /* Check for proper number of arguments */
  
   if (nrhs != Nargs_rhs) {
     mexErrMsgTxt(PROGNAME " requires " Nargs_rhs_str " input arguments.");
   } else if (nlhs != Nargs_lhs) {
     mexErrMsgTxt(PROGNAME " requires " Nargs_lhs_str " output argument.");
   }
  
   /* Check dimensions for prhs[0], prhs[2], prhs[4]], prhs[5] */

	m = mxGetM(prhs[0]);
	n = mxGetN(prhs[0]); 
	qu = m;
	N = n;
	if (qu<2) mexErrMsgTxt(PROGNAME " requires that qu = size(L_u,1) > 1");		
	if (!mxIsNumeric(prhs[0]) || mxIsComplex(prhs[0]) || 
		 mxIsSparse(prhs[0])  || !mxIsDouble(prhs[0]) ) {
		mexErrMsgTxt(PROGNAME " requires that L_u be a real qu x N matrix.");
	}
  
	m = mxGetM(prhs[4]);
   n = mxGetN(prhs[4]);
   itp = m-3;
   qhK = n; qhL = qhK / qu;
   if ( (qhK%qu) != 0 ) mexErrMsgTxt(PROGNAME " requires that size(State_Table,2) > qu^K.");
	if ( qhK < (qu*qu) ) mexErrMsgTxt(PROGNAME " requires that size(State_Table,2) > qu^2.");
	if ( itp < 1 ) mexErrMsgTxt(PROGNAME " requires that size(State_Table,1) > 2.");
   if ( !mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) || 
        mxIsSparse(prhs[4])  || !mxIsDouble(prhs[4]) ) {
      mexErrMsgTxt(PROGNAME " requires that State_Table be a matrix.");
   }
     
	m = mxGetM(prhs[1]);
	n = mxGetN(prhs[1]); 
	qc = m;
	if (qc<2) mexErrMsgTxt(PROGNAME " requires that qc = size(L_c,1) > 1");		
	if (n != N*itp) mexErrMsgTxt(PROGNAME " requires that size(L_c,2) = N*itp");	
	if (!mxIsNumeric(prhs[1]) || mxIsComplex(prhs[1]) || 
		 mxIsSparse(prhs[1])  || !mxIsDouble(prhs[1]) ) {
		mexErrMsgTxt(PROGNAME " requires that L_c be a real qc x N*itp matrix.");
	}

   m = mxGetM(prhs[2]);
   n = mxGetN(prhs[2]);
   if (    !mxIsNumeric(prhs[2]) || mxIsComplex(prhs[2]) || 
       mxIsSparse(prhs[2])  || !mxIsDouble(prhs[2]) ||
       (  ( (m!=1) || (n!=qhL) ) && ( (n!=1) || (m!=qhL) )  )     ){
      mexErrMsgTxt(PROGNAME " requires that L_Anfangszustand be a qu^L vector.");
   }

   m = mxGetM(prhs[3]);
   n = mxGetN(prhs[3]);
   if (    !mxIsNumeric(prhs[3]) || mxIsComplex(prhs[3]) || 
       mxIsSparse(prhs[3])  || !mxIsDouble(prhs[3]) ||
       (  ( (m!=1) || (n!=qhL) ) && ( (n!=1) || (m!=qhL) )  )     ) {
      mexErrMsgTxt(PROGNAME " requires that L_Endzustand be a qu^L vector.");
	}

	/* Assign pointers and calculate Parameters: r, chips, q, qhL, sigmaq  */
	L_u = mxGetPr(prhs[0]);
	L_c = mxGetPr(prhs[1]);
	L_Anfangszustand = mxGetPr(prhs[2]);
	L_Endzustand = mxGetPr(prhs[3]);
	State_Table = mxGetPr(prhs[4]);

	State_Table_int_dims[0]=itp+3; State_Table_int_dims[1]=qhK;
	State_Table_int_array_ptr = mxCreateNumericArray(2, State_Table_int_dims, mxUINT32_CLASS, mxREAL);
	State_Table_int = (int *)mxGetData(State_Table_int_array_ptr);
    
	/* Test all elements of State_Table to be valid and fill State_Table_int */
	for (n=0; n<qhK; n++)
	{
		m=0;
			i = (int)floor( State_Table[m+(itp+3)*n] + 0.5 );
			if ( (i<1) || (i>qu) )
				{mexPrintf(" Zeile %i, Spalte %i, qc: %i, qu: %i, qhL: %i, Value: %i, Index: %i  \n", m, n,qc,qu,qhL,i,m+qhK*n );
				mexErrMsgTxt("invalid State_Table-elements.");}
			else State_Table_int[m+(itp+3)*n] = i - 1;
		for (m=1; m<(itp+1); m++)
		{
			i = (int)floor( State_Table[m+(itp+3)*n] + 0.5 );
			if ( (i<1) || (i>qc) )
				{mexPrintf(" Zeile %i, Spalte %i, qc: %i, qu: %i, qhL: %i, Value: %i, Index: %i \n", m, n,qc,qu,qhL,i,m+qhK*n );
				mexErrMsgTxt("invalid State_Table-elements.");}
			else State_Table_int[m+(itp+3)*n] = i - 1;
		}
		for (m=(itp+1); m<(itp+3); m++)
		{
			i = (int)floor( State_Table[m+(itp+3)*n] + 0.5 );
			if ( (i<1) || (i>qhL) )
				{mexPrintf(" Zeile %i, Spalte %i, qc: %i, qu: %i, qhL: %i, Value: %i, Index: %i  \n", m, n,qc,qu,qhL,i,m+qhK*n );
				mexErrMsgTxt("invalid State_Table-elements.");}
			else State_Table_int[m+(itp+3)*n] = i - 1;
		}
	}
	
	/* Create a matrix for the return argument */
	plhs[0] = mxCreateDoubleMatrix(qu, N, mxREAL);
	plhs[1] = mxCreateDoubleMatrix(qc, N*itp, mxREAL);

	L_ud = mxGetPr(plhs[0]);
	L_cd = mxGetPr(plhs[1]);
   
	/* Do the actual computations in a subroutine */
	MAP_c( L_ud, L_cd, L_u, L_c, L_Anfangszustand, L_Endzustand, State_Table_int, qu, qc, qhL, qhK, itp, N );

	mxDestroyArray( State_Table_int_array_ptr );
   return;
}


#endif

/******************************************************************************/
/*  MAP-Algorithmus                                                           */
/******************************************************************************/

void MAP_c( double L_ud[], double L_cd[],
            double L_u[], double L_c[], double L_Anfangszustand[], double L_Endzustand[],
            int ST[], const int qu, const int qc,
            const int qhL, const int qhK, const int itp, const int N )

{
	mxArray *alpha_array_ptr, *beta_array_ptr, *beta_s_array_ptr;
	double *alpha, *alpha_s, *beta, *beta_s, *buffer;
	double *alpha_ptr, *alpha_s_ptr, *beta_ptr, *beta_s_ptr;
	double *L_ud_ptr, *L_cd_ptr;
	double L_ud_Summe, L_cd_Summe, f_z;
	int *ST_ptr;
	int j;
	int mu, k, z, i;
	int Nqu = N*qu;
	int Nqcitp = N*qc*itp;
	int NmqhL = N-qhL;
	int Nm1 = N-1;
	/* alpha-, beta-Register */
	alpha_array_ptr = mxCreateDoubleMatrix(qhL, N, mxREAL);	
	alpha = mxGetPr(alpha_array_ptr);
	alpha_s = alpha;	/* Vorg鋘gerknoten */
	alpha = alpha_s + qhL;	/* Nachfolgeknoten */
	beta_array_ptr = mxCreateDoubleMatrix(qhL, 1, mxREAL);
	beta = mxGetPr(beta_array_ptr);
	beta_s_array_ptr = mxCreateDoubleMatrix(qhL, 1, mxREAL);
	beta_s = mxGetPr(beta_s_array_ptr);

	/* Initialisierung der alpha_s- und beta_s-Register */
	alpha_ptr = alpha_s;
	beta_ptr = beta_s;

	for (k=0; k<qhL; k++)
	{
		*(alpha_ptr++) = - *(L_Anfangszustand++);
		*(beta_ptr++) = - *(L_Endzustand++);
	}

	/* L_ud initialisieren (-INF) */
	L_ud_ptr = L_ud;
	for (i=0; i<(Nqu); i++) *(L_ud_ptr++) = -INF;

	/* L_cd initialisieren (-INF) */
	L_cd_ptr = L_cd;
	for (i=0; i<(Nqcitp); i++) *(L_cd_ptr++) = -INF;

	/* alpha initialisieren (INF) */
	alpha_ptr = alpha;
	for (i=0; i<(Nm1*qhL); i++) *(alpha_ptr++) = INF;

	/* MAP-Algorithmus: Vorw鋜tsrekursion */
	/* F黵 alle Symbole */
	for (mu=0; mu<Nm1; mu++)
	{
		ST_ptr = ST;

		/* f黵 alle Zweige z */
		for (z=0; z<qhK; z++)
		{
			/* Apriori-Wahrscheinlichkeit der Inputsymbole ber點ksichtigen */
			f_z = - L_u[ mu*qu + *(ST_ptr++) ];			/* Metrikinkrement pro Zweig */

			/* f黵 alle Codesymbole: Summe der log-躡ergangswahrscheinlichkeiten bestimmen */
			/* (*ST_ptr) = dem Zweig zugeh鰎iges Codesymbol */ 
			for (i=0; i<itp; i++)
				f_z = f_z - L_c[ (mu*itp+i)*qc + *(ST_ptr++) ]; 

			/* Neues alpha aufakkumulieren: Nachfolgeknoten = *(ST_ptr+1), Vorg鋘gerknoten = *(ST_ptr) */
			alpha_s_ptr = alpha_s + *(ST_ptr++);	/* Jetzt zeigt ST_ptr auf die erste Zeile der n鋍hsten Spalte (=n鋍hster Zweig) */
			alpha_ptr = alpha + *(ST_ptr++);

			f_z = f_z + *alpha_s_ptr;
			if ( *alpha_ptr > f_z )
					*alpha_ptr = f_z - log( 1 + exp( f_z - *alpha_ptr ) );
				else
					*alpha_ptr = *alpha_ptr - log( 1 + exp( *alpha_ptr - f_z ) );
/* F黵 Testzwecke! */
/*mexPrintf("mu: %i    alpha: %e      ST_ptr: %i        *ST_ptr: %i\n", mu, exp(- *alpha_ptr), ST_ptr - ST, *ST_ptr );*/
		}
		/*for (z=0;z<qhL;z++) {
			mexPrintf("-%lf, ",*(alpha + z));
        }
		mexPrintf("\n");*/
		/* alpha- und alpha_s-Pointer verschieben */
		alpha_s = alpha;
		alpha = alpha + qhL;
	}

	/*mexPrintf("\n");*/

	/* MAP-Algorithmus: R點kw鋜tsrekursion */
	/* F黵 alle Symbole */
	for (mu=Nm1; mu>=0; mu--)
	{
		/* beta initialisieren/zur點ksetzen (INF) */
		beta_ptr = beta;
		for (k=0; k<qhL; k++) *(beta_ptr++) = INF;

		ST_ptr = ST;

		/* f黵 alle Zweige z */
		for (z=0; z<qhK; z++)
		{
			/* Apriori-Wahrscheinlichkeit der Inputsymbole ber點ksichtigen */
			f_z = - L_u[ mu*qu + *(ST_ptr++) ];		/* Metrikinkrement pro Zweig */

			/* f黵 alle Codesymbole: Summe der log-躡ergangswahrscheinlichkeiten bestimmen */
			/* (*ST_ptr) = dem Zweig zugeh鰎iges Codesymbol */ 
			for (i=0; i<itp; i++)
				f_z = f_z - L_c[ (mu*itp+i)*qc + *(ST_ptr++) ]; 

			/* neues beta aufakkumulieren */
			alpha_ptr = alpha_s + *(ST_ptr);
			beta_ptr = beta + *(ST_ptr++);
			beta_s_ptr = beta_s + *(ST_ptr++);	/* Jetzt zeigt ST_ptr auf die erste Zeile der n鋍hsten Spalte (=n鋍hster Zweig) */
			
			f_z = f_z + *beta_s_ptr;
			if ( *beta_ptr > f_z )
					*beta_ptr = f_z - log( 1 + exp( f_z - *beta_ptr ) );
				else
					*beta_ptr = *beta_ptr - log( 1 + exp( *beta_ptr - f_z ) );
/* F黵 Testzwecke! */
/*mexPrintf("mu: %i    beta: %e\n", mu, exp(- *beta_ptr) );*/

			ST_ptr = ST_ptr - itp - 3;		/* State_Table-Pointer nochmal zur點k um die zugeh鰎igen input und output Symbole zu adressieren */
			f_z = f_z + *alpha_ptr;		/* -f_z ist jetzt proportional der Wahrscheinlichkeit des zugeh鰎igen Zweiges */
			
			/* log-Wahrscheinlichkeit der Inputsymbole L_ud akkumulieren */
			j = *(ST_ptr++);
			if ( L_ud[ mu*qu + j ] > -f_z)
					L_ud[ mu*qu + j ] = L_ud[ mu*qu + j ] + log( 1 + exp( -f_z - L_ud[ mu*qu + j ] ) );
				else
					L_ud[ mu*qu + j ] = -f_z + log( 1 + exp( L_ud[ mu*qu + j ] + f_z ) );

			/* log-Wahrscheinlichkeit der Codesymbole L_cd akkumulieren */
			for (i=0; i<itp; i++)
			{
				j = *(ST_ptr++);
				if ( L_cd[ (mu*itp+i)*qc + j ] > -f_z)
						L_cd[ (mu*itp+i)*qc + j ] = L_cd[ (mu*itp+i)*qc + j ] + log( 1 + exp( -f_z - L_cd[ (mu*itp+i)*qc + j ] ) );
					else
						L_cd[ (mu*itp+i)*qc + j ] = -f_z + log( 1 + exp( L_cd[ (mu*itp+i)*qc + j ] + f_z ) );
			}

			ST_ptr = ST_ptr + 2;	/* State_Table-Pointer auf die naechste Spalte (Zweig) setzen */
		}
		/*
		for (z=0;z<qhL;z++) {
			mexPrintf("-%lf, ",*(beta + z));
        }
		mexPrintf("\n");
		*/
		/* beta- und beta_s-Pointer vertauschen (double-buffer) und alpha_s neu setzen*/
		buffer = beta_s;
		beta_s = beta;
		beta = buffer;
		alpha_s = alpha_s - qhL;



		/* L_ud normieren, damit sum(exp(L_ud))=1 */
		L_ud_Summe = L_ud[ mu*qu ];
		for (z=1;z<qu;z++)
			if ( L_ud_Summe > L_ud[ mu*qu + z ] )
				L_ud_Summe = L_ud_Summe + log( 1 + exp( L_ud[ mu*qu + z ] - L_ud_Summe ) );
			else
				L_ud_Summe = L_ud[ mu*qu + z ] + log( 1 + exp( L_ud_Summe - L_ud[ mu*qu + z ] ) );
		for (z=0;z<qu;z++)
			L_ud[ mu*qu + z ] = L_ud[ mu*qu + z ] - L_ud_Summe;

		/* L_cd normieren, damit sum(exp(L_cd))=1 */
		for (i=0; i<itp; i++)
		{
			L_cd_Summe = L_cd[ (mu*itp+i)*qc ];
			for (z=1;z<qc;z++)
				if ( L_cd_Summe > L_cd[ (mu*itp+i)*qc + z ] )
					L_cd_Summe = L_cd_Summe + log( 1 + exp( L_cd[ (mu*itp+i)*qc + z ] - L_cd_Summe ) );
				else
					L_cd_Summe = L_cd[ (mu*itp+i)*qc + z ] + log( 1 + exp( L_cd_Summe - L_cd[ (mu*itp+i)*qc + z ] ) );
			for (z=0;z<qc;z++)
				L_cd[ (mu*itp+i)*qc + z ] = L_cd[ (mu*itp+i)*qc + z ] - L_cd_Summe;
		}
	}

	
	mxDestroyArray(beta_s_array_ptr);
	mxDestroyArray(beta_array_ptr);
	mxDestroyArray(alpha_array_ptr);

}




/* TEST ****************************************************************************************