power spectral density for UW-OFDM systems

1 次查看（过去 30 天）
Shan Sha
Shan Sha 2019-2-8
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 alexander.onic@aau.at
* 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);
}

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Modulation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by