do not print CPU value.

dlsch_compensation_core.
Minor changes in openair1/PHY/LTE_UE_TRANSPORT/linear_preprocessing_rec.c

cmake_targets/CMakeLists.txt

@@ -1215,6 +1215,12 @@ set(PHY_MEX_UE
   ${OPENAIR1_DIR}/PHY/TOOLS/signal_energy.c
   ${OPENAIR1_DIR}/PHY/LTE_ESTIMATION/lte_ue_measurements.c
   ${OPENAIR2_DIR}/UTIL/LOG/log.c
+  ${OPENAIR_DIR}/common/utils/T/T.c
+  ${OPENAIR_DIR}/common/utils/T/local_tracer.c
+  ${OPENAIR_DIR}/common/config/config_cmdline.c
+  ${OPENAIR_DIR}/common/config/config_userapi.c
+  ${OPENAIR_DIR}/common/config/config_load_configmodule.c
+  ${OPENAIR1_DIR}/PHY/TOOLS/time_meas.c
   )
 add_library(PHY_MEX ${PHY_MEX_UE})
 
@@ -2182,7 +2188,7 @@ if (${T_TRACER})
         oai_eth_transpro
         FLPT_MSG ASYNC_IF FLEXRAN_AGENT HASHTABLE MSC UTIL OMG_SUMO SECU_OSA
         SECU_CN SCHED_LIB PHY L2 default_sched remote_sched RAL CN_UTILS
-        GTPV1U SCTP_CLIENT UDP LIB_NAS_UE LFDS LFDS7 SIMU OPENAIR0_LIB)
+        GTPV1U SCTP_CLIENT UDP LIB_NAS_UE LFDS LFDS7 SIMU OPENAIR0_LIB PHY_MEX)
     if (TARGET ${i})
       add_dependencies(${i} generate_T)
     endif()

openair1/PHY/LTE_UE_TRANSPORT/linear_preprocessing_rec.c

@@ -59,42 +59,26 @@ void transpose (int N, float complex *A, float complex *Result)
  }
 
 
-void conjugate_transpose (int N, float complex *A, float complex *Result)
+
+void conjugate_transpose (int rows_A, int col_A, 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;
+  float complex alpha = 1.0;
+  float complex beta = 0.0;
   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));
+  B = (float complex*)calloc(rows_A*rows_A,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);
+  for (i=0; i<rows_A*rows_A; i+=rows_A+1)
+    B[i]=1.0;
 
-  memcpy(Result, C, N*N*sizeof(float complex));
+  cblas_cgemm(CblasColMajor, transa, transb, col_A, rows_A, rows_A, &alpha, A, rows_A, B, rows_A, &beta, Result, col_A);
 
   free(B);
-  free(C);
  }
 
 void H_hermH_plus_sigma2I (int N, int M, float complex *A, float sigma2, float complex *Result)
@@ -121,30 +105,23 @@ void H_hermH_plus_sigma2I (int N, int M, float complex *A, float sigma2, float c
 
   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)
+
+ void HH_herm_plus_sigma2I (int rows_A, int col_A, 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 < rows_A*rows_A; i += rows_A+1)
+    Result[i]=1.0+I*0;
 
-  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);
+  cblas_cgemm(CblasColMajor, transa, transb, rows_A, rows_A, col_A, &alpha, A, rows_A, A, rows_A, &sigma2, Result, rows_A);
 
 }
 
