ZF/MMSE equalizer for single antenna and diversity working correctly

master
ismagom 10 years ago
parent 4f653dc3ab
commit 68a6169164

@ -126,10 +126,6 @@ int find_cell(void *uhd, ue_celldetect_result_t *found_cell, uint32_t N_id_2)
goto free_and_exit;
}
ue_celldetect_set_nof_frames_detected(&cd, 50);
ue_celldetect_set_nof_frames_total(&cd, 500);
INFO("Setting sampling frequency 960 KHz for PSS search\n", 0);
cuhd_set_rx_srate(uhd, 960000.0);
INFO("Starting receiver...\n", 0);
@ -223,6 +219,10 @@ int cell_search(void *uhd, int force_N_id_2, lte_cell_t *cell)
ue_celldetect_result_t found_cells[3];
bzero(found_cells, 3*sizeof(ue_celldetect_result_t));
ue_celldetect_set_nof_frames_detected(&cd, 50);
ue_celldetect_set_nof_frames_total(&cd, 500);
if (force_N_id_2 >= 0) {
ret = find_cell(uhd, &found_cells[force_N_id_2], force_N_id_2);
} else {

@ -40,6 +40,9 @@ typedef _Complex float cf_t;
typedef struct {
cf_t *h_mod;
cf_t *tmp1;
cf_t *tmp2;
cf_t *tmp3;
float *y_mod;
float *z_real;
float *z_imag;
@ -74,25 +77,20 @@ LIBLTE_API int precoding_type(precoding_t *q,
/* Estimates the vector "x" based on the received signal "y" and the channel estimates "h"
*/
LIBLTE_API int predecoding_single_zf(precoding_t *q,
LIBLTE_API int predecoding_single(precoding_t *q,
cf_t *y,
cf_t *h,
cf_t *x,
int nof_symbols,
float noise_estimate);
LIBLTE_API int predecoding_diversity(precoding_t *q,
cf_t *y,
cf_t *h,
cf_t *x,
int nof_symbols);
LIBLTE_API int predecoding_single_mmse(precoding_t *q,
cf_t *y,
cf_t *h,
cf_t *x,
int nof_symbols,
float noise_estimate);
LIBLTE_API int predecoding_diversity_zf(precoding_t *q,
cf_t *y,
cf_t *h[MAX_PORTS],
cf_t *x[MAX_LAYERS],
int nof_ports,
int nof_symbols);
cf_t *h[MAX_PORTS],
cf_t *x[MAX_LAYERS],
int nof_ports,
int nof_symbols,
float noise_estimate);
LIBLTE_API int predecoding_type(precoding_t *q,
cf_t *y,
@ -101,6 +99,7 @@ LIBLTE_API int predecoding_type(precoding_t *q,
int nof_ports,
int nof_layers,
int nof_symbols,
lte_mimo_type_t type);
lte_mimo_type_t type,
float noise_estimate);
#endif /* PRECODING_H_ */

@ -107,6 +107,7 @@ LIBLTE_API float vec_dot_prod_fff(float *x, float *y, uint32_t len);
/* z=x/y vector division (element-wise) */
LIBLTE_API void vec_div_ccc(cf_t *x, cf_t *y, float *y_mod, cf_t *z, float *z_real, float *z_imag, uint32_t len);
void vec_div_cfc(cf_t *x, float *y, cf_t *z, float *z_real, float *z_imag, uint32_t len);
LIBLTE_API void vec_div_fff(float *x, float *y, float *z, uint32_t len);
/* conjugate */

@ -173,7 +173,7 @@ int main(int argc, char **argv) {
gettimeofday(&t[1], NULL);
for (int j=0;j<100;j++) {
predecoding_single_zf(&cheq, input, ce, output, num_re);
predecoding_single(&cheq, input, ce, output, num_re, 0);
}
gettimeofday(&t[2], NULL);
get_time_interval(t);
@ -188,7 +188,7 @@ int main(int argc, char **argv) {
gettimeofday(&t[1], NULL);
for (int j=0;j<100;j++) {
predecoding_single_mmse(&cheq, input, ce, output, num_re, chest_dl_get_noise_estimate(&est));
predecoding_single(&cheq, input, ce, output, num_re, chest_dl_get_noise_estimate(&est));
}
gettimeofday(&t[2], NULL);
get_time_interval(t);

@ -30,10 +30,10 @@
#define CELLID prhs[0]
#define PORTS prhs[1]
#define INPUT prhs[2]
#define NOF_INPUTS 3
#define SFIDX prhs[3]
#define FREQ_FILTER prhs[3]
#define TIME_FILTER prhs[4]
#define NOF_INPUTS 5
#define SFIDX prhs[5]
void help()
{
@ -56,7 +56,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
lte_cell_t cell;
chest_dl_t chest;
precoding_t cheq;
cf_t *input_signal = NULL, *output_signal = NULL;
cf_t *input_signal = NULL, *output_signal[MAX_LAYERS];
cf_t *output_signal2 = NULL;
cf_t *ce[MAX_PORTS];
double *outr0=NULL, *outi0=NULL;
double *outr1=NULL, *outi1=NULL;
@ -107,32 +108,37 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
return;
}
sf_idx = (uint32_t) *((double*) mxGetPr(SFIDX));
} else {
if (nrhs != NOF_INPUTS) {
help();
return;
}
}
uint32_t filter_len = 0;
float *filter;
double *f;
if (nrhs > NOF_INPUTS && nsubframes == 10) {
if (nrhs >= NOF_INPUTS + 1) {
filter_len = mxGetNumberOfElements(FREQ_FILTER);
filter = malloc(sizeof(float) * filter_len);
f = (double*) mxGetPr(FREQ_FILTER);
for (i=0;i<filter_len;i++) {
filter[i] = (float) f[i];
}
chest_dl_set_filter_freq(&chest, filter, filter_len);
}
if (nrhs >= NOF_INPUTS + 2) {
filter_len = mxGetNumberOfElements(TIME_FILTER);
filter = malloc(sizeof(float) * filter_len);
f = (double*) mxGetPr(TIME_FILTER);
for (i=0;i<filter_len;i++) {
filter[i] = (float) f[i];
}
chest_dl_set_filter_time(&chest, filter, filter_len);
}
filter_len = mxGetNumberOfElements(FREQ_FILTER);
filter = malloc(sizeof(float) * filter_len);
f = (double*) mxGetPr(FREQ_FILTER);
for (i=0;i<filter_len;i++) {
filter[i] = (float) f[i];
}
chest_dl_set_filter_freq(&chest, filter, filter_len);
filter_len = mxGetNumberOfElements(TIME_FILTER);
filter = malloc(sizeof(float) * filter_len);
f = (double*) mxGetPr(TIME_FILTER);
for (i=0;i<filter_len;i++) {
filter[i] = (float) f[i];
}
chest_dl_set_filter_time(&chest, filter, filter_len);
double *inr=(double *)mxGetPr(INPUT);
double *ini=(double *)mxGetPi(INPUT);
@ -143,14 +149,16 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
ce[i] = vec_malloc(nof_re * sizeof(cf_t));
}
input_signal = vec_malloc(nof_re * sizeof(cf_t));
output_signal = vec_malloc(nof_re * sizeof(cf_t));
for (i=0;i<MAX_PORTS;i++) {
output_signal[i] = vec_malloc(nof_re * sizeof(cf_t));
}
output_signal2 = vec_malloc(nof_re * sizeof(cf_t));
precoding_init(&cheq, nof_re);
/* Create output values */
if (nlhs >= 1) {
plhs[0] = mxCreateDoubleMatrix(1,nof_re * nsubframes, mxCOMPLEX);
plhs[0] = mxCreateDoubleMatrix(nof_re * nsubframes, cell.nof_ports, mxCOMPLEX);
outr0 = mxGetPr(plhs[0]);
outi0 = mxGetPi(plhs[0]);
}
@ -160,7 +168,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
outi1 = mxGetPi(plhs[1]);
}
if (nlhs >= 3) {
plhs[2] = mxCreateDoubleMatrix(1,nof_re * nsubframes, mxCOMPLEX);
plhs[2] = mxCreateDoubleMatrix(nof_re * nsubframes, 1, mxCOMPLEX);
outr2 = mxGetPr(plhs[2]);
outi2 = mxGetPi(plhs[2]);
}
@ -185,17 +193,23 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt("Error running channel estimator\n");
return;
}
//predecoding_single_zf(&cheq, input_signal, ce[0], output_signal, nof_re);
predecoding_single_mmse(&cheq, input_signal, ce[0], output_signal, nof_re, chest_dl_get_noise_estimate(&chest));
if (cell.nof_ports == 1) {
predecoding_single(&cheq, input_signal, ce[0], output_signal2, nof_re, chest_dl_get_noise_estimate(&chest));
} else {
predecoding_diversity(&cheq, input_signal, ce, output_signal, cell.nof_ports, nof_re, chest_dl_get_noise_estimate(&chest));
layerdemap_diversity(output_signal, output_signal2, cell.nof_ports, nof_re/cell.nof_ports);
}
if (nlhs >= 1) {
for (i=0;i<nof_re;i++) {
*outr0 = (double) crealf(ce[0][i]);
if (outi0) {
*outi0 = (double) cimagf(ce[0][i]);
for (int j=0;j<cell.nof_ports;j++) {
for (i=0;i<nof_re;i++) {
*outr0 = (double) crealf(ce[j][i]);
if (outi0) {
*outi0 = (double) cimagf(ce[j][i]);
}
outr0++;
outi0++;
}
outr0++;
outi0++;
}
}
if (nlhs >= 2) {
@ -212,9 +226,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
if (nlhs >= 3) {
for (i=0;i<nof_re;i++) {
*outr2 = (double) crealf(output_signal[i]);
*outr2 = (double) crealf(output_signal2[i]);
if (outi2) {
*outi2 = (double) cimagf(output_signal[i]);
*outi2 = (double) cimagf(output_signal2[i]);
}
outr2++;
outi2++;

@ -51,6 +51,21 @@ int precoding_init(precoding_t *q, uint32_t max_frame_len) {
perror("malloc");
goto clean_exit;
}
q->tmp1 = vec_malloc(sizeof(cf_t) * max_frame_len);
if (!q->tmp1) {
perror("malloc");
goto clean_exit;
}
q->tmp2 = vec_malloc(sizeof(cf_t) * max_frame_len);
if (!q->tmp2) {
perror("malloc");
goto clean_exit;
}
q->tmp3 = vec_malloc(sizeof(cf_t) * max_frame_len);
if (!q->tmp3) {
perror("malloc");
goto clean_exit;
}
q->y_mod = vec_malloc(sizeof(float) * max_frame_len);
if (!q->y_mod) {
perror("malloc");
@ -77,6 +92,16 @@ clean_exit:
}
void precoding_free(precoding_t *q) {
if (q->tmp1) {
free(q->tmp1);
}
if (q->tmp2) {
free(q->tmp2);
}
if (q->tmp3) {
free(q->tmp3);
}
if (q->h_mod) {
free(q->h_mod);
}
@ -92,54 +117,83 @@ void precoding_free(precoding_t *q) {
bzero(q, sizeof(precoding_t));
}
/* ZF SISO equalizer: x=y/h */
int predecoding_single_zf(precoding_t *q, cf_t *y, cf_t *h, cf_t *x, int nof_symbols) {
if (nof_symbols <= q->max_frame_len) {
vec_div_ccc(y, h, q->y_mod, x, q->z_real, q->z_imag, nof_symbols);
return nof_symbols;
} else {
return LIBLTE_ERROR;
}
}
/* MMSE SISO equalizer x=y*h'/(h*h'+no) */
int predecoding_single_mmse(precoding_t *q, cf_t *y, cf_t *h, cf_t *x, int nof_symbols, float noise_estimate) {
/* ZF/MMSE SISO equalizer x=y(h'h+no)^(-1)h' (ZF if n0=0.0)*/
int predecoding_single(precoding_t *q, cf_t *y, cf_t *h, cf_t *x, int nof_symbols, float noise_estimate) {
if (nof_symbols <= q->max_frame_len) {
// h*h'
vec_prod_conj_ccc(h, h, q->h_mod, nof_symbols);
// real(h*h')
vec_deinterleave_real_cf(q->h_mod, q->y_mod, nof_symbols);
// (h*h' + n0)
vec_sc_add_fff(q->y_mod, noise_estimate, q->y_mod, nof_symbols);
// h'h
vec_abs_square_cf(h, q->y_mod, nof_symbols);
if (noise_estimate > 0.0) {
// (h'h + n0)
vec_sc_add_fff(q->y_mod, noise_estimate, q->y_mod, nof_symbols);
}
// y*h'
vec_prod_conj_ccc(y, h, x, nof_symbols);
// decompose real/imag parts
vec_deinterleave_cf(x, q->z_real, q->z_imag, nof_symbols);
// real and imag division
vec_div_fff(q->z_real, q->y_mod, q->z_real, nof_symbols);
vec_div_fff(q->z_imag, q->y_mod, q->z_imag, nof_symbols);
// interleave again
vec_interleave_cf(q->z_real, q->z_imag, x, nof_symbols);
// divide by (h'h+no)
vec_div_cfc(x,q->y_mod,x,q->z_real,q->z_imag, nof_symbols);
return nof_symbols;
} else {
return LIBLTE_ERROR;
}
}
/* ZF STBC equalizer */
int predecoding_diversity_zf(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *x[MAX_LAYERS],
int nof_ports, int nof_symbols) {
/* ZF/MMSE STBC equalizer x=y(H'H+n0·I)^(-1)H' (ZF is n0=0.0)
*/
int predecoding_diversity(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *x[MAX_LAYERS], int nof_ports, int nof_symbols, float noise_estimate) {
int i;
if (nof_ports == 2) {
#define new
#ifdef new
// reuse buffers
cf_t *r0 = q->tmp3;
cf_t *r1 = &q->tmp3[nof_symbols/2];
cf_t *h0 = q->h_mod;
cf_t *h1 = &q->h_mod[nof_symbols/2];
float *modhh = q->y_mod;
float *modh0 = q->z_real;
float *modh1 = q->z_imag;
// prepare buffers
for (i=0;i<nof_symbols/2;i++) {
h0[i] = h[0][2*i]; // h0
h1[i] = h[1][2*i+1]; // h1
r0[i] = y[2*i]; // r0
r1[i] = y[2*i+1]; // r1
}
// Compute common dividend and store in y_mod
vec_abs_square_cf(h0, modh0, nof_symbols/2);
vec_abs_square_cf(h1, modh1, nof_symbols/2);
vec_sum_fff(modh0, modh1, modhh, nof_symbols/2);
if (noise_estimate > 0.0) {
// (H'H + n0)
vec_sc_add_fff(modhh, noise_estimate, modhh, nof_symbols/2);
}
vec_sc_prod_fff(modhh, 1.0/sqrt(2), modhh, nof_symbols/2);
// x[0] = r0·h0*/(|h0|+|h1|)+r1*·h1/(|h0|+|h1|)
vec_prod_conj_ccc(r0,h0,q->tmp1, nof_symbols/2);
vec_prod_conj_ccc(h1,r1,q->tmp2, nof_symbols/2);
vec_sum_ccc(q->tmp1, q->tmp2, x[0], nof_symbols/2);
vec_div_cfc(x[0], modhh, x[0], q->z_real, q->z_imag, nof_symbols/2);
// x[1] = r1·h0*/(|h0|+|h1|)-r0*·h1/(|h0|+|h1|)
vec_prod_conj_ccc(r1,h0,q->tmp1, nof_symbols/2);
vec_prod_conj_ccc(h1,r0,q->tmp2, nof_symbols/2);
vec_sub_ccc(q->tmp1, q->tmp2, x[1], nof_symbols/2);
vec_div_cfc(x[1], modhh, x[1], q->z_real, q->z_imag, nof_symbols/2);
#else
cf_t h0, h1, h2, h3, r0, r1, r2, r3;
float hh, hh02, hh13;
if (nof_ports == 2) {
/* TODO: Use VOLK here */
for (i = 0; i < nof_symbols / 2; i++) {
h0 = h[0][2 * i];
h1 = h[1][2 * i];
hh = crealf(h0) * crealf(h0) + cimagf(h0) * cimagf(h0)
+ crealf(h1) * crealf(h1) + cimagf(h1) * cimagf(h1);
+ crealf(h1) * crealf(h1) + cimagf(h1) * cimagf(h1) + noise_estimate;
r0 = y[2 * i];
r1 = y[2 * i + 1];
if (hh == 0) {
@ -148,8 +202,11 @@ int predecoding_diversity_zf(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *
x[0][i] = (conjf(h0) * r0 + h1 * conjf(r1)) / hh * sqrt(2);
x[1][i] = (-h1 * conj(r0) + conj(h0) * r1) / hh * sqrt(2);
}
#endif
return i;
} else if (nof_ports == 4) {
cf_t h0, h1, h2, h3, r0, r1, r2, r3;
float hh02, hh13;
int m_ap = (nof_symbols % 4) ? ((nof_symbols - 2) / 4) : nof_symbols / 4;
for (i = 0; i < m_ap; i++) {
@ -181,7 +238,7 @@ int predecoding_diversity_zf(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *
/* 36.211 v10.3.0 Section 6.3.4 */
int predecoding_type(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *x[MAX_LAYERS],
int nof_ports, int nof_layers, int nof_symbols, lte_mimo_type_t type) {
int nof_ports, int nof_layers, int nof_symbols, lte_mimo_type_t type, float noise_estimate) {
if (nof_ports > MAX_PORTS) {
fprintf(stderr, "Maximum number of ports is %d (nof_ports=%d)\n", MAX_PORTS,
@ -197,7 +254,7 @@ int predecoding_type(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *x[MAX_LA
switch (type) {
case SINGLE_ANTENNA:
if (nof_ports == 1 && nof_layers == 1) {
return predecoding_single_zf(q, y, h[0], x[0], nof_symbols);
return predecoding_single(q, y, h[0], x[0], nof_symbols, noise_estimate);
} else {
fprintf(stderr,
"Number of ports and layers must be 1 for transmission on single antenna ports\n");
@ -206,7 +263,7 @@ int predecoding_type(precoding_t *q, cf_t *y, cf_t *h[MAX_PORTS], cf_t *x[MAX_LA
break;
case TX_DIVERSITY:
if (nof_ports == nof_layers) {
return predecoding_diversity_zf(q, y, h, x, nof_ports, nof_symbols);
return predecoding_diversity(q, y, h, x, nof_ports, nof_symbols, noise_estimate);
} else {
fprintf(stderr,
"Error number of layers must equal number of ports in transmit diversity\n");
@ -239,13 +296,15 @@ int precoding_diversity(precoding_t *q, cf_t *x[MAX_LAYERS], cf_t *y[MAX_PORTS],
int nof_symbols) {
int i;
if (nof_ports == 2) {
/* FIXME: Use VOLK here */
for (i = 0; i < nof_symbols; i++) {
y[0][2 * i] = x[0][i] / sqrtf(2);
y[1][2 * i] = -conjf(x[1][i]) / sqrtf(2);
y[0][2 * i + 1] = x[1][i] / sqrtf(2);
y[1][2 * i + 1] = conjf(x[0][i]) / sqrtf(2);
y[0][2 * i] = x[0][i];
y[1][2 * i] = -conjf(x[1][i]);
y[0][2 * i + 1] = x[1][i];
y[1][2 * i + 1] = conjf(x[0][i]);
}
// normalize
vec_sc_prod_cfc(y[0], 1.0/sqrtf(2), y[0], 2*nof_symbols);
vec_sc_prod_cfc(y[1], 1.0/sqrtf(2), y[1], 2*nof_symbols);
return 2 * i;
} else if (nof_ports == 4) {
//int m_ap = (nof_symbols%4)?(nof_symbols*4-2):nof_symbols*4;

@ -97,47 +97,46 @@ int main(int argc, char **argv) {
}
for (i = 0; i < nof_layers; i++) {
x[i] = malloc(sizeof(cf_t) * nof_symbols);
x[i] = vec_malloc(sizeof(cf_t) * nof_symbols);
if (!x[i]) {
perror("malloc");
perror("vec_malloc");
exit(-1);
}
xr[i] = malloc(sizeof(cf_t) * nof_symbols);
xr[i] = calloc(1,sizeof(cf_t) * nof_symbols);
if (!xr[i]) {
perror("malloc");
perror("vec_malloc");
exit(-1);
}
}
for (i = 0; i < nof_ports; i++) {
y[i] = malloc(sizeof(cf_t) * nof_symbols * nof_layers);
y[i] = vec_malloc(sizeof(cf_t) * nof_symbols * nof_layers);
// TODO: The number of symbols per port is different in spatial multiplexing.
if (!y[i]) {
perror("malloc");
perror("vec_malloc");
exit(-1);
}
h[i] = malloc(sizeof(cf_t) * nof_symbols * nof_layers);
h[i] = vec_malloc(sizeof(cf_t) * nof_symbols * nof_layers);
if (!h[i]) {
perror("malloc");
perror("vec_malloc");
exit(-1);
}
}
/* only 1 receiver antenna supported now */
r[0] = malloc(sizeof(cf_t) * nof_symbols * nof_layers);
r[0] = vec_malloc(sizeof(cf_t) * nof_symbols * nof_layers);
if (!r[0]) {
perror("malloc");
perror("vec_malloc");
exit(-1);
}
/* generate random data */
for (i = 0; i < nof_layers; i++) {
for (j = 0; j < nof_symbols; j++) {
x[i][j] = 100
* ((float) rand() / RAND_MAX + (float) I * rand() / RAND_MAX);
x[i][j] = (float) rand()/RAND_MAX+((float) rand()/RAND_MAX)*_Complex_I;
}
}
if (precoding_init(&precoding, nof_symbols)) {
if (precoding_init(&precoding, nof_symbols * nof_layers)) {
fprintf(stderr, "Error initializing precoding\n");
exit(-1);
}
@ -150,9 +149,12 @@ int main(int argc, char **argv) {
/* generate channel */
for (i = 0; i < nof_ports; i++) {
for (j = 0; j < nof_symbols * nof_layers; j++) {
float hc = -1 + (float) i / nof_ports;
h[i][j] = (3 + hc) * cexpf(I * hc);
for (j = 0; j < nof_symbols; j++) {
h[i][nof_layers*j] = (float) rand()/RAND_MAX+((float) rand()/RAND_MAX)*_Complex_I;
// assume the channel is time-invariant in nlayer consecutive symbols
for (int k=0;k<nof_layers;k++) {
h[i][nof_layers*j+k] = h[i][nof_layers*j];
}
}
}
@ -168,18 +170,32 @@ int main(int argc, char **argv) {
}
}
/* predecoding / equalization */
if (predecoding_type(&precoding, r[0], h, xr, nof_ports, nof_layers,
nof_symbols * nof_layers, type) < 0) {
fprintf(stderr, "Error layer mapper encoder\n");
exit(-1);
}
float noise_estimate[2] = {0, 0.0000001};
const char *algs[2] = {"ZF", "MMSE"};
/* check errors */
mse = 0;
for (i = 0; i < nof_layers; i++) {
for (j = 0; j < nof_symbols; j++) {
mse += cabsf(xr[i][j] - x[i][j]) / nof_layers / nof_symbols;
/* predecoding / equalization */
struct timeval t[3];
for (int a=0;a<2;a++) {
gettimeofday(&t[1], NULL);
if (predecoding_type(&precoding, r[0], h, xr, nof_ports, nof_layers,
nof_symbols * nof_layers, type, noise_estimate[a]) < 0) {
fprintf(stderr, "Error layer mapper encoder\n");
exit(-1);
}
gettimeofday(&t[2], NULL);
get_time_interval(t);
printf("Execution Time %s: %d us\n", algs[a], t[0].tv_usec);
/* check errors */
mse = 0;
for (i = 0; i < nof_layers; i++) {
for (j = 0; j < nof_symbols; j++) {
mse += cabsf(xr[i][j] - x[i][j]);
}
}
printf("MSE %s: %f\n", algs[a], mse);
if (mse / nof_layers / nof_symbols > MSE_THRESHOLD) {
exit(-1);
}
}
@ -196,11 +212,6 @@ int main(int argc, char **argv) {
precoding_free(&precoding);
if (mse > MSE_THRESHOLD) {
printf("MSE: %f\n", mse);
exit(-1);
} else {
printf("Ok\n");
exit(0);
}
printf("Ok\n");
exit(0);
}

@ -377,11 +377,11 @@ int pbch_decode(pbch_t *q, cf_t *slot1_symbols, cf_t *ce_slot1[MAX_PORTS], float
/* in conctrol channels, only diversity is supported */
if (nant == 1) {
/* no need for layer demapping */
predecoding_single_mmse(&q->precoding, q->pbch_symbols[0], q->ce[0], q->pbch_d,
predecoding_single(&q->precoding, q->pbch_symbols[0], q->ce[0], q->pbch_d,
q->nof_symbols, noise_estimate);
} else {
predecoding_diversity_zf(&q->precoding, q->pbch_symbols[0], q->ce, x, nant,
q->nof_symbols);
predecoding_diversity(&q->precoding, q->pbch_symbols[0], q->ce, x, nant,
q->nof_symbols, noise_estimate);
layerdemap_diversity(x, q->pbch_d, nant, q->nof_symbols / nant);
}

@ -190,11 +190,11 @@ int pcfich_decode(pcfich_t *q, cf_t *slot_symbols, cf_t *ce[MAX_PORTS], float no
/* in control channels, only diversity is supported */
if (q->cell.nof_ports == 1) {
/* no need for layer demapping */
predecoding_single_mmse(&q->precoding, q->pcfich_symbols[0], q->ce[0], q->pcfich_d,
predecoding_single(&q->precoding, q->pcfich_symbols[0], q->ce[0], q->pcfich_d,
q->nof_symbols, noise_estimate);
} else {
predecoding_diversity_zf(&q->precoding, q->pcfich_symbols[0], ce_precoding, x,
q->cell.nof_ports, q->nof_symbols);
predecoding_diversity(&q->precoding, q->pcfich_symbols[0], ce_precoding, x,
q->cell.nof_ports, q->nof_symbols, noise_estimate);
layerdemap_diversity(x, q->pcfich_d, q->cell.nof_ports,
q->nof_symbols / q->cell.nof_ports);
}

@ -402,9 +402,9 @@ int pdcch_extract_llr(pdcch_t *q, cf_t *sf_symbols, cf_t *ce[MAX_PORTS], float n
/* in control channels, only diversity is supported */
if (q->cell.nof_ports == 1) {
/* no need for layer demapping */
predecoding_single_mmse(&q->precoding, q->pdcch_symbols[0], q->ce[0], q->pdcch_d, nof_symbols, noise_estimate);
predecoding_single(&q->precoding, q->pdcch_symbols[0], q->ce[0], q->pdcch_d, nof_symbols, noise_estimate);
} else {
predecoding_diversity_zf(&q->precoding, q->pdcch_symbols[0], q->ce, x, q->cell.nof_ports, nof_symbols);
predecoding_diversity(&q->precoding, q->pdcch_symbols[0], q->ce, x, q->cell.nof_ports, nof_symbols, noise_estimate);
layerdemap_diversity(x, q->pdcch_d, q->cell.nof_ports, nof_symbols / q->cell.nof_ports);
}

@ -664,11 +664,11 @@ int pdsch_decode(pdsch_t *q, cf_t *sf_symbols, cf_t *ce[MAX_PORTS], float noise_
/* TODO: only diversity is supported */
if (q->cell.nof_ports == 1) {
/* no need for layer demapping */
predecoding_single_mmse(&q->precoding, q->pdsch_symbols[0], q->ce[0], q->pdsch_d,
predecoding_single(&q->precoding, q->pdsch_symbols[0], q->ce[0], q->pdsch_d,
nof_symbols, noise_estimate);
} else {
predecoding_diversity_zf(&q->precoding, q->pdsch_symbols[0], q->ce, x, q->cell.nof_ports,
nof_symbols);
predecoding_diversity(&q->precoding, q->pdsch_symbols[0], q->ce, x, q->cell.nof_ports,
nof_symbols, noise_estimate);
layerdemap_diversity(x, q->pdsch_d, q->cell.nof_ports,
nof_symbols / q->cell.nof_ports);
}

@ -208,11 +208,11 @@ int phich_decode(phich_t *q, cf_t *slot_symbols, cf_t *ce[MAX_PORTS], float nois
/* in control channels, only diversity is supported */
if (q->cell.nof_ports == 1) {
/* no need for layer demapping */
predecoding_single_mmse(&q->precoding, q->phich_symbols[0], q->ce[0], q->phich_d0,
predecoding_single(&q->precoding, q->phich_symbols[0], q->ce[0], q->phich_d0,
PHICH_MAX_NSYMB, noise_estimate);
} else {
predecoding_diversity_zf(&q->precoding, q->phich_symbols[0], ce_precoding, x,
q->cell.nof_ports, PHICH_MAX_NSYMB);
predecoding_diversity(&q->precoding, q->phich_symbols[0], ce_precoding, x,
q->cell.nof_ports, PHICH_MAX_NSYMB, noise_estimate);
layerdemap_diversity(x, q->phich_d0, q->cell.nof_ports,
PHICH_MAX_NSYMB / q->cell.nof_ports);
}

@ -377,9 +377,21 @@ void vec_div_ccc(cf_t *x, cf_t *y, float *y_mod, cf_t *z, float *z_real, float *
#ifdef DIV_USE_VEC
vec_prod_conj_ccc(x,y,z,len);
vec_abs_square_cf(y,y_mod,len);
vec_deinterleave_cf(z, z_real, z_imag, len);
vec_div_fff(z_real, y_mod, z_real, len);
vec_div_fff(z_imag, y_mod, z_imag, len);
vec_div_cfc(z,y_mod,z,z_real,z_imag,len);
#else
int i;
for (i=0;i<len;i++) {
z[i] = x[i] / y[i];
}
#endif
}
/* Complex division by float z=x/y */
void vec_div_cfc(cf_t *x, float *y, cf_t *z, float *z_real, float *z_imag, uint32_t len) {
#ifdef DIV_USE_VEC
vec_deinterleave_cf(x, z_real, z_imag, len);
vec_div_fff(z_real, y, z_real, len);
vec_div_fff(z_imag, y, z_imag, len);
vec_interleave_cf(z_real, z_imag, z, len);
#else
int i;

@ -0,0 +1,285 @@
%% PDSCH Transmit Diversity Throughput Conformance Test
% This example demonstrates the throughput performance under conformance
% test conditions as defined in TS36.101[ <#9 1> ]: single codeword,
% transmit diversity 4Tx-2Rx with medium correlation, EPA5 (Extended
% Pedestrian A) channel. The example also introduces the use of Parallel
% Computing Toolbox(TM) to provide improvements in the simulation time.
% Copyright 2009-2013 The MathWorks, Inc.
%% Introduction
% In this example, Hybrid Automatic Repeat Request(HARQ) is used in line
% with conformance test requirements. A total of 8 HARQ processes are used
% with a maximum of 4 retransmissions permitted. This example uses the R.12
% Reference Measurement Channel (RMC).
%
% This example also uses <matlab:doc('parfor') parfor> loop instead of the
% <matlab:doc('for') for> loop for SNR calculation. <matlab:doc('parfor')
% parfor>, as part of the Parallel Computing Toolbox, executes the SNR loop
% in parallel to reduce the total simulation time.
%% Simulation Settings
% The default simulation length is set to 10 frames at a number of |SNR|
% values including 0.2dB (as per TS36.101 (Section 8.2.1.2.2-Test1)[ <#9
% 1> ]).
NFrames = 50; % Number of frames
SNRdB = [4 5 6]; % SNR range
% eNodeB Configuration
enb = struct; % eNodeB config structure
enb.TotSubframes = 1; % Total subframes RMC will generate
enb.RC = 'R.10'; % RMC number
% Channel Configuration
channel = struct; % Channel config structure
channel.Seed = 2; % Random channel seed
channel.NRxAnts = 1; % 2 receive antennas
channel.DelayProfile ='EVA'; % Delay profile
channel.DopplerFreq = 5; % Doppler frequency
channel.MIMOCorrelation = 'Low'; % Multi-antenna correlation
channel.NTerms = 16; % Oscillators used in fading model
channel.ModelType = 'GMEDS'; % Rayleigh fading model type
channel.InitPhase = 'Random'; % Random initial phases
channel.NormalizePathGains = 'On'; % Normalize delay profile power
channel.NormalizeTxAnts = 'On'; % Normalize for transmit antennas
% Channel Estimator Configuration
cec = struct; % Channel estimation config structure
cec.PilotAverage = 'UserDefined'; % Type of pilot symbol averaging
cec.FreqWindow = 9; % Frequency window size
cec.TimeWindow = 9; % Time window size
cec.InterpType = 'Linear'; % 2D interpolation type
cec.InterpWindow = 'Centered'; % Interpolation window type
cec.InterpWinSize = 1; % Interpolation window size
% PDSCH Configuration
enb.PDSCH.TxScheme = 'TxDiversity'; % Transmission scheme
enb.PDSCH.RNTI = 1; % 16-bit User Equipment (UE) mask
enb.PDSCH.Rho = 0; % PDSCH RE power adjustment factor
enb.PDSCH.CSI = 'Off'; % No CSI scaling of soft bits
% Simulation Variables
totalBLKCRC = []; % Define total block CRC vector
bitThroughput = []; % Define total bit throughput vector
%% System Processing
% Working on a subframe by subframe basis and using the LTE System
% Toolbox(TM) a populated resource grid is generated and OFDM modulated to
% create a transmit waveform. The generated waveform is transmitted through
% a propagation channel and AWGN is added. Channel estimation, equalization
% and the inverse of transmission chain are performed at receiver. The
% throughput performance of the PDSCH is determined using the block CRC
% result.
% Generate the RMC configuration structure for RMC R.12
rmc = lteRMCDL(enb);
rvSeq = rmc.PDSCH.RVSeq;
% Transport block sizes for each subframe in a frame
trBlkSizes = rmc.PDSCH.TrBlkSizes;
codedTrBlkSizes = rmc.PDSCH.CodedTrBlkSizes;
% Determine resource grid dimensions
dims = lteDLResourceGridSize(rmc);
p = dims(3);
% Set up channel model sampling rate
ofdmInfo = lteOFDMInfo(rmc);
channel.SamplingRate = ofdmInfo.SamplingRate;
% Generation HARQ table for 8-HARQ processes
harqTable = hHARQTable();
% Initializing state of all HARQ processes
for i=1:9
harqProcess_init(i) = hTxDiversityNewHARQProcess ...
(trBlkSizes(i),codedTrBlkSizes(i),rvSeq); %#ok<SAGROW>
end
% Display the SNR points being simulated
for s=1:numel(SNRdB)
fprintf('\nSimulating at %gdB SNR for a total %d Frame(s)\n', ...
SNRdB(s),NFrames);
end
% The temporary variables 'rmc_init' and 'channel_init' are used to create
% the temporary variables 'rmc' and 'channel' within the SNR loop to create
% independent simulation loops for the parfor loop
rmc_init = rmc;
channel_init = channel;
% 'parfor' will default to the normal 'for' when executed without Parallel
% Computing Toolbox.
parfor index = 1:numel(SNRdB)
% Set the random number generator seed depending to the loop variable
% to ensure independent random streams
rng(index,'combRecursive');
% Set up variables for the SNR loop
offsets = 0; % Initialize overall frame offset value for the SNR
offset = 0; % Initialize frame offset value for the radio frame
rmc = rmc_init; % Initialize RMC configuration
channel = channel_init; % Initialize channel configuration
blkCRC = []; % Define intermediate block CRC vector
bitTput = []; % Intermediate bit throughput vector
% Initializing state of all HARQ processes
harqProcess = harqProcess_init;
for subframeNo = 0:(NFrames*10-1)
% Updating subframe number
rmc.NSubframe = subframeNo;
% HARQ index table
harqIdx = harqTable(mod(subframeNo,length(harqTable))+1); %#ok<PFBNS>
% Update HARQ process
harqProcess(harqIdx) = hTxDiversityHARQScheduling( ...
harqProcess(harqIdx));
% Updating the RV value for correct waveform generation
rmc.PDSCH.RV = harqProcess(harqIdx).rvSeq ...
(harqProcess(harqIdx).rvIdx);
rmc.PDSCH.RVSeq = harqProcess(harqIdx).rvSeq ...
(harqProcess(harqIdx).rvIdx);
[txWaveform txGrid] = lteRMCDLTool(rmc, ...
{harqProcess(harqIdx).dlschTransportBlk});
txWaveform = [txWaveform; zeros(25,p)];
% Initialize at time zero
channel.InitTime = subframeNo/1000;
% Pass data through the fading channel model
rxWaveform = lteFadingChannel(channel,txWaveform);
% Noise setup including compensation for downlink power allocation
SNR = 10^((SNRdB(index)-rmc.PDSCH.Rho)/20); % Linear SNR
% Normalize noise power to take account of sampling rate, which is
% a function of the IFFT size used in OFDM modulation, and the
% number of antennas
N0 = 1/(sqrt(2.0*rmc.CellRefP*double(ofdmInfo.Nfft))*SNR);
% Create additive white Gaussian noise
noise = N0*complex(randn(size(rxWaveform)), ...
randn(size(rxWaveform)));
% Add AWGN to the received time domain waveform
rxWaveform = rxWaveform + noise;
% Receiver
% Perform synchronization
% An offset within the range of delays expected from the channel
% modeling(a combination of implementation delay and channel delay
% spread) indicates success
if (mod(subframeNo,10)==0)
[offset] = lteDLFrameOffset(rmc,rxWaveform);
if (offset > 25)
offset = offsets(end);
end
offsets = [offsets offset];
end
rxWaveform = rxWaveform(1+offset:end,:);
% Perform OFDM demodulation on the received data to recreate the
% resource grid
rxSubframe = lteOFDMDemodulate(rmc,rxWaveform);
% Equalization and channel estimation
[estChannelGrid,noiseEst] = lteDLChannelEstimate(rmc,cec, ...
rxSubframe);
addpath('../../debug/lte/phy/lib/ch_estimation/test')
[est, ~, output] = liblte_chest(rmc.NCellID,rmc.CellRefP,rxSubframe,[0.25 0.5 0.25],[0.1 0.9],mod(rmc.NSubframe,10));
%estChannelGrid=reshape(est,size(estChannelGrid));
% Perform deprecoding, layer demapping, demodulation and
% descrambling on the received data using the estimate of
% the channel
rxEncodedBits = ltePDSCHDecode2(rmc,rmc.PDSCH,rxSubframe,estChannelGrid,noiseEst);
% Decode DownLink Shared Channel (DL-SCH)
[decbits,harqProcess(harqIdx).crc,harqProcess(harqIdx).decState] = ...
lteDLSCHDecode(rmc,rmc.PDSCH,harqProcess(harqIdx).trBlkSize, ...
rxEncodedBits{1},harqProcess(harqIdx).decState);
if(harqProcess(harqIdx).trBlkSize ~= 0)
blkCRC = [blkCRC harqProcess(harqIdx).crc];
bitTput = [bitTput harqProcess(harqIdx).trBlkSize.*(1- ...
harqProcess(harqIdx).crc)];
end
end
% Record the block CRC and bit throughput for the total number of
% frames simulated at a particular SNR
totalBLKCRC(index,:) = blkCRC;
bitThroughput(index,:) = bitTput;
end
%%
% |totalBLKCRC| is a matrix where each row contains the results of decoding
% the block CRC for a defined value of SNR. |bitThroughput| is a matrix
% containing the total number of bits per subframe at the different SNR
% points that have been successfully received and decoded.
%% Results
% First graph shows the throughput as total bits per second against the
% range of SNRs
% figure;
% plot(SNRdB,mean(bitThroughput,2),'-*');
% %axis([-5 3 200 400])
% title(['Throughput for ', num2str(NFrames) ' frame(s)'] );
% xlabel('SNRdB'); ylabel('Throughput (kbps)');
% grid on;
% hold on;
% plot(SNRdB,mean([trBlkSizes(1:5) trBlkSizes(7:10)])*0.7*ones ...
% (1,numel(SNRdB)),'--rs');
% legend('Simulation Result','70 Percent Throughput','Location','SouthEast');
%
% % Second graph shows the total throughput as a percentage of CRC passes
% % against SNR range
% figure;
plot(SNRdB,100*(1-mean(totalBLKCRC,2)),'-*');
%axis([-5 3 50 110])
title(['Throughput for ', num2str(NFrames) ' frame(s)'] );
xlabel('SNRdB'); ylabel('Throughput (%)');
grid on;
hold on;
plot(SNRdB,70*ones(1,numel(SNRdB)),'--rs');
legend('Simulation Result','70 Percent Throughput','Location','SouthEast');
%% Further Exploration
%
% You can modify parts of this example to experiment with different number
% of |NFrames| and different values of SNR. SNR can be a vector of
% values or a single value. Following scenarios can be simulated.
%%
% * Allows control over the total number of frames to run the demo at an
% SNR of 0.2dB (as per TS 36.101).
%
% * Allows control over the total number of frames to run the demo, as well
% as defining a set of desired SNR values. |SNRIn| can be a single value
% or a vector containing a range of values.
%
% * For simulations of multiple SNR points over a large number of frames,
% the use of Parallel Computing Toolbox provides significant improvement in
% the simulation time. This can be easily verified by changing the |parfor|
% in the SNR loop to |for| and re-running the example.
%% Appendix
% This example uses the following helper functions:
%
% * <matlab:edit('hHARQTable.m') hHARQTable.m>
% * <matlab:edit('hTxDiversityHARQScheduling.m') hTxDiversityHARQScheduling.m>
% * <matlab:edit('hTxDiversityNewHARQProcess.m') hTxDiversityNewHARQProcess.m>
%% Selected Bibliography
% # 3GPP TS 36.101
displayEndOfDemoMessage(mfilename)

@ -4,7 +4,7 @@
clear
SNR_values_db=linspace(0,20,8);
SNR_values_db=[0 1 2];%15;%[5 10 15];%linspace(0,20,8);
Nrealizations=1;
preEVM = zeros(length(SNR_values_db),Nrealizations);
@ -22,10 +22,10 @@ enb.DuplexMode = 'FDD'; % FDD
%% Channel Model Configuration
rng(1); % Configure random number generators
cfg.Seed = 1; % Random channel seed
cfg.Seed = 2; % Random channel seed
cfg.NRxAnts = 1; % 1 receive antenna
cfg.DelayProfile = 'EVA'; % EVA delay spread
cfg.DopplerFreq = 120; % 120Hz Doppler frequency
cfg.DopplerFreq = 5; % 120Hz Doppler frequency
cfg.MIMOCorrelation = 'Low'; % Low (no) MIMO correlation
cfg.InitTime = 0; % Initialize at time zero
cfg.NTerms = 16; % Oscillators used in fading model
@ -35,24 +35,13 @@ cfg.NormalizePathGains = 'On'; % Normalize delay profile power
cfg.NormalizeTxAnts = 'On'; % Normalize for transmit antennas
%% Channel Estimator Configuration
cec.FreqWindow = 9; % Frequency averaging window in
% Resource Elements (REs)
cec.TimeWindow = 9; % Time averaging window in REs
cec.InterpType = 'Cubic'; % Cubic interpolation
cec.PilotAverage = 'UserDefined'; % Pilot averaging method
cec.InterpWinSize = 3; % Interpolate up to 3 subframes
% simultaneously
cec.InterpWindow = 'Centred'; % Interpolation windowing method
cec2.FreqWindow = 9; % Frequency averaging window in
% Resource Elements (REs)
cec2.TimeWindow = 9; % Time averaging window in REs
cec2.InterpType = 'Linear'; % Cubic interpolation
cec2.PilotAverage = 'UserDefined'; % Pilot averaging method
cec2.InterpWinSize = 3; % Interpolate up to 3 subframes
% simultaneously
cec2.InterpWindow = 'Causal'; % Interpolation windowing method
cec = struct; % Channel estimation config structure
cec.PilotAverage = 'UserDefined'; % Type of pilot symbol averaging
cec.FreqWindow = 9; % Frequency window size
cec.TimeWindow = 9; % Time window size
cec.InterpType = 'Linear'; % 2D interpolation type
cec.InterpWindow = 'Centered'; % Interpolation window type
cec.InterpWinSize = 1; % Interpolation window size
%% Subframe Resource Grid Size
@ -151,14 +140,16 @@ rxWaveform = rxWaveform(1+offset:end,:);
rxGrid = lteOFDMDemodulate(enb,rxWaveform);
addpath('../../debug/lte/phy/lib/ch_estimation/test')
%% Channel Estimation
[estChannel, noiseEst] = lteDLChannelEstimate(enb,cec,rxGrid);
[estChannel_lin, noiseEst_lin] = lteDLChannelEstimate(enb,cec2,rxGrid);
[d, a, output] = liblte_chest(enb.NCellID,enb.CellRefP,rxGrid,[0.15 0.7 0.15],[0.1 0.9]);
output=[];
for i=0:9
[d, a, out] = liblte_chest(enb.NCellID,enb.CellRefP,rxGrid(:,i*14+1:(i+1)*14),[0.15 0.7 0.15],[0.1 0.9],i);
output = [output out];
end
%% MMSE Equalization
eqGrid_mmse = lteEqualizeMMSE(rxGrid, estChannel, noiseEst);
eqGrid_mmse_lin = lteEqualizeMMSE(rxGrid, estChannel_lin, noiseEst_lin);
eqGrid_liblte = reshape(output,size(eqGrid_mmse));
@ -170,24 +161,19 @@ preEqualisedEVM = lteEVM(txGrid,rxGrid);
fprintf('%d-%d: Pre-EQ: %0.3f%%\n', ...
snr_idx,nreal,preEqualisedEVM.RMS*100);
%EVM of post-equalized receive signal
postEqualisedEVM_mmse = lteEVM(txGrid,eqGrid_mmse);
postEqualisedEVM_mmse = lteEVM(txGrid,reshape(eqGrid_mmse,size(txGrid)));
fprintf('%d-%d: MMSE: %0.3f%%\n', ...
snr_idx,nreal,postEqualisedEVM_mmse.RMS*100);
postEqualisedEVM_mmse_lin = lteEVM(txGrid,eqGrid_mmse_lin);
fprintf('%d-%d: MMSE-LIN: %0.3f%%\n', ...
snr_idx,nreal,postEqualisedEVM_mmse_lin.RMS*100);
postEqualisedEVM_liblte = lteEVM(txGrid,eqGrid_liblte);
postEqualisedEVM_liblte = lteEVM(txGrid,reshape(eqGrid_liblte,size(txGrid)));
fprintf('%d-%d: liblte: %0.3f%%\n', ...
snr_idx,nreal,postEqualisedEVM_liblte.RMS*100);
preEVM(snr_idx,nreal) =preEqualisedEVM.RMS;
postEVM_mmse(snr_idx,nreal) = postEqualisedEVM_mmse.RMS;
postEVM_mmse_lin(snr_idx,nreal) = postEqualisedEVM_mmse_lin.RMS;
postEVM_liblte(snr_idx,nreal) = postEqualisedEVM_liblte.RMS;
preEVM(snr_idx,nreal) = preEqualisedEVM.RMS;
postEVM_mmse(snr_idx,nreal) = mean([postEqualisedEVM_mmse.RMS]);
postEVM_liblte(snr_idx,nreal) = mean([postEqualisedEVM_liblte.RMS]);
end
end
@ -195,7 +181,6 @@ end
% legend('real','seu','meu')
plot(SNR_values_db, mean(preEVM,2), ...
SNR_values_db, mean(postEVM_mmse,2), ...
SNR_values_db, mean(postEVM_mmse_lin,2), ...
SNR_values_db, mean(postEVM_liblte,2))
legend('No Eq','MMSE-cubic','MMSE-lin','MMSE-liblte')
grid on

Loading…
Cancel
Save