1. Adding cblast and lapacke-based functions required for mmse computation

to linear_preprocessing_rec.c.
2. Adding linear_preprocessing_rec.c to PHY_SRC_UE library.
3. Adding mmse_flag and mmse functionalities to dlsch_demodulation.
4. For now, dlsim_tm4 will not compile.

cmake_targets/CMakeLists.txt

@@ -1177,6 +1177,7 @@ set(PHY_SRC_UE
   ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/sss_ue.c
   ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/dlsch_demodulation.c
   ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/dlsch_llr_computation.c
+  ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/linear_preprocessing_rec.c
   ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/dlsch_decoding.c
   ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/dci_tools_ue.c
   ${OPENAIR1_DIR}/PHY/LTE_UE_TRANSPORT/uci_tools_ue.c
@@ -1817,11 +1818,13 @@ endif()
 # So, here are some hacks here. Hope this gets fixed in future!
 if(EXISTS "/usr/include/atlas/cblas.h" OR EXISTS "/usr/include/cblas.h")
   include_directories("/usr/include/atlas")
+  LINK_DIRECTORIES("/usr/lib/lapack")
   LINK_DIRECTORIES("/usr/lib64")
   LINK_DIRECTORIES("/usr/lib64/atlas") #Added because atlas libraries in CentOS 7 are here!
   
   if(EXISTS "/usr/lib64/libblas.so" OR EXISTS "/usr/lib/libblas.so") #Case for CentOS7
      list(APPEND ATLAS_LIBRARIES blas)
+
   else() # Case for Ubuntu
      list(APPEND ATLAS_LIBRARIES cblas)
   endif()
@@ -1847,6 +1850,8 @@ else()
   message("No Blas/Atlas libs found, some targets will fail")
 endif()
 
+list(APPEND ATLAS_LIBRARIES lapack lapacke)
+
 if (${XFORMS})
   include_directories ("/usr/include/X11")
   set(XFORMS_SOURCE

openair1/PHY/LTE_UE_TRANSPORT/linear_preprocessing_rec.c

+/* These functions compute linear preprocessing for
+the UE using LAPACKE and CBLAS modules of
+LAPACK libraries.
+MMSE and MMSE whitening filters are available.
+Functions are using RowMajor storage of the
+matrices, like in conventional C. Traditional
+Fortran functions of LAPACK employ ColumnMajor
+data storage. */
+
+#include<stdio.h>
+#include<math.h>
+#include<complex.h>
+#include <stdlib.h>
+#include <cblas.h>
+#include <string.h>
+#include <lapacke_utils.h>
+#include <lapacke.h>
+
+//#define DEBUG_PREPROC
+
+
+void transpose (int N, float complex *A, float complex *Result)
+{
+  // COnputes C := alpha*op(A)*op(B) + beta*C,
+  enum CBLAS_TRANSPOSE transa = CblasTrans;
+  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
+  int rows_opA = N; // number of rows in op(A) and in C
+  int col_opB = N; //number of columns of op(B) and in C
+  int col_opA = N; //number of columns in op(A) and rows in op(B)
+  int col_B; //number of columns in B
+  float complex alpha = 1.0+I*0;
+  int lda  = rows_opA;
+  float complex beta = 0.0+I*0;
+  int ldc = rows_opA;
+  int i;
+  float complex* B;
+
+  int ldb = col_opB;
+
+  if (transb == CblasNoTrans) {
+    B = (float complex*)calloc(ldb*col_opB,sizeof(float complex));
+    col_B= col_opB;
+  }
+  else {
+    B = (float complex*)calloc(ldb*col_opA, sizeof(float complex));
+    col_B = col_opA;
+  }
+  float complex* C = (float complex*)malloc(ldc*col_opB*sizeof(float complex));
+
+  for (i=0; i<lda*col_B; i+=N+1)
+    B[i]=1.0+I*0;
+
+  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+  memcpy(Result, C, N*N*sizeof(float complex));
+
+  free(B);
+  free(C);
+ }
+
+
+void conjugate_transpose (int N, float complex *A, float complex *Result)
+{
+  // Computes C := alpha*op(A)*op(B) + beta*C,
+  enum CBLAS_TRANSPOSE transa = CblasConjTrans;
+  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
+  int rows_opA = N; // number of rows in op(A) and in C
+  int col_opB = N; //number of columns of op(B) and in C
+  int col_opA = N; //number of columns in op(A) and rows in op(B)
+  int col_B; //number of columns in B
+  float complex alpha = 1.0+I*0;
+  int lda  = rows_opA;
+  float complex beta = 0.0+I*0;
+  int ldc = rows_opA;
+  int i;
+  float complex* B;
+  int ldb = col_opB;
+
+  if (transb == CblasNoTrans) {
+    B = (float complex*)calloc(ldb*col_opB,sizeof(float complex));
+    col_B= col_opB;
+  }
+  else {
+    B = (float complex*)calloc(ldb*col_opA, sizeof(float complex));
+    col_B = col_opA;
+  }
+  float complex* C = (float complex*)malloc(ldc*col_opB*sizeof(float complex));
+
+  for (i=0; i<lda*col_B; i+=N+1)
+    B[i]=1.0+I*0;
+
+  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+  memcpy(Result, C, N*N*sizeof(float complex));
+
+  free(B);
+  free(C);
+ }
+
+void H_hermH_plus_sigma2I (int N, int M, float complex *A, float sigma2, float complex *Result)
+{
+  //C := alpha*op(A)*op(B) + beta*C,
+  enum CBLAS_TRANSPOSE transa = CblasConjTrans;
+  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
+  int rows_opA = N; // number of rows in op(A) and in C
+  int col_opB = N; //number of columns of op(B) and in C
+  int col_opA = N; //number of columns in op(A) and rows in op(B)
+  int col_C = N; //number of columns in B
+  float complex alpha = 1.0+I*0;
+  int lda  = col_opA;
+  float complex beta = 1.0 + I*0;
+  int ldc = col_opA;
+  int i;
+
+  float complex* C = (float complex*)calloc(ldc*col_opB, sizeof(float complex));
+
+  for (i=0; i<lda*col_C; i+=N+1)
+    C[i]=sigma2*(1.0+I*0);
+
+  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, A, lda, &beta, C, ldc);
+
+  memcpy(Result, C, N*M*sizeof(float complex));
+  free(C);
+
+ }
+
+ void HH_herm_plus_sigma2I (int M, int N, float complex *A, float sigma2, float complex *Result)
+{
+  //C := alpha*op(A)*op(B) + beta*C,
+  enum CBLAS_TRANSPOSE transa = CblasNoTrans;
+  enum CBLAS_TRANSPOSE transb = CblasConjTrans;
+  int k = N; //number of columns in op(A) and rows in op(B),k
+  float complex alpha = 1.0+I*0;
+  int lda  = N;
+  int ldb  = N;
+  int ldc = M;
+  int i;
+
+  float complex* C = (float complex*)calloc(M*M, sizeof(float complex));
+
+  for (i=0; i<M*M; i+=M+1)
+    C[i]=1.0+I*0;
+
+  cblas_cgemm(CblasRowMajor, transa, transb, M, M, k, &alpha, A, lda, A, ldb, &sigma2, C, ldc);
+
+  memcpy(Result, C, M*M*sizeof(float complex));
+  free(C);
+
+}
+
+void eigen_vectors_values (int N, float complex *A, float complex *Vectors, float *Values_Matrix)
+{
+  // This function computes ORTHONORMAL eigenvectors and eigenvalues of matrix A,
+  // where Values_Matrix is a diagonal matrix of eigenvalues.
+  // A=Vectors*Values_Matrix*Vectors'
+  char jobz = 'V';
+  char uplo = 'U';
+  int order_A = N;
+  int lda = N;
+  int i;
+  float* Values = (float*)malloc(sizeof(float)*1*N);
+
+  LAPACKE_cheev(LAPACK_ROW_MAJOR, jobz, uplo, order_A, A, lda, Values);
+
+  memcpy(Vectors, A, N*N*sizeof(float complex));
+
+  for (i=0; i<lda; i+=1)
+    Values_Matrix[i*(lda+1)]=Values[i];
+
+  free(Values);
+}
+
+ void lin_eq_solver (int N, float complex* A, float complex* B, float complex* Result)
+{
+  int n = N;
+  int lda = N;
+  int ldb = N;
+  int nrhs = N;
+
+  char transa = 'N';
+  int* IPIV = malloc(N*N*sizeof(int));
+
+  // Compute LU-factorization
+  LAPACKE_cgetrf(LAPACK_ROW_MAJOR, n, nrhs, A, lda, IPIV);
+
+  // Solve AX=B
+  LAPACKE_cgetrs(LAPACK_ROW_MAJOR, transa, n, nrhs, A, lda, IPIV, B, ldb);
+
+  // cgetrs( "N", N, 4, A, lda, IPIV, B, ldb, INFO )
+
+  memcpy(Result, B, N*N*sizeof(float complex));
+
+  free(IPIV);
+
+}
+
+void mutl_matrix_matrix_row_based(float complex* M0, float complex* M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex* Result ){
+  enum CBLAS_TRANSPOSE transa = CblasNoTrans;
+  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
+  int rows_opA = rows_M0; // number of rows in op(A) and in C
+  int col_opB = col_M1; //number of columns of op(B) and in C
+  int col_opA = col_M0; //number of columns in op(A) and rows in op(B)
+  float complex alpha =1.0;
+  int lda  = col_M0;
+  float complex beta = 0.0;
+  int ldc = col_M1;
+  int ldb = col_M1;
+
+#ifdef DEBUG_PREPROC
+  int i=0;
+  printf("rows_M0 %d, col_M0 %d, rows_M1 %d, col_M1 %d\n", rows_M0, col_M0, rows_M1, col_M1);
+
+  for(i=0; i<rows_M0*col_M0; ++i)
+    printf(" rows_opA = %d, col_opB = %d, W_MMSE[%d] = (%f + i%f)\n", rows_opA, col_opB,  i , creal(M0[i]), cimag(M0[i]));
+
+  for(i=0; i<rows_M1*col_M1; ++i)
+    printf(" M1[%d] = (%f + i%f)\n", i , creal(M1[i]), cimag(M1[i]));
+#endif
+
+  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, M0, lda, M1, ldb, &beta, Result, ldc);
+
+#ifdef DEBUG_PREPROC
+  for(i=0; i<rows_opA*col_opB; ++i)
+    printf(" result[%d] = (%f + i%f)\n", i , creal(Result[i]), cimag(Result[i]));
+#endif
+
+}
+void mutl_matrix_matrix_col_based(float complex* M0, float complex* M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex* Result ){
+  enum CBLAS_TRANSPOSE transa = CblasNoTrans;
+  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
+  int rows_opA = rows_M0; // number of rows in op(A) and in C
+  int col_opB = col_M1; //number of columns of op(B) and in C
+  int col_opA = col_M0; //number of columns in op(A) and rows in op(B)
+  float complex alpha =1.0;
+  int lda  = col_M0;
+  float complex beta = 0.0;
+  int ldc = rows_M1;
+  int ldb = rows_M1;
+
+#ifdef DEBUG_PREPROC
+  int i = 0;
+  printf("rows_M0 %d, col_M0 %d, rows_M1 %d, col_M1 %d\n", rows_M0, col_M0, rows_M1, col_M1);
+
+  for(i=0; i<rows_M0*col_M0; ++i)
+    printf(" rows_opA = %d, col_opB = %d, W_MMSE[%d] = (%f + i%f)\n", rows_opA, col_opB,  i , creal(M0[i]), cimag(M0[i]));
+
+
+  for(i=0; i<rows_M1*col_M1; ++i)
+    printf(" M1[%d] = (%f + i%f)\n", i , creal(M1[i]), cimag(M1[i]));
+#endif
+
+  cblas_cgemm(CblasColMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, M0, lda, M1, ldb, &beta, Result, ldc);
+
+#ifdef DEBUG_PREPROC
+  for(i=0; i<rows_opA*col_opB; ++i)
+    printf(" result[%d] = (%f + i%f)\n", i , creal(Result[i]), cimag(Result[i]));
+#endif
+}
+
+
+/*FILTERS */
+void compute_MMSE(float complex* H, int order_H, float sigma2, float complex* W_MMSE)
+{
+  int N = order_H;
+  float complex* H_hermH_sigmaI = malloc(N*N*sizeof(float complex));
+  float complex* H_herm =  malloc(N*N*sizeof(float complex));
+
+  H_hermH_plus_sigma2I(N, N, H, sigma2, H_hermH_sigmaI);
+
+#ifdef DEBUG_PREPROC
+  int i =0;
+  for(i=0;i<N*N;i++)
+    printf(" H_hermH_sigmaI[%d] = (%f + i%f)\n", i , creal(H_hermH_sigmaI[i]), cimag(H_hermH_sigmaI[i]));
+#endif
+
+  conjugate_transpose (N, H, H_herm); //equals H_herm
+
+#ifdef DEBUG_PREPROC
+  for(i=0;i<N*N;i++)
+    printf(" H_herm[%d] = (%f + i%f)\n", i , creal(H_herm[i]), cimag(H_herm[i]));
+#endif
+
+  lin_eq_solver(N, H_hermH_sigmaI, H_herm, W_MMSE);
+
+#ifdef DEBUG_PREPROC
+  for(i=0;i<N*N;i++)
+    printf(" W_MMSE[%d] = (%f + i%f)\n", i , creal(W_MMSE[i]), cimag(W_MMSE[i]));
+#endif
+
+  free(H_hermH_sigmaI);
+  free(H_herm);
+}
+
+#if 0
+void compute_white_filter(float complex* H_re,
+                          int order_H,
+                          float sigma2,
+                          float complex* W_Wh_0_re,
+                          float complex* W_Wh_1_re){
+
+  int aatx, aarx, re;
+  int i,j;
+  int M =n_rx;
+  int N = n_tx;
+  int sigma2=noise_power;
+
+  float complex *H0_re = malloc(n_rx*(n_tx>>2)*sizeof(float complex));
+  float complex *H1_re = malloc(n_rx*(n_tx>>2)*sizeof(float complex));
+  float complex *R_corr_col_n_0_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *R_corr_col_n_1_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *U_0_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *U_1_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *U_0_herm_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *U_1_herm_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *D_0_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *D_1_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *W_Wh_0_re = malloc(n_rx*n_tx*sizeof(float complex));
+  float complex *W_Wh_1_re = malloc(n_rx*n_tx*sizeof(float complex));
+
+  for (aatx=0; aatx<n_tx/2; aatx++){
+    for (aarx=0; aarx<n_rx; aarx++) {
+      H0_re[aatx*n_rx + aarx] = H_re[aatx*n_rx + aarx][re]; // H0 gets [0 1 2 3; 4,5,6,7].' coefficients of H
+      H1_re[aatx*n_rx + aarx] = H_re[aatx*n_rx + aarx + 8][re]; // H1 gets [8 9 10 11; 12, 13, 14, 15].' coefficients of H
+        if (re == 0)
+        printf("ant %d, H_re = (%f + i%f) \n", aatx*n_rx + aarx, creal(H[aatx*n_rx + aarx][re]), cimag(H[aatx*n_rx + aarx][re]));
+      }
+   }
+
+
+    //HH_herm_plus_sigma2I(n_rx, (n_tx>>2), H1_re, sigma2, R_corr_col_n_0_re);
+    HH_herm_plus_sigma2I(n_rx, (n_tx>>2), H0_re, sigma2, R_corr_col_n_1_re);
+
+    eigen_vectors_values(n_rx, R_corr_col_n_0_re, U_0_re, D_0_re);
+    eigen_vectors_values(n_rx, R_corr_col_n_1_re, U_1_re, D_1_re);
+
+    transpose (n_rx, U_0_re, U_0_herm_re);
+    transpose (n_rx, U_1_re, U_1_herm_re);
+
+    sigma = (float)(sqrt((double)(sigma2)));
+
+    /*The inverse of a diagonal matrix is obtained by replacing each element in the diagonal with its reciprocal.
+    A square root of a diagonal matrix is given by the diagonal matrix, whose diagonal entries are just the square
+     roots of the original matrix.*/
+
+
+    D_0_re_inv_sqrt[0] = sqrt_float(1/D_0_re_inv[0]);
+    D_0_re_inv_sqrt[5] = sqrt_float(1/D_0_re_inv[5]);
+    D_0_re_inv_sqrt[10] = sqrt_float(1/D_0_re_inv[10]);
+    D_0_re_inv_sqrt[15] = sqrt_float(1/D_0_re_inv[15]);
+
+    D_1_re_inv[0] = sqrt_float(1/D_1_re_inv[0]);
+    D_1_re_inv[5] = sqrt_float(1/D_1_re_inv[5]);
+    D_1_re_inv[10] = sqrt_float(1/D_1_re_inv[10]);
+    D_1_re_inv[15] = sqrt_float(1/D_1_re_inv[15]);
+
+    now only to multiply
+
+  free(H0);
+  free(H1);
+  free(R_corr_col_n_0);
+  free(R_corr_col_n_1);
+}
+#endif
+
+float sqrt_float(float x, float sqrt_x)
+{
+  sqrt_x = (float)(sqrt((double)(x)));
+  return sqrt_x;
+}
\ No newline at end of file

