power spectral density for UW-OFDM systems

I am having code in C language. How to convert this C program into .m file? I have attached my C program with this. please help me to convert this code. Thanks in advance.

#include <math.h>

#include "mex.h"

#include "matrix.h"

#include "complexd.c"

#include <omp.h>

#define INF 1e300

#define DEFAULT_MAX_NODEVISITS 1e11

#define DEFAULT_L_MAX 5

#define ERRMSGPATH UWOFDM:sd_worker_sts


/**
 
* SOFT-OUTPUT SPHERE DECODER
 
*
 
* Alexander Onic [e-mail]
 
* Alpen-Adria-Universität Klagenfurt, Austria
 
*

* Matlab call:
 
* [d_hat, L] = sd_worker_sts(r, H, con, max_nv, L_max)
 
*
 
* d_hat	... decoding hypotheses, one OFDM symbol per column (symbol numbers acc. to alphabet vector 'con')
 
* L		... LLRs for each bit, one OFDM symbol per column
 
* r		... receive vectors r = y' = Q^H y [1]

 * H		... channel propagation matrix H = tilde{H} G [1]
 
* con		... constellation vector = transmit alphabet
 
* [max_nv]	... maximum number of node visits per OFDM symbol, if reached detection of current OFDM symbol is stopped (optimum: inf)

* [L_max]	... LLR clipping level (optimum: inf)
 
*/
 

struct qam_metric {
	
        int symbol;					/* symbol number */
	
struct complexd diff;		/* difference of symbol from receive value */
	
double metric;				/* metric = diff^2, Euclidean distance */

};


void decision_feedback(const struct complexd *e, const struct complexd *y, const struct complexd *H, const unsigned int k, 
struct complexd *out) {

	unsigned int r;

	for (r=0; r<k; r++) {
		
             out[r].re = e[r].re - y->re*H[r].re + y->im*H[r].im;
		
             out[r].im = e[r].im - y->re*H[r].im - y->im*H[r].re;
	
}

}


/**

 * Matrix multiplication with a vector H*x

 * dimensions (Nd x Nd) and (Nd x 1)
 
*/

void mul_matvec(const struct complexd *H, const struct complexd *x, const unsigned int Nd, struct complexd *out) {

	unsigned int r, c;
	
for (r=0; r<Nd; r++) {
	
	out[r].re = 0;
	
	out[r].im = 0;
		
for (c=0; c<Nd; c++) {
	
		out[r].re += H[r+c*Nd].re * x[c].re - H[r+c*Nd].im * x[c].im;

	        out[r].im += H[r+c*Nd].re * x[c].im + H[r+c*Nd].im * x[c].re;
		
}
	
}

}


/**
 
* Translates the decimal number x into a binary of length nBits.
 
* LSB last.

 */

void dec2bin(const int x, const unsigned int nBits, char
 *out) {

	unsigned int bit;
	
for (bit=0; bit<nBits; bit++)
		
out[nBits-bit-1] = x>>bit & 1;

}


/**

 * Creates, sorts and returns a list of the symbols of con minus

 * receive value e with its difference and the metric. The list is

 * sorted in ascending order by the difference.

 */
void initList(

		const struct complexd con[], const struct complexd e, const struct complexd H, const int M,

	struct qam_metric list_k[])

{
	struct complexd curdist;
	
int sym;

	
/* insert each symbol number 'sym' */
	
for (sym=0; sym<M; sym++) {
		
/* Determine distance of symbol to equalized receive symbol */
		
curdist.re = ((e.re - (con+sym)->re)*H.re + (e.im - (con+sym)->im)*H.im) / cmplx_abs2(&H);

curdist.im = ((e.im - (con+sym)->im)*H.re - (e.re - (con+sym)->re)*H.im) / cmplx_abs2(&H);

		
/* insert symbol into unsorted list */
	
	list_k[sym].symbol = sym;
		
list_k[sym].diff = curdist;
		
list_k[sym].metric = cmplx_abs2(&curdist);

	}

}


/**
 
* Return the next symbol from the sorted list. NULL if list
 
* is empty.
 
*/
void getNextSymbol(struct qam_metric list_k[], const int M, struct qam_metric **best_sym)

{
	
double best_metric = INF;
	
int dum;
	
*best_sym = NULL;

	
/* Search for entry with lowest metric in list */
	
for (dum=0; dum<M; dum++)
		
if (list_k[dum].metric < best_metric) 
{
			
*best_sym = list_k+dum;
			
best_metric = (*best_sym)->metric;
	
}
}