@@ -153,14 +130,14 @@ void eigen_vectors_values (int N, float complex *A, float complex *Vectors, floa
   // 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 jobz = 'V'; // compute both eigenvectors and eigenvalues
   char uplo = 'U';
   int order_A = N;
   int lda = N;
   int i;
-  float* Values = (float*)malloc(sizeof(float)*1*N);
+  float* Values = (float*)calloc(1*N, sizeof(float));
 
-  LAPACKE_cheev(LAPACK_ROW_MAJOR, jobz, uplo, order_A, A, lda, Values);
+  LAPACKE_cheev(LAPACK_COL_MAJOR, jobz, uplo, order_A, A, lda, Values);
 
   memcpy(Vectors, A, N*N*sizeof(float complex));
 
@@ -178,7 +155,7 @@ void eigen_vectors_values (int N, float complex *A, float complex *Vectors, floa
   int nrhs = N;
 
   char transa = 'N';
-  int* IPIV = malloc(N*N*sizeof(int));
+  int* IPIV = calloc(N*N, sizeof(int));
 
   // Compute LU-factorization
   LAPACKE_cgetrf(LAPACK_ROW_MAJOR, n, nrhs, A, lda, IPIV);
@@ -208,23 +185,24 @@ void mutl_matrix_matrix_row_based(float complex* M0, float complex* M1, int rows
 
 #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);
+  printf("mutl_matrix_matrix_row_based: 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]));
+    printf("mutl_matrix_matrix_row_based: 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]));
+    printf("mutl_matrix_matrix_row_based: 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]));
+    printf("mutl_matrix_matrix_row_based: 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;
@@ -239,131 +217,221 @@ void mutl_matrix_matrix_col_based(float complex* M0, float complex* M1, int rows
 
 #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);
+  printf("mutl_matrix_matrix_col_based: 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_M0*col_M0; ++i)
+    printf("mutl_matrix_matrix_col_based: rows_opA = %d, col_opB = %d, filter[%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]));
+  for(i = 0; i < rows_M1*col_M1; ++i)
+    printf("mutl_matrix_matrix_col_based: 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]));
+  for(i = 0; i < rows_opA*col_opB; ++i)
+    printf("mutl_matrix_matrix_col_based: result[%d] = (%f + i%f)\n", i , creal(Result[i]), cimag(Result[i]));
 #endif
 }
 
+void mutl_scal_matrix_matrix_col_based(float *M0, float complex *M1, float alpha, 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;
+  float complex beta = 0.0;
+  int i;
+
+  // Convert float M0 into complex float D_0_complex required by cblas_cgemm
+  float complex *D_0_complex = calloc(rows_M0*col_M0, sizeof(float complex));
+  for(i = 0; i < rows_M0*col_M0; ++i)
+  {
+    D_0_complex[i] = M0[i] + I*0.00001;
+#ifdef DEBUG_PREPROC
+    printf("mutl_scal_matrix_matrix_col_based: D_0_complex[%d] = (%f, %f)\n", i , creal(D_0_complex[i]), cimag(D_0_complex[i]));
+#endif
+  }
+
+#ifdef DEBUG_PREPROC
+  printf("mutl_scal_matrix_matrix_col_based: alpha = %f\n", alpha);
+
+  for(i = 0; i < rows_M0*col_M0; ++i){
+    printf("mutl_scal_matrix_matrix_col_based M0[%d] = %f\n", i , M0[i]);
+  }
+
+  for(i = 0; i < rows_M1*col_M1; ++i)
+    printf("mutl_scal_matrix_matrix_col_based: M1[%d] = (%f + i%f)\n", i , creal(M1[i]), cimag(M1[i]));
+#endif
+
+  cblas_cgemm(CblasColMajor, transa, transb, rows_M0, col_M1, col_M0, &alpha, D_0_complex, rows_M0, M1, col_M0, &beta, Result, rows_M0);
+
+#ifdef DEBUG_PREPROC
+  for(i = 0; i < rows_M0*col_M1; ++i)
+    printf("mutl_scal_matrix_matrix_col_based: result[%d] = (%f + i%f)\n", i , creal(Result[i]), cimag(Result[i]));
+#endif
+
+  free(D_0_complex);
+}
+
 
 /*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));
+  float complex* H_hermH_sigmaI = calloc(N*N, sizeof(float complex));
+  float complex* H_herm =  calloc(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]));
+  int i = 0;
+  for(i = 0;i < N*N; ++i)
+    printf("compute_MMSE: 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
+  conjugate_transpose (N, 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]));
+  for(i = 0;i < N*N;i++)
+    printf("compute_MMSE: 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]));
+  for(i = 0;i < N*N; ++i)
+    printf("compute_MMSE: 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 sqrt_float(float x)
+{
+  float sqrt_x = 0.0;
+  sqrt_x = (float)(sqrt((double)(x)));
+  return sqrt_x;
+}
+
+void compute_white_filter(float complex* H0_re,
+                          float complex* H1_re,
                           float sigma2,
+                          int n_rx,
+                          int n_tx,
                           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 sigma = 0.0;
+  int i;
 
-  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));
+  float complex *R_corr_col_n_0_re = calloc(n_rx*n_tx, sizeof(float complex));
+  float complex *R_corr_col_n_1_re = calloc(n_rx*n_tx, sizeof(float complex));
+  float complex *U_0_re = calloc(n_rx*n_tx, sizeof(float complex));
+  float complex *U_1_re = calloc(n_rx*n_tx, sizeof(float complex));
+  float complex *U_0_herm_re = calloc(n_rx*n_tx, sizeof(float complex));
+  float complex *U_1_herm_re = calloc(n_rx*n_tx, sizeof(float complex));
+  float *D_0_re = calloc(n_rx*n_tx, sizeof(float));
+  float *D_1_re = calloc(n_rx*n_tx, sizeof(float));
+  float *D_0_re_inv_sqrt = calloc(n_rx*n_tx, sizeof(float));
+  float *D_1_re_inv_sqrt = calloc(n_rx*n_tx, sizeof(float));
 
-  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]));
-      }
-   }
+  // Whitening filter can be computed using the following algorithm:
+  // 1. Compute covariance of the colored noise: R = HH' + sigma2I.
+  // 2. Compute eigen value decomposition of R = UDU'.
+  // 3. W_wh = sigma sqrt(inv(D))U'.
+  // 4. This function computes W_wh for both branches.
 
-
-    //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);
+#ifdef DEBUG_PREPROC
+  printf("compute_white_filter: sigma2 = %f\n", sigma2);
+  for(i=0; i<n_rx*n_tx/2; i++){
+    printf("compute_white_filter: H1_re[%d] = (%f + i%f)\n", i , creal(H1_re[i]), cimag(H1_re[i]));
+    printf("compute_white_filter: H0_re[%d] = (%f + i%f)\n", i , creal(H0_re[i]), cimag(H0_re[i]));
 }
 #endif
 
-float sqrt_float(float x, float sqrt_x)
-{
-  sqrt_x = (float)(sqrt((double)(x)));
-  return sqrt_x;
-}
+
+  // 1. Compute covariance of the colored noise: R = HH' + sigma2I.
+  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);
+
+#ifdef DEBUG_PREPROC
+  for(i=0;i<n_rx*n_tx;i++){
+    printf("compute_white_filter: R_corr_col_n_0_re[%d] = (%f + i%f)\n", i , creal(R_corr_col_n_0_re[i]), cimag(R_corr_col_n_0_re[i]));
+    printf("compute_white_filter: R_corr_col_n_1_re[%d] = (%f + i%f)\n", i , creal(R_corr_col_n_1_re[i]), cimag(R_corr_col_n_1_re[i]));
+  }
+
+#endif
+  // 2. Compute eigen value decomposition of R = UDU'.
+  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);
+
+#ifdef DEBUG_PREPROC
+  for(i=0;i<n_rx*n_tx;i++){
+    printf("compute_white_filter: U_0_re[%d] = (%f + i%f)\n", i , creal(U_0_re[i]), cimag(U_0_re[i]));
+    printf("compute_white_filter: D_0_re[%d] = (%f + i%f)\n", i , creal(D_0_re[i]), cimag(D_0_re[i]));
+}
+#endif
+
+  // 3. Compute eigen value decomposition of R = UDU'.
+
+  conjugate_transpose(n_rx, n_tx, U_0_re, U_0_herm_re);
+  conjugate_transpose(n_rx, n_tx, U_1_re, U_1_herm_re);
+
+
+#ifdef DEBUG_PREPROC
+  for(i = 0;i < n_rx*n_tx; i++){
+    printf("compute_white_filter: U_0_herm_re[%d] = (%f + i%f)\n", i , creal(U_0_herm_re[i]), cimag(U_0_herm_re[i]));
+  }
+#endif
+
+
+  sigma = (float)(sqrt((double)(sigma2)));
+  if (sigma <= 0.0001){
+      sigma = 0.0001;
+    }
+
+  //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. However, if SNR is high,
+  //the diagonal elements of D are very small, and inverse is not always possible. We thus appy a threshold to avoid too low values.
+
+  for (i = 0; i < n_rx*n_tx; i += (n_rx + 1)){
+
+    if (D_0_re[i] <= 0.0001){
+      D_0_re[i] = 0.0001;
+    }
+    if (D_1_re[i] <= 0.0001){
+      D_1_re[i] = 0.0001;
+    }
+
+    D_0_re_inv_sqrt[i] = sqrt_float(1/D_0_re[i]);
+    D_1_re_inv_sqrt[i] = sqrt_float(1/D_1_re[i]);
+  }
+
+
+#ifdef DEBUG_PREPROC
+  for(i = 0;i <n_rx*n_tx; i++){
+    printf("compute_white_filter: D_0_re_inv_sqrt[%d] = %f\n", i , D_0_re_inv_sqrt[i]);
+  }
+#endif
+
+
+  mutl_scal_matrix_matrix_col_based(D_0_re_inv_sqrt, U_0_herm_re, sigma, n_rx, n_tx, n_rx, n_tx, W_Wh_0_re);
+
+#ifdef DEBUG_PREPROC
+  for(i = 0;i < n_rx*n_tx; i++){
+    printf("compute_white_filter: W_Wh_0_re[%d] = (%f + i%f)\n", i , creal(W_Wh_0_re[i]), cimag(W_Wh_0_re[i]));
+  }
+#endif
+
+  mutl_scal_matrix_matrix_col_based(D_1_re_inv_sqrt, U_1_herm_re, sigma, n_rx, n_tx, n_rx, n_tx, W_Wh_1_re);
+
+  free(R_corr_col_n_0_re);
+  free(R_corr_col_n_1_re);
+  free(U_0_herm_re);
+  free(U_1_herm_re);
+  free(D_0_re_inv_sqrt);
+  free(D_1_re_inv_sqrt);
+  free(U_0_re);
+  free(U_1_re);
+  free(D_0_re);
+  free(D_1_re);
+}

openair1/PHY/LTE_UE_TRANSPORT/linear_preprocessing_rec.h

@@ -9,26 +9,35 @@
 /* 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 conjugate_transpose (int rows_A, int col_A, float complex *A, float complex *Result);
 
-void H_hermH_plus_sigma2I(int N, int M, float complex *A, float sigma2, float complex *Result);
+void H_hermH_plus_sigma2I (int row_A, int col_A, 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 HH_herm_plus_sigma2I (int rows_A, int col_A, 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);
+void 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 );
+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 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 mutl_scal_matrix_matrix_col_based(float complex *M0, float complex *M1, float complex alpha, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex *Result);
 
-void compute_white_filter(float complex* H, int order_H, float sigma2, float complex* U_1, float complex* D_1);
+void compute_MMSE(float complex *H, int order_H, float sigma2, float complex *W_MMSE);
+
+float sqrt_float(float x);
+
+void compute_white_filter(float complex *H0_re,
+                          float complex *H1_re,
+                          float sigma2,
+                          int n_rx,
+                          int n_tx,
+                          float complex *W_Wh_0_re,
+                          float complex *W_Wh_1_re);
 
 void mmse_processing_oai(LTE_UE_PDSCH *pdsch_vars,
                          LTE_DL_FRAME_PARMS *frame_parms,
@@ -57,65 +66,71 @@ void rxdataF_to_float(int32_t **rxdataF_ext,
 
 void chan_est_to_float(int32_t **dl_ch_estimates_ext,
                        float complex **dl_ch_estimates_ext_f,
+                       uint8_t n_tx,
+                       uint8_t n_rx,
+                       int32_t length,
+                       int32_t start_point);
+
+void float_to_chan_est(float complex **chan_est_flp,
+                       int32_t **result,
                        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(float complex **rxdataF_flp,
+                      int32_t **result,
+                      uint8_t n_tx,
+                      uint8_t n_rx,
+                      int32_t length,
+                      int32_t 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_filter_chan_est(float complex **W,
+                          float complex **chan_est_flp,
+                          float complex **result,
+                          uint8_t n_tx,
+                          uint8_t n_rx,
+                          int32_t n_col_chan_est_flp,
+                          int32_t length,
+                          int32_t 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_filter_rxdataF(float complex **W,
+                         float complex **rxdataF_flp,
+                         float complex **result,
+                         uint8_t n_tx,
+                         uint8_t n_rx,
+                         int32_t length,
+                         int32_t 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(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_flp,
+                              float complex **chan_est_flp,
+                              int32_t **rxdataF_filt_fp,
+                              int32_t **chan_est_eff_fp,
+                              float noise_power,
+                              uint8_t n_tx,
+                              uint8_t n_rx,
+                              int32_t length,
+                              int32_t 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,
+void whitening_processing_core_flp(float complex **rxdataF_flp,
+                                   float complex **chan_est_flp_0,
+                                   float complex **chan_est_flp_1,
+                                   int32_t **rxdataF_filt_fp_0,
+                                   int32_t **rxdataF_filt_fp_1,
+                                   int32_t **chan_est_eff_fp_0,
+                                   int32_t **chan_est_eff_fp_1,
                                    float sigma2,
-                                   int n_tx,
-                                   int n_rx,
-                                   int length,
-                                   int start_point);
+                                   uint8_t n_tx,
+                                   uint8_t n_rx,
+                                   int32_t length,
+                                   int32_t start_point);
+
 
-float sqrt_float(float x, float sqrt_x);

openair1/PHY/LTE_UE_TRANSPORT/transport_proto_ue.h

@@ -691,6 +691,19 @@ void dlsch_detection_mrc(LTE_DL_FRAME_PARMS *frame_parms,
                          uint16_t nb_rb,
                          uint8_t dual_stream_UE);
 
+void dlsch_detection_mrc_core(int **rxdataF_comp,
+                              int **rxdataF_comp_i,
+                              int **rho,
+                              int **rho_i,
+                              int **dl_ch_mag,
+                              int **dl_ch_magb,
+                              int **dl_ch_mag_i,
+                              int **dl_ch_magb_i,
+                              unsigned char n_tx,
+                              unsigned char n_rx,
+                              int length,
+                              int start_point);
+
 void dlsch_detection_mrc_TM34(LTE_DL_FRAME_PARMS *frame_parms,
                               LTE_UE_PDSCH *lte_ue_pdsch_vars,
                               int harq_pid,
@@ -853,6 +866,15 @@ void dlsch_dual_stream_correlation(LTE_DL_FRAME_PARMS *frame_parms,
                                    int **dl_ch_rho_ext,
                                    unsigned char output_shift);
 
+void dlsch_dual_stream_correlation_core(int **dl_ch_estimates_ext,
+                                        int **dl_ch_estimates_ext_i,
+                                        int **dl_ch_rho_ext,
+                                        unsigned char n_tx,
+                                        unsigned char n_rx,
+                                        unsigned char output_shift,
+                                        int length,
+                                        int start_point);
+
 void dlsch_dual_stream_correlationTM34(LTE_DL_FRAME_PARMS *frame_parms,
                                    unsigned char symbol,
                                    unsigned short nb_rb,

openair1/PHY/TOOLS/time_meas.c

@@ -37,8 +37,8 @@ double get_cpu_freq_GHz(void) {
   sleep(1);
   ts.diff = (rdtsc_oai()-ts.in);
   cpu_freq_GHz = (double)ts.diff/1000000000;
-  printf("CPU Freq is %f \n", cpu_freq_GHz);
-  return cpu_freq_GHz; 
+  //printf("CPU Freq is %f \n", cpu_freq_GHz);
+  return cpu_freq_GHz;
 }
 
 void print_meas_now(time_stats_t *ts, const char* name, FILE* file_name){
@@ -54,7 +54,7 @@ void print_meas_now(time_stats_t *ts, const char* name, FILE* file_name){
 
       //fprintf(file_name,"Name %25s: Processing %15.3f ms for SF %d, diff_now %15.3f \n", name,(ts->diff_now/(cpu_freq_GHz*1000000.0)),subframe,ts->diff_now);
       fprintf(file_name,"%15.3f us, diff_now %15.3f \n",(ts->diff_now/(cpu_freq_GHz*1000.0)),(double)ts->diff_now);
-      
+
     }
   }
 }