openair1/PHY/LTE_UE_TRANSPORT/linear_preprocessing_rec.h

+
+#include<stdio.h>
+#include<math.h>
+#include<complex.h>
+#include <stdlib.h>
+#include "PHY/defs_UE.h"
+
+
+/* FUNCTIONS FOR LINEAR PREPROCESSING: MMSE, WHITENNING, etc*/
+void transpose(int N, float complex *A, float complex *Result);
+
+void conjugate_transpose(int N, float complex *A, float complex *Result);
+
+void H_hermH_plus_sigma2I(int N, int M, float complex *A, float sigma2, float complex *Result);
+
+void HH_herm_plus_sigma2I(int M, int N, float complex *A, float sigma2, float complex *Result);
+
+void eigen_vectors_values(int N, float complex *A, float complex *Vectors, float *Values_Matrix);
+
+void lin_eq_solver(int N, float complex *A, float complex* B);
+//float complex* lin_eq_solver (int N, float complex* A, float complex* B);
+
+/* mutl_matrix_matrix_row_based performs multiplications when matrix is row-oriented H[0], H[1]; H[2], H[3]*/
+void mutl_matrix_matrix_row_based(float complex* M0, float complex* M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex* Result );
+
+/* mutl_matrix_matrix_col_based performs multiplications matrix is column-oriented H[0], H[2]; H[1], H[3]*/
+void mutl_matrix_matrix_col_based(float complex* M0, float complex* M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex* Result );
+
+void compute_MMSE(float complex* H, int order_H, float sigma2, float complex* W_MMSE);
+
+void compute_white_filter(float complex* H, int order_H, float sigma2, float complex* U_1, float complex* D_1);
+
+void mmse_processing_oai(LTE_UE_PDSCH *pdsch_vars,
+                         LTE_DL_FRAME_PARMS *frame_parms,
+                         PHY_MEASUREMENTS *measurements,
+                         unsigned char first_symbol_flag,
+                         MIMO_mode_t mimo_mode,
+                         unsigned short mmse_flag,
+                         int noise_power,
+                         unsigned char symbol,
+                         unsigned short nb_rb);
+
+
+void precode_channel_est(int32_t **dl_ch_estimates_ext,
+                        LTE_DL_FRAME_PARMS *frame_parms,
+                        LTE_UE_PDSCH *pdsch_vars,
+                        unsigned char symbol,
+                        unsigned short nb_rb,
+                        MIMO_mode_t mimo_mode);
+
+
+void rxdataF_to_float(int32_t **rxdataF_ext,
+                      float complex **rxdataF_f,
+                      int n_rx,
+                      int length,
+                      int start_point);
+
+void chan_est_to_float(int32_t **dl_ch_estimates_ext,
+                       float complex **dl_ch_estimates_ext_f,
+                       int n_tx,
+                       int n_rx,
+                       int length,
+                       int start_point);
+
+void float_to_chan_est(int32_t **dl_ch_estimates_ext,
+                      float complex **dl_ch_estimates_ext_f,
+                      int n_tx,
+                      int n_rx,
+                      int length,
+                      int start_point);
+
+void float_to_rxdataF(int32_t **rxdataF_ext,
+                      float complex **rxdataF_f,
+                      int n_tx,
+                      int n_rx,
+                      int length,
+                      int start_point);
+
+void mult_mmse_rxdataF(float complex** Wmmse,
+                       float complex** rxdataF_ext_f,
+                       int n_tx,
+                       int n_rx,
+                       int length,
+                       int start_point);
+
+void mult_mmse_chan_est(float complex** Wmmse,
+                        float complex** dl_ch_estimates_ext_f,
+                        int n_tx,
+                        int n_rx,
+                        int length,
+                        int start_point);
+
+void  mmse_processing_core(int32_t **rxdataF_ext,
+                           int32_t **dl_ch_estimates_ext,
+                           int sigma2,
+                           int n_tx,
+                           int n_rx,
+                           int length,
+                           int start_point);
+
+void mmse_processing_core_flp(float complex** rxdataF_ext_flcpx,
+                              float complex **H,
+                              int32_t **rxdataF_ext,
+                              int32_t **dl_ch_estimates_ext,
+                              float sigma2,
+                              int n_tx,
+                              int n_rx,
+                              int length,
+                              int start_point);
+
+void whitening_processing_core_flp(float complex** rxdataF_ext_flcpx,
+                                   float complex **H,
+                                   int32_t **rxdataF_ext,
+                                   int32_t **dl_ch_estimates_ext,
+                                   float sigma2,
+                                   int n_tx,
+                                   int n_rx,
+                                   int length,
+                                   int start_point);
+
+float sqrt_float(float x, float sqrt_x);

openair1/PHY/LTE_UE_TRANSPORT/transport_proto_ue.h

@@ -925,6 +925,7 @@ void dlsch_channel_compensation_TM34(LTE_DL_FRAME_PARMS *frame_parms,
                                     int round,
                                     MIMO_mode_t mimo_mode,
                                     unsigned short nb_rb,
+                                    unsigned short mmse_flag,
                                     unsigned char output_shift0,
                                     unsigned char output_shift1);
 
@@ -957,6 +958,7 @@ void dlsch_channel_level_TM34(int **dl_ch_estimates_ext,
                               int *avg_1,
                               uint8_t symbol,
                               unsigned short nb_rb,
+                              unsigned int mmse_flag,
                               MIMO_mode_t mimo_mode);