/**
 
* MAIN FUNCTION, SOFT-OUTPUT SPHERE DECODER
 
*/
void sd_decode(
		
const struct complexd x[], const struct complexd H[],
		
const struct complexd constellation[],
		
const int M, const unsigned int Nd, const unsigned long nSym,


const unsigned long max_nodevisits, const double L_max,
	
	int *rxSymbols, double *rxLLRs)

{
	
unsigned int		k, r;
	
double				metric_best, radius, *LLRs, *metric;

	struct qam_metric	*new_symbol = NULL, *metriclist;

	struct complexd		*e;

	unsigned long		iSym, nodevisits;

	int					*un, bps = ceil(log2(M));

	char				*un_bits, *un_bits_best;


#pragma omp parallel private(iSym, k, r, metric_best, radius, nodevisits, un, un_bits, un_bits_best, LLRs, metric, new_symbol, e, metriclist)

{
	
un = calloc(Nd, sizeof(int));
	
un_bits = calloc(Nd*bps, sizeof(char));

un_bits_best = calloc(Nd*bps, sizeof(char));
	
metriclist  = calloc(M*Nd, sizeof(struct qam_metric));
	
e = calloc(Nd*Nd, sizeof(struct complexd));
	
metric = calloc(Nd+1, sizeof(double));
	
LLRs = calloc(Nd*bps, sizeof(double));

	
#pragma omp for
	for (iSym = 0; iSym<nSym; iSym++) 
{
		
metric_best = INF;
		
k = Nd-1;
		
nodevisits = max_nodevisits;

		
for (r = 0; r<Nd*bps; r++) 
{
			
un_bits[r] = 0;
			
un_bits_best[r] = -1;
			
LLRs[r] = INF;
		
}

		
/* Matrix multiplication: Determine e = R^-1 * y */
	
mul_matvec(H, x+iSym*Nd, Nd, e+k*Nd);

		
/* Initial symbol and metrics list for top level */
	
	initList(constellation, e[k+k*Nd], H[k+k*Nd], M, metriclist+k*M);


		while (k < Nd && nodevisits > 0) {

			getNextSymbol(metriclist+k*M, M, &new_symbol);

			if (new_symbol == NULL) /* all symbols tried -> up
 */
				
k++;
		
	
else {
			
	nodevisits--;

		
		/* Update current vector and bit pattern. */
		
		un[k] = new_symbol->symbol;
				
for (r=k; r<Nd; r++)
	
				dec2bin(un[r], bps, un_bits + r*bps);

	
			metric[k] = metric[k+1] + new_symbol->metric;
	
			new_symbol->metric = 2*INF;

			
	/* update search radius = tree pruning */
		
		radius = 0;
				
for (r=0; r<Nd*bps; r++)
	
if (radius < LLRs[r] && (r<(k+1)*bps || un_bits[r] != un_bits_best[r]))
		
radius = LLRs[r];

			
if (metric[k] <= radius) 
{

	
if (k > 0) {
					
	/* e - u/h
						* [2] alg Decode() line 13 */
						decision_feedback(e+k*Nd, &new_symbol->diff, H+k*Nd, k, e+(k-1)*Nd);
						k--;

						/* Create new symbol and metrics list for the current level k */
						initList(constellation, e[k+k*Nd], H[k+k*Nd], M, metriclist+k*M);
					} else {	/* reached a leaf */
						if (metric[k] < metric_best) {
							if (un_bits_best[0] == -1)
								for (r=0; r<Nd*bps; r++)
									un_bits_best[r] = un_bits[r];

							for (r=0; r<Nd*bps; r++) {
								if (un_bits[r] != un_bits_best[r])
									LLRs[r] = metric[k];
								un_bits_best[r] = un_bits[r];
							}
							for (r=0; r<Nd; r++)
								rxSymbols[r+iSym*Nd] = un[r];
							metric_best = metric[k];
						} else {
							for (r=0; r<Nd*bps; r++)
								if (un_bits[r] != un_bits_best[r])
									if (metric[k] < LLRs[r])
										LLRs[r] = metric[k];
						}

						/* LLR clipping */
						for (r=0; r<Nd*bps; r++)
							if (metric_best + L_max < LLRs[r])
								LLRs[r] = metric_best + L_max;
					}
				} else {
					/* only check number of node visits when going up a level */
					if (nodevisits <= 0)
						k = Nd;
					else
						k++;
				}
			}
		}

		/* Store result in output variables */
		for (r=0; r<Nd*bps; r++)
			rxLLRs[iSym*Nd*bps + r] = (metric_best - LLRs[r] - 1e-9) * (un_bits_best[r]*2-1);
		/** Subtraction of 1e-9:
		 * We need to make sure rxLLRs is not 0, especially for the case
		 * L_max = 0.
		 */
	}

	free(un);
	free(un_bits);
	free(un_bits_best);
	free(metric);
	free(metriclist);
	free(e);
	free(LLRs);
}
}

