diff --git a/lte/examples/cell_search_utils.c b/lte/examples/cell_search_utils.c
index e9576276c..331b4a6ca 100644
--- a/lte/examples/cell_search_utils.c
+++ b/lte/examples/cell_search_utils.c
@@ -125,11 +125,7 @@ int find_cell(void *uhd, ue_celldetect_result_t *found_cell, uint32_t N_id_2)
     fprintf(stderr, "Error initiating UE cell detect\n");
     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);
@@ -222,6 +218,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);
diff --git a/lte/phy/include/liblte/phy/mimo/precoding.h b/lte/phy/include/liblte/phy/mimo/precoding.h
index adecc03b8..d0e5abf09 100644
--- a/lte/phy/include/liblte/phy/mimo/precoding.h
+++ b/lte/phy/include/liblte/phy/mimo/precoding.h
@@ -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_ */
diff --git a/lte/phy/include/liblte/phy/utils/vector.h b/lte/phy/include/liblte/phy/utils/vector.h
index 79bdaa819..eaa8150b3 100644
--- a/lte/phy/include/liblte/phy/utils/vector.h
+++ b/lte/phy/include/liblte/phy/utils/vector.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 */
diff --git a/lte/phy/lib/ch_estimation/test/chest_test_dl.c b/lte/phy/lib/ch_estimation/test/chest_test_dl.c
index bf5e216e9..8b60a23c0 100644
--- a/lte/phy/lib/ch_estimation/test/chest_test_dl.c
+++ b/lte/phy/lib/ch_estimation/test/chest_test_dl.c
@@ -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);
diff --git a/lte/phy/lib/ch_estimation/test/chest_test_dl_mex.c b/lte/phy/lib/ch_estimation/test/chest_test_dl_mex.c
index 8c9e320a6..6d6a4768e 100644
--- a/lte/phy/lib/ch_estimation/test/chest_test_dl_mex.c
+++ b/lte/phy/lib/ch_estimation/test/chest_test_dl_mex.c
@@ -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,18 +193,24 @@ 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]);
-        }
-        outr0++;
-        outi0++;
-      } 
+      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++;
+        } 
+      }
     }
     if (nlhs >= 2) {    
       for (int j=0;j<cell.nof_ports;j++) {
@@ -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++;
diff --git a/lte/phy/lib/mimo/src/precoding.c b/lte/phy/lib/mimo/src/precoding.c
index e501a5d42..2dd39e34d 100644
--- a/lte/phy/lib/mimo/src/precoding.c
+++ b/lte/phy/lib/mimo/src/precoding.c
@@ -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;
diff --git a/lte/phy/lib/mimo/test/precoding_test.c b/lte/phy/lib/mimo/test/precoding_test.c
index 7efb3e0b1..5528fb8fb 100644
--- a/lte/phy/lib/mimo/test/precoding_test.c
+++ b/lte/phy/lib/mimo/test/precoding_test.c
@@ -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];
+      }
     }
   }
 
@@ -167,20 +169,34 @@ int main(int argc, char **argv) {
       r[0][j] += y[i][j] * h[i][j];
     }
   }
-
+  
+  float noise_estimate[2] = {0, 0.0000001};
+  const char *algs[2] = {"ZF", "MMSE"};
+  
   /* 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);
-  }
-
-  /* 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;
+  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);
+    } 
   }
 
   for (i = 0; i < nof_layers; i++) {
@@ -195,12 +211,7 @@ int main(int argc, char **argv) {
   free(r[0]);
   
   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); 
 }
diff --git a/lte/phy/lib/phch/src/pbch.c b/lte/phy/lib/phch/src/pbch.c
index bedbfb669..c32daf979 100644
--- a/lte/phy/lib/phch/src/pbch.c
+++ b/lte/phy/lib/phch/src/pbch.c
@@ -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);
         }
 
diff --git a/lte/phy/lib/phch/src/pcfich.c b/lte/phy/lib/phch/src/pcfich.c
index 8db9fbd60..5dd68ea05 100644
--- a/lte/phy/lib/phch/src/pcfich.c
+++ b/lte/phy/lib/phch/src/pcfich.c
@@ -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);
     }
diff --git a/lte/phy/lib/phch/src/pdcch.c b/lte/phy/lib/phch/src/pdcch.c
index d39a282d5..9a5351f7d 100644
--- a/lte/phy/lib/phch/src/pdcch.c
+++ b/lte/phy/lib/phch/src/pdcch.c
@@ -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);
       }
 
diff --git a/lte/phy/lib/phch/src/pdsch.c b/lte/phy/lib/phch/src/pdsch.c
index 4d4b91d2f..a9d223cf0 100644
--- a/lte/phy/lib/phch/src/pdsch.c
+++ b/lte/phy/lib/phch/src/pdsch.c
@@ -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);
     }
diff --git a/lte/phy/lib/phch/src/phich.c b/lte/phy/lib/phch/src/phich.c
index 8bb45a695..d039a7de9 100644
--- a/lte/phy/lib/phch/src/phich.c
+++ b/lte/phy/lib/phch/src/phich.c
@@ -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);
   }
diff --git a/lte/phy/lib/utils/src/vector.c b/lte/phy/lib/utils/src/vector.c
index c264ce5cb..a807f10bf 100644
--- a/lte/phy/lib/utils/src/vector.c
+++ b/lte/phy/lib/utils/src/vector.c
@@ -377,10 +377,7 @@ 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_interleave_cf(z_real, z_imag, z, len);
+  vec_div_cfc(z,y_mod,z,z_real,z_imag,len);  
 #else 
   int i; 
   for (i=0;i<len;i++) {
@@ -389,6 +386,21 @@ void vec_div_ccc(cf_t *x, cf_t *y, float *y_mod, cf_t *z, float *z_real, float *
 #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; 
+  for (i=0;i<len;i++) {
+    z[i] = x[i] / y[i]; 
+  }
+#endif
+}
+
 void vec_div_fff(float *x, float *y, float *z, uint32_t len) {
 #ifdef HAVE_VOLK_DIVIDE_FUNCTION
   volk_32f_x2_divide_32f(z, x, y, len);
diff --git a/matlab/tests/diversity_test_rate.m b/matlab/tests/diversity_test_rate.m
new file mode 100644
index 000000000..0cefaca5a
--- /dev/null
+++ b/matlab/tests/diversity_test_rate.m
@@ -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) 
diff --git a/matlab/tests/equalizer_test.m b/matlab/tests/equalizer_test.m
index 7c16a86c5..573493d7b 100644
--- a/matlab/tests/equalizer_test.m
+++ b/matlab/tests/equalizer_test.m
@@ -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));
 
@@ -169,25 +160,20 @@ eqGrid_liblte = reshape(output,size(eqGrid_mmse));
 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