void mexFunction(
	int nlhs, mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	const mxArray	*m_y, *m_H, *m_con;
	unsigned int	r, c, Nd;
	int				*rxSymbols, bps, M;
	unsigned long	nSym, iSym, max_nodevisits;
	double			*m_real, *m_imag, *rxLLRs, L_max;
	struct complexd	*y, *H, *con;

	if (nrhs < 3 || nrhs > 5)
		mexErrMsgIdAndTxt("ERRMSGPATH:nrhs", "Three to five inputs required.");
	if (nlhs > 2)
		mexErrMsgIdAndTxt("ERRMSGPATH:nlhs", "One or two outputs required.");

	m_y = prhs[0];		/* receive vectors in columns */
	m_H = prhs[1];		/* channel propagation matrix */
	m_con = prhs[2];	/* transmit symbol alphabet */
	Nd = mxGetM(m_H);	/* number of values per OFDM symbol */
	nSym = mxGetN(m_y);	/* number of OFDM symbols per burst */

	if (mxGetM(m_y) != mxGetM(m_H) || mxGetM(m_H) != mxGetN(m_H))
		mexErrMsgIdAndTxt("ERRMSGPATH:input_dim", "Matrix dimensions of input y and H must agree. H must be square.");

	/*** copy Matlab data ***/
	/* Process receive vectors */
	y = calloc(Nd*nSym, sizeof(struct complexd));
	m_real = mxGetPr(m_y);
	m_imag = mxGetPi(m_y);
	for (iSym=0; iSym<nSym; iSym++) {
		for (r=0; r<Nd; r++) {
			y[r+iSym*Nd].re = m_real[r+iSym*Nd];
			y[r+iSym*Nd].im = m_imag[r+iSym*Nd];
		}
	}

	/* Process channel propagation matrix */
	m_real = mxGetPr(m_H);
	m_imag = mxGetPi(m_H);
	if (m_imag == NULL) {
		free(y);
		mexErrMsgIdAndTxt("ERRMSGPATH:H_noimag", "Matrix H has no imaginary part or could not be retrieved from Matlab.");
	}
	H = calloc(Nd*Nd, sizeof(struct complexd));
	for (c=0; c<Nd; c++) {
		for (r=0; r<Nd; r++) {
			H[r+c*Nd].re = m_real[r+c*Nd];
			H[r+c*Nd].im = m_imag[r+c*Nd];
		}
	}
	/* Ensure upper triangular structure of H */
	for (c=0; c<Nd; c++) {
		for (r=c+1; r<Nd; r++) {
			if (!cmplx_iszero(H+(r + Nd*c))) {
				free(y);
				free(H);
				mexErrMsgIdAndTxt("ERRMSGPATH:H_not_triangular", "H must be upper triangular.");
			}
		}
	}

	/* Process constellation vector */
	M = mxGetNumberOfElements(m_con);	/* constellation size */
	bps = ceil(log2(M));				/* bits per symbol */
	con = calloc(M, sizeof(struct complexd));
	m_real = mxGetPr(m_con);
	m_imag = mxGetPi(m_con);
	for (r = 0; r<M; r++) {
		con[r].re = m_real[r];
		con[r].im = m_imag[r];
	}

	/* maximum number of node visits and L_max */
	if (nrhs >= 4)
		max_nodevisits = (unsigned long) mxGetScalar(prhs[3]);
	else
		max_nodevisits = DEFAULT_MAX_NODEVISITS;
	if (nrhs >= 5)
		L_max = mxGetScalar(prhs[4]);
	else
		L_max = DEFAULT_L_MAX;

	/* prepare output variables */
	plhs[0] = mxCreateNumericMatrix(Nd, nSym, mxINT32_CLASS, mxREAL);
	rxSymbols = mxGetData(plhs[0]);

	plhs[1] = mxCreateNumericMatrix(Nd*bps, nSym, mxDOUBLE_CLASS, mxREAL);
	rxLLRs = mxGetData(plhs[1]);

	sd_decode(y, H, con, M, Nd, nSym, max_nodevisits, L_max, rxSymbols, rxLLRs);

	free(y);
	free(H);
	free(con);
}