Added PSS/SSS mex. Added new SSS correlation algorithms. Use PSR in PSS

lte/examples/synch_file.c

@@ -221,7 +221,7 @@ int main(int argc, char **argv) {
 
       sss_idx = peak_pos[N_id_2]-2*(symbol_sz+CP(symbol_sz,CPNORM_LEN));
       if (sss_idx >= 0) {
-        sss_synch_m0m1(&sss[N_id_2], &input[sss_idx],
+        sss_synch_m0m1_diff(&sss[N_id_2], &input[sss_idx],
             &m0, &m0_value, &m1, &m1_value);
 
         cfo[frame_cnt] = pss_synch_cfo_compute(&pss[N_id_2], &input[peak_pos[N_id_2]-128]);

lte/phy/include/liblte/phy/common/phy_common.h

@@ -93,10 +93,11 @@ typedef enum {CPNORM, CPEXT} lte_cp_t;
 #define SLOT_LEN_RE(nof_prb, cp)        (nof_prb*RE_X_RB*CP_NSYMB(cp))
 #define SF_LEN_RE(nof_prb, cp)          (2*SLOT_LEN_RE(nof_prb, cp))
 
-#define SLOT_IDX_CPNORM(idx, symbol_sz) (idx==0?(CP(symbol_sz, CPNORM_0_LEN)):(CP(symbol_sz, CPNORM_0_LEN)+idx*(symbol_sz+CP(symbol_sz, CPNORM_LEN))))
+#define SLOT_IDX_CPNORM(symbol_idx, symbol_sz) (symbol_idx==0?0:(symbol_sz + CP(symbol_sz, CPNORM_0_LEN) + \
+                                                (symbol_idx-1)*(symbol_sz+CP(symbol_sz, CPNORM_LEN))))
 #define SLOT_IDX_CPEXT(idx, symbol_sz) (idx*(symbol_sz+CP(symbol_sz, CPEXT_LEN)))
 
-#define SAMPLE_IDX(nof_prb, symbol_idx, sample_idx) ((symbol_idx)*(nof_prb)*(RE_X_RB) + sample_idx)
+#define RE_IDX(nof_prb, symbol_idx, sample_idx) ((symbol_idx)*(nof_prb)*(RE_X_RB) + sample_idx)
 
 #define RS_VSHIFT(cell_id) (cell_id%6)
 

lte/phy/include/liblte/phy/sync/pss.h

@@ -46,6 +46,13 @@ typedef _Complex float cf_t; /* this is only a shortcut */
 #define PSS_RE      6*12
 
 
+/* PSS processing options */
+
+//#define PSS_ACCUMULATE_ABS   // If enabled, accumulates the correlation absolute value on consecutive calls to pss_synch_find_pss
+
+#define PSS_ABS_SQUARE   // If enabled, compute abs square, otherwise computes absolute value only 
+
+#define PSS_RETURN_PSR  // If enabled returns peak to side-lobe ratio, otherwise returns absolute peak value
 
 /**
  * The pss_synch_t object provides functions for fast computation of the crosscorrelation
@@ -62,6 +69,8 @@ typedef _Complex float cf_t; /* this is only a shortcut */
 /* Low-level API */
 typedef struct LIBLTE_API {
   
+  dft_plan_t dftp_input; 
+  
 #ifdef CONVOLUTION_FFT
   conv_fft_cc_t conv_fft;
 #endif
@@ -70,10 +79,14 @@ typedef struct LIBLTE_API {
   uint32_t N_id_2;
   uint32_t fft_size;
 
+  cf_t pss_signal_time[3][PSS_LEN];
   cf_t *pss_signal_freq[3]; // One sequence for each N_id_2
   cf_t *tmp_input;
   cf_t *conv_output;
-    
+#ifdef PSS_ACCUMULATE_ABS
+  float *conv_output_abs;
+#endif
+  float *conv_output_avg;
 }pss_synch_t;
 
 typedef enum { PSS_TX, PSS_RX } pss_direction_t;
@@ -88,6 +101,8 @@ LIBLTE_API int pss_synch_init(pss_synch_t *q,
 
 LIBLTE_API void pss_synch_free(pss_synch_t *q);
 
+LIBLTE_API void pss_synch_reset(pss_synch_t *q); 
+
 LIBLTE_API int pss_generate(cf_t *signal, 
                             uint32_t N_id_2);
 
@@ -103,6 +118,10 @@ LIBLTE_API int pss_synch_find_pss(pss_synch_t *q,
                                   cf_t *input, 
                                   float *corr_peak_value);
 
+LIBLTE_API int pss_synch_chest(pss_synch_t *q, 
+                               cf_t *input, 
+                               cf_t ce[PSS_LEN]); 
+
 LIBLTE_API float pss_synch_cfo_compute(pss_synch_t* q, 
                                        cf_t *pss_recv);
 

lte/phy/include/liblte/phy/sync/sss.h

@@ -42,8 +42,6 @@ typedef _Complex float cf_t; /* this is only a shortcut */
 #define N_SSS      31
 #define SSS_LEN    2*N_SSS
 
-#define SSS_MAX_FFT_LEN 2048
-
 struct sss_tables{
   int z1[N_SSS][N_SSS];
   int c[2][N_SSS];
@@ -54,9 +52,10 @@ struct sss_tables{
  * Should use vec_malloc() to make it platform agnostic.
  */
 struct fc_tables{
-  cf_t z1[N_SSS+1][N_SSS+1];
+  float z1[N_SSS][N_SSS];
-  cf_t c[2][N_SSS+1];
+  float c[2][N_SSS];
-  cf_t s[N_SSS+1][N_SSS+1];
+  float s[N_SSS][N_SSS];
+  float sd[N_SSS][N_SSS-1];
 };
 
 
@@ -75,6 +74,9 @@ typedef struct LIBLTE_API {
   uint32_t N_id_1_table[30][30];
   struct fc_tables fc_tables[3]; // one for each N_id_2
 
+  float corr_output_m0[N_SSS];
+  float corr_output_m1[N_SSS];
+
 }sss_synch_t;
 
 
@@ -99,13 +101,31 @@ LIBLTE_API void sss_put_slot(float *sss,
 LIBLTE_API int sss_synch_set_N_id_2(sss_synch_t *q, 
                                     uint32_t N_id_2);
 
-LIBLTE_API int sss_synch_m0m1(sss_synch_t *q, 
+LIBLTE_API int sss_synch_m0m1_partial(sss_synch_t *q, 
+                                      cf_t *input, 
+                                      uint32_t M, 
+                                      cf_t ce[2*N_SSS],
+                                      uint32_t *m0, 
+                                      float *m0_value, 
+                                      uint32_t *m1, 
+                                      float *m1_value);
+
+LIBLTE_API int sss_synch_m0m1_diff_coh(sss_synch_t *q, 
                                        cf_t *input, 
+                                       cf_t ce[2*N_SSS],
                                        uint32_t *m0, 
                                        float *m0_value, 
                                        uint32_t *m1, 
                                        float *m1_value);
 
+LIBLTE_API int sss_synch_m0m1_diff(sss_synch_t *q, 
+                                   cf_t *input, 
+                                   uint32_t *m0, 
+                                   float *m0_value, 
+                                   uint32_t *m1, 
+                                   float *m1_value);
+
+
 LIBLTE_API uint32_t sss_synch_subframe(uint32_t m0, 
                                        uint32_t m1);
 

lte/phy/lib/ch_estimation/src/chest_dl.c

@@ -254,7 +254,7 @@ static void average_pilots(chest_dl_t *q, uint32_t port_id)
   }
 }
 
-#define cesymb(i) ce[SAMPLE_IDX(q->cell.nof_prb,i,0)]
+#define cesymb(i) ce[RE_IDX(q->cell.nof_prb,i,0)]
 
 static void interpolate_pilots(chest_dl_t *q, cf_t *ce, uint32_t port_id) 
 {

lte/phy/lib/ch_estimation/src/refsignal_dl.c

@@ -209,7 +209,7 @@ int refsignal_cs_put_sf(lte_cell_t cell, uint32_t port_id, cf_t *pilots, cf_t *s
       /* Compute offset frequency index */
       fidx = ((refsignal_cs_v(port_id, l) + (cell.id % 6)) % 6); 
       for (i = 0; i < 2*cell.nof_prb; i++) {
-        sf_symbols[SAMPLE_IDX(cell.nof_prb, nsymbol, fidx)] = pilots[REFSIGNAL_PILOT_IDX(i,l,cell)];
+        sf_symbols[RE_IDX(cell.nof_prb, nsymbol, fidx)] = pilots[REFSIGNAL_PILOT_IDX(i,l,cell)];
         fidx += RE_X_RB/2;       // 1 reference every 6 RE        
       }    
     }
@@ -238,7 +238,7 @@ int refsignal_cs_get_sf(lte_cell_t cell, uint32_t port_id, cf_t *sf_symbols, cf_
       /* Compute offset frequency index */
       fidx = ((refsignal_cs_v(port_id, l) + (cell.id % 6)) % 6); 
       for (i = 0; i < 2*cell.nof_prb; i++) {
-        pilots[REFSIGNAL_PILOT_IDX(i,l,cell)] = sf_symbols[SAMPLE_IDX(cell.nof_prb, nsymbol, fidx)];
+        pilots[REFSIGNAL_PILOT_IDX(i,l,cell)] = sf_symbols[RE_IDX(cell.nof_prb, nsymbol, fidx)];
         fidx += RE_X_RB/2;       // 2 references per PRB
       }
     }

lte/phy/lib/sync/src/find_sss.c

@@ -27,30 +27,82 @@
 
 
 #include <string.h>
+#include <complex.h>
+#include <math.h>
 
 #include "liblte/phy/utils/vector.h"
 #include "liblte/phy/sync/sss.h"
 
-void corr_all_zs(cf_t *z, cf_t s[N_SSS+1][N_SSS+1], cf_t *output) {
+#define MAX_M 3
+
+
+static void corr_all_zs(cf_t z[N_SSS], float s[N_SSS][N_SSS-1], float output[N_SSS]) {
   uint32_t m;
+  cf_t tmp[N_SSS]; 
+  
   for (m = 0; m < N_SSS; m++) {
-    output[m] = vec_dot_prod_ccc(z, s[m], N_SSS - 1);
+    tmp[m] = vec_dot_prod_cfc(z, s[m], N_SSS - 1);    
   }
+  vec_abs_square_cf(tmp, output, N_SSS);
 }
 
+static void corr_all_sz_partial(cf_t z[N_SSS], float s[N_SSS][N_SSS], uint32_t M, float output[N_SSS]) {
+  uint32_t Nm = N_SSS/M; 
+  cf_t tmp[N_SSS];
+  float tmp_abs[MAX_M-1][N_SSS];
+  int j, m; 
+  float *ptr; 
   
-/* Assumes input points to the beginning of the SSS symbol. The SSS symbol start is
+  for (j=0;j<M;j++) {
- * given by SSS_SYMBOL_ST() macro in sss.h.
+    for (m = 0; m < N_SSS; m++) {      
- * Estimates the m0 and m1 values and saves in m0_value and m1_value
+      tmp[m] = vec_dot_prod_cfc(&z[j*Nm], &s[m][j*Nm], Nm);        
- * the resulted correlation (higher is more likely)
+    }
- *
+    if (j == 0) {
+      ptr = output; 
+    } else {
+      ptr = tmp_abs[j-1];
+    }
+    vec_abs_square_cf(tmp, ptr, N_SSS);   
+  }
+  for (j=1;j<M;j++) {
+    vec_sum_fff(tmp_abs[j-1], output, output, N_SSS);
+  }    
+}
+
+static void extract_pair_sss(sss_synch_t *q, cf_t *input, cf_t *ce, cf_t y[2][N_SSS]) {
+  cf_t input_fft[SYMBOL_SZ_MAX];
+
+  
+  dft_run_c(&q->dftp_input, &input[CP_NORM(5, q->fft_size)], input_fft);
+  
+  if (ce) {
+    vec_prod_conj_ccc(&input_fft[q->fft_size/2-N_SSS], ce, &input_fft[q->fft_size/2-N_SSS], 2*N_SSS);
+  }
+  
+  for (int i = 0; i < N_SSS; i++) {
+    y[0][i] = input_fft[q->fft_size/2-N_SSS + 2 * i];
+    y[1][i] = input_fft[q->fft_size/2-N_SSS + 2 * i + 1];
+  }
+
+  vec_prod_cfc(y[0], q->fc_tables[q->N_id_2].c[0], y[0], N_SSS);
+  vec_prod_cfc(y[1], q->fc_tables[q->N_id_2].c[1], y[1], N_SSS);
+
+}    
+
+int sss_synch_m0m1_diff(sss_synch_t *q, cf_t *input, uint32_t *m0, float *m0_value,
+    uint32_t *m1, float *m1_value) {
+  return sss_synch_m0m1_diff_coh(q, input, NULL, m0, m0_value, m1, m1_value);
+}
+
+/* Differential SSS estimation. 
+ * Returns m0 and m1 estimates 
  *
  * Source: "SSS Detection Method for Initial Cell Search in 3GPP LTE FDD/TDD Dual Mode Receiver"
  *       Jung-In Kim, Jung-Su Han, Hee-Jin Roh and Hyung-Jin Choi
 
  *
  */
-int sss_synch_m0m1(sss_synch_t *q, cf_t *input, uint32_t *m0, float *m0_value,
+int sss_synch_m0m1_diff_coh(sss_synch_t *q, cf_t *input, cf_t ce[2*N_SSS], uint32_t *m0, float *m0_value,
     uint32_t *m1, float *m1_value) 
 {
 
@@ -62,43 +114,63 @@ int sss_synch_m0m1(sss_synch_t *q, cf_t *input, uint32_t *m0, float *m0_value,
       m1                != NULL)
   {
     
-    /* Consider moving to sss_synch_t?? */
+    cf_t yprod[N_SSS];
-    cf_t zdelay[N_SSS+1],zconj[N_SSS+1],zprod[N_SSS+1];
+    cf_t y[2][N_SSS];
-    cf_t y[2][N_SSS+1], z[N_SSS+1], tmp[N_SSS+1];
-    float tmp_real[N_SSS+1];
-    cf_t input_fft[SSS_MAX_FFT_LEN];
-    uint32_t i;
 
-    dft_run_c(&q->dftp_input, input, input_fft);
+    extract_pair_sss(q, input, ce, y);
     
-    for (i = 0; i < N_SSS; i++) {
+    vec_prod_conj_ccc(&y[0][1], y[0], yprod, N_SSS - 1);    
-      y[0][i] = input_fft[q->fft_size/2-N_SSS + 2 * i];
+    corr_all_zs(yprod, q->fc_tables[q->N_id_2].sd, q->corr_output_m0);
-      y[1][i] = input_fft[q->fft_size/2-N_SSS + 2 * i + 1];
+    *m0 = vec_max_fi(q->corr_output_m0, N_SSS);
+    if (m0_value) {
+      *m0_value = q->corr_output_m0[*m0];
     }    
     
-    vec_prod_ccc(y[0], q->fc_tables[q->N_id_2].c[0], z, N_SSS);
+    vec_prod_cfc(y[1], q->fc_tables[q->N_id_2].z1[*m0], y[1], N_SSS);
-    memcpy(zdelay, &z[1], (N_SSS - 1) * sizeof(cf_t));
+    vec_prod_conj_ccc(&y[1][1], y[1], yprod, N_SSS - 1);
-    vec_conj_cc(z, zconj, N_SSS - 1);
+    corr_all_zs(yprod, q->fc_tables[q->N_id_2].sd, q->corr_output_m1);
-    vec_prod_ccc(zdelay, zconj, zprod, N_SSS - 1);
+    *m1 = vec_max_fi(q->corr_output_m1, N_SSS);
-
+    if (m1_value) {
-    corr_all_zs(zprod, q->fc_tables[q->N_id_2].s, tmp);
+      *m1_value = q->corr_output_m1[*m1];
-    vec_abs_cf(tmp, tmp_real, N_SSS);
+    }    
-    *m0 = vec_max_fi(tmp_real, N_SSS);
+    ret = LIBLTE_SUCCESS;
-    if (m0_value) {
+  } 
-      *m0_value = tmp_real[*m0];
+  return ret; 
 }
 
-    vec_prod_ccc(y[1], q->fc_tables[q->N_id_2].c[1], tmp, N_SSS);
+/* Partial correlation SSS estimation. 
-    vec_prod_ccc(tmp, q->fc_tables[q->N_id_2].z1[*m0], z, N_SSS);
+ * Returns m0 and m1 estimates 
-    memcpy(zdelay, &z[1], (N_SSS - 1) * sizeof(cf_t));
+ *
-    vec_conj_cc(z, zconj, N_SSS - 1);
+ * Source: "SSS Detection Method for Initial Cell Search in 3GPP LTE FDD/TDD Dual Mode Receiver"
-    vec_prod_ccc(zdelay, zconj, zprod, N_SSS - 1);
+ *       Jung-In Kim, Jung-Su Han, Hee-Jin Roh and Hyung-Jin Choi
+
+ */
+int sss_synch_m0m1_partial(sss_synch_t *q, cf_t *input, uint32_t M, cf_t ce[2*N_SSS], uint32_t *m0, float *m0_value,
+    uint32_t *m1, float *m1_value) 
+{
+
+  int ret = LIBLTE_ERROR_INVALID_INPUTS; 
+
+  if (q                 != NULL   &&
+      input             != NULL   &&
+      m0                != NULL   && 
+      m1                != NULL   && 
+      M                 <= MAX_M)
+  {
+    cf_t y[2][N_SSS];
+    
+    extract_pair_sss(q, input, ce, y);
     
-    corr_all_zs(zprod, q->fc_tables[q->N_id_2].s, tmp);
+    corr_all_sz_partial(y[0], q->fc_tables[q->N_id_2].s, M, q->corr_output_m0);
-    vec_abs_cf(tmp, tmp_real, N_SSS);
+    *m0 = vec_max_fi(q->corr_output_m0, N_SSS);
-    *m1 = vec_max_fi(tmp_real, N_SSS);
+    if (m0_value) {
+      *m0_value = q->corr_output_m0[*m0];
+    }
+    vec_prod_cfc(y[1], q->fc_tables[q->N_id_2].z1[*m0], y[1], N_SSS);
+    corr_all_sz_partial(y[1], q->fc_tables[q->N_id_2].s, M, q->corr_output_m1);
+    *m1 = vec_max_fi(q->corr_output_m1, N_SSS);
     if (m1_value) {
-      *m1_value = tmp_real[*m1];
+      *m1_value = q->corr_output_m1[*m1];
     }    
     ret = LIBLTE_SUCCESS;
   } 
@@ -110,17 +182,22 @@ void convert_tables(struct fc_tables *fc_tables, struct sss_tables *in) {
 
   for (i = 0; i < N_SSS; i++) {
     for (j = 0; j < N_SSS; j++) {
-      __real__ fc_tables->z1[i][j] = (float) in->z1[i][j];
+      fc_tables->z1[i][j] = (float) in->z1[i][j];
+    }
+  }
+  for (i = 0; i < N_SSS; i++) {
+    for (j = 0; j < N_SSS; j++) {
+      fc_tables->s[i][j] = (float) in->s[i][j];
     }
   }
   for (i = 0; i < N_SSS; i++) {
     for (j = 0; j < N_SSS - 1; j++) {
-      __real__ fc_tables->s[i][j] = (float) in->s[i][j + 1] * in->s[i][j];
+      fc_tables->sd[i][j] = (float) in->s[i][j + 1] * in->s[i][j];
     }
   }
   for (i = 0; i < 2; i++) {
     for (j = 0; j < N_SSS; j++) {
-      __real__ fc_tables->c[i][j] = (float) in->c[i][j];
+      fc_tables->c[i][j] = (float) in->c[i][j];
     }
   }
 }

lte/phy/lib/sync/src/gen_sss.c

@@ -121,7 +121,7 @@ void generate_sss_all_tables(struct sss_tables *tables, uint32_t N_id_2) {
   generate_s_all(tables->s, s_t);
   generate_z_all(tables->z1, z_t);
   for (i = 0; i < 2; i++) {
-    generate_c(tables->c[i], c_t, N_id_2, i > 0);
+    generate_c(tables->c[i], c_t, N_id_2, i != 0);
   }
 }
 

lte/phy/lib/sync/src/pss.c

@@ -39,10 +39,9 @@
 #include "liblte/phy/utils/debug.h"
 
 
-int pss_synch_init_N_id_2(cf_t *pss_signal_freq, uint32_t N_id_2, uint32_t fft_size) {
+int pss_synch_init_N_id_2(cf_t *pss_signal_time, cf_t *pss_signal_freq, uint32_t N_id_2, uint32_t fft_size) {
   dft_plan_t plan;
   cf_t pss_signal_pad[2048];
-  cf_t pss_signal_time[PSS_LEN];
   int ret = LIBLTE_ERROR_INVALID_INPUTS;
   
   if (lte_N_id_2_isvalid(N_id_2)    && 
@@ -65,7 +64,7 @@ int pss_synch_init_N_id_2(cf_t *pss_signal_freq, uint32_t N_id_2, uint32_t fft_s
     dft_run_c(&plan, pss_signal_pad, pss_signal_freq);
 
     vec_conj_cc(pss_signal_freq, pss_signal_freq, fft_size);
-    vec_sc_prod_cfc(pss_signal_freq, 1.0/62.0, pss_signal_freq, fft_size);
+    vec_sc_prod_cfc(pss_signal_freq, 1.0/PSS_LEN, pss_signal_freq, fft_size);
 
     dft_plan_free(&plan);
     
@@ -99,16 +98,39 @@ int pss_synch_init_fft(pss_synch_t *q, uint32_t frame_size, uint32_t fft_size) {
     
     buffer_size = fft_size + frame_size + 1;
     
+    
+    if (dft_plan(&q->dftp_input, fft_size, FORWARD, COMPLEX)) {
+      fprintf(stderr, "Error creating DFT plan \n");
+      goto clean_and_exit;
+    }
+    dft_plan_set_mirror(&q->dftp_input, true);
+    dft_plan_set_dc(&q->dftp_input, true);
+
     q->tmp_input = vec_malloc(buffer_size * sizeof(cf_t));
     if (!q->tmp_input) {
       fprintf(stderr, "Error allocating memory\n");
       goto clean_and_exit;
     }
+
+    bzero(&q->tmp_input[q->frame_size], q->fft_size * sizeof(cf_t));
+
     q->conv_output = vec_malloc(buffer_size * sizeof(cf_t));
     if (!q->conv_output) {
       fprintf(stderr, "Error allocating memory\n");
       goto clean_and_exit;
     }
+    q->conv_output_avg = vec_malloc(buffer_size * sizeof(float));
+    if (!q->conv_output_avg) {
+      fprintf(stderr, "Error allocating memory\n");
+      goto clean_and_exit;
+    }
+#ifdef PSS_ACCUMULATE_ABS
+    q->conv_output_abs = vec_malloc(buffer_size * sizeof(float));
+    if (!q->conv_output_abs) {
+      fprintf(stderr, "Error allocating memory\n");
+      goto clean_and_exit;
+    }
+#endif
     for (N_id_2=0;N_id_2<3;N_id_2++) {
       q->pss_signal_freq[N_id_2] = vec_malloc(buffer_size * sizeof(cf_t));
       if (!q->pss_signal_freq[N_id_2]) {
@@ -116,10 +138,12 @@ int pss_synch_init_fft(pss_synch_t *q, uint32_t frame_size, uint32_t fft_size) {
         goto clean_and_exit;
       }
       /* The PSS is translated into the frequency domain for each N_id_2  */
-      if (pss_synch_init_N_id_2(q->pss_signal_freq[N_id_2], N_id_2, fft_size)) {
+      if (pss_synch_init_N_id_2(q->pss_signal_time[N_id_2], q->pss_signal_freq[N_id_2], N_id_2, fft_size)) {
         fprintf(stderr, "Error initiating PSS detector for N_id_2=%d fft_size=%d\n", N_id_2, fft_size);
         goto clean_and_exit;
       }      
+      bzero(&q->pss_signal_freq[N_id_2][q->fft_size], q->frame_size * sizeof(cf_t));
+
     }    
     #ifdef CONVOLUTION_FFT
     if (conv_fft_cc_init(&q->conv_fft, frame_size, fft_size)) {
@@ -128,6 +152,8 @@ int pss_synch_init_fft(pss_synch_t *q, uint32_t frame_size, uint32_t fft_size) {
     }
     #endif
     
+    pss_synch_reset(q);
+    
     ret = LIBLTE_SUCCESS;
   }
 
@@ -157,11 +183,27 @@ void pss_synch_free(pss_synch_t *q) {
     if (q->conv_output) {
       free(q->conv_output);
     }
+#ifdef PSS_ACCUMULATE_ABS
+    if (q->conv_output_abs) {
+      free(q->conv_output_abs);
+    }
+#endif
+    if (q->conv_output_avg) {
+      free(q->conv_output_avg);
+    }
 
     bzero(q, sizeof(pss_synch_t));    
   }
 }
 
+void pss_synch_reset(pss_synch_t *q) {
+#ifdef PSS_ACCUMULATE_ABS
+  uint32_t buffer_size = q->fft_size + q->frame_size + 1;
+  bzero(q->conv_output_avg, sizeof(cf_t) * buffer_size);
+#endif
+  
+}
+
 /**
  * This function calculates the Zadoff-Chu sequence.
  * @param signal Output array.
@@ -219,7 +261,8 @@ int pss_synch_set_N_id_2(pss_synch_t *q, uint32_t N_id_2) {
   }
 }
 
-/** Returns the index of the PSS correlation peak in a subframe.
+/** Performs time-domain PSS correlation. 
+ * Returns the index of the PSS correlation peak in a subframe.
  * The frame starts at corr_peak_pos-subframe_size/2.
  * The value of the correlation is stored in corr_peak_value.
  *
@@ -241,9 +284,7 @@ int pss_synch_find_pss(pss_synch_t *q, cf_t *input, float *corr_peak_value)
       return LIBLTE_ERROR;
     }
     
-    bzero(&q->pss_signal_freq[q->N_id_2][q->fft_size], q->frame_size * sizeof(cf_t));
     memcpy(q->tmp_input, input, q->frame_size * sizeof(cf_t));
-    bzero(&q->tmp_input[q->frame_size], q->fft_size * sizeof(cf_t));
 
     /* Correlate input with PSS sequence */
   #ifdef CONVOLUTION_FFT
@@ -253,16 +294,85 @@ int pss_synch_find_pss(pss_synch_t *q, cf_t *input, float *corr_peak_value)
     conv_output_len = conv_cc(input, q->pss_signal_freq[q->N_id_2], q->conv_output, q->frame_size, q->fft_size);
   #endif
 
+   
+#ifdef PSS_ACCUMULATE_ABS
+#ifdef PSS_ABS_SQUARE
+    vec_abs_square_cf(q->conv_output, q->conv_output_abs, conv_output_len-1);
+#else
+    vec_abs_cf(q->conv_output, q->conv_output_abs, conv_output_len-1);
+#endif
     /* Find maximum of the absolute value of the correlation */
-    corr_peak_pos = vec_max_abs_ci(q->conv_output, conv_output_len-1);
+    corr_peak_pos = vec_max_fi(q->conv_output_abs, conv_output_len-1);
+    
+    // Normalize correlation output
+    vec_sc_prod_fff(q->conv_output_abs, 1/q->conv_output_abs[corr_peak_pos], q->conv_output_abs, conv_output_len-1);
+    
+    vec_sum_fff(q->conv_output_abs, q->conv_output_avg, q->conv_output_avg, conv_output_len-1);
+#else
+
+#ifdef PSS_ABS_SQUARE
+    vec_abs_square_cf(q->conv_output, q->conv_output_avg, conv_output_len-1);
+#else
+    vec_abs_cf(q->conv_output, q->conv_output_avg, conv_output_len-1);
+#endif
+
+#endif
+    
+    /* Find maximum of the absolute value of the correlation */
+    corr_peak_pos = vec_max_fi(q->conv_output_avg, conv_output_len-1);
+    
+#ifdef PSS_RETURN_PSR    
+    // Find second side lobe
+    float tmp = q->conv_output_avg[corr_peak_pos];
+    q->conv_output_avg[corr_peak_pos] = 0; 
+    int side_lobe_pos = vec_max_fi(q->conv_output_avg, conv_output_len-1);
+    q->conv_output_avg[corr_peak_pos] = tmp;
+    if (corr_peak_value) {
+      *corr_peak_value = tmp/q->conv_output_avg[side_lobe_pos];
+    }
+#else
     if (corr_peak_value) {
-      *corr_peak_value = cabsf(q->conv_output[corr_peak_pos]);
+      *corr_peak_value = q->conv_output_avg[corr_peak_pos];
     }
+#endif
+
     ret = (int) corr_peak_pos;                
   } 
   return ret;
 }
 
+LIBLTE_API cf_t *tmp2; 
+
+/* Computes frequency-domain channel estimation of the PSS symbol 
+ * input signal is in the time-domain. 
+ * ce is the returned frequency-domain channel estimates. 
+ */
+int pss_synch_chest(pss_synch_t *q, cf_t *input, cf_t ce[PSS_LEN]) {
+  int ret = LIBLTE_ERROR_INVALID_INPUTS;
+  cf_t input_fft[SYMBOL_SZ_MAX];
+
+  if (q                 != NULL  && 
+      input             != NULL)
+  {
+
+    if (!lte_N_id_2_isvalid(q->N_id_2)) {
+      fprintf(stderr, "Error finding PSS peak, Must set N_id_2 first\n");
+      return LIBLTE_ERROR;
+    }
+    
+    tmp2 = input_fft; 
+    
+    /* Transform to frequency-domain */
+    dft_run_c(&q->dftp_input, input, input_fft);
+    
+    /* Compute channel estimate taking the PSS sequence as reference */
+    vec_prod_conj_ccc(q->pss_signal_time[q->N_id_2], &input_fft[(q->fft_size-PSS_LEN)/2], ce, PSS_LEN);
+      
+    ret = LIBLTE_SUCCESS;
+  }
+  return ret; 
+}
+
 /* Returns the CFO estimation given a PSS received sequence
  *
  * Source: An Efficient CFO Estimation Algorithm for the Downlink of 3GPP-LTE

lte/phy/lib/sync/src/sss.c

@@ -55,11 +55,12 @@ int sss_synch_init(sss_synch_t *q, uint32_t fft_size) {
       sss_synch_free(q);
       return LIBLTE_ERROR;
     }
+    dft_plan_set_mirror(&q->dftp_input, true);
+    dft_plan_set_dc(&q->dftp_input, true);
+
     q->fft_size = fft_size; 
     
     generate_N_id_1_table(q->N_id_1_table);
-    dft_plan_set_mirror(&q->dftp_input, true);
-    dft_plan_set_dc(&q->dftp_input, true);
     
     for (N_id_2=0;N_id_2<3;N_id_2++) {
       generate_sss_all_tables(&sss_tables, N_id_2);
@@ -81,6 +82,7 @@ int sss_synch_realloc(sss_synch_t *q, uint32_t fft_size) {
       return LIBLTE_ERROR;
     }
     dft_plan_set_mirror(&q->dftp_input, true);
+    dft_plan_set_norm(&q->dftp_input, true);
     dft_plan_set_dc(&q->dftp_input, true);
     
     q->fft_size = fft_size;

lte/phy/lib/sync/src/sync.c

@@ -214,7 +214,7 @@ int sync_sss(sync_t *q, cf_t *input, uint32_t peak_pos) {
     return LIBLTE_ERROR;
   }
       
-  sss_synch_m0m1(&q->sss, &input[sss_idx], &q->m0, &q->m0_value, &q->m1, &q->m1_value);
+  sss_synch_m0m1_diff(&q->sss, &input[sss_idx], &q->m0, &q->m0_value, &q->m1, &q->m1_value);
 
   q->sf_idx = sss_synch_subframe(q->m0, q->m1);
   ret = sss_synch_N_id_1(&q->sss, q->m0, q->m1);
@@ -315,7 +315,8 @@ int sync_find(sync_t *q, cf_t *input, uint32_t find_offset, uint32_t *peak_posit
       ret = 0;
     }
     
-    INFO("SYNC ret=%d N_id_2=%d pos=%d peak=%.2f/%.2f=%.2f mean_energy=%.2f threshold=%.2f sf_idx=%d, CFO=%.3f KHz\n",
+    INFO("SYNC ret=%d N_id_2=%d pos=%d peak=%.2f/%.2f=%.2f mean_energy=%.2f"
+         "threshold=%.2f sf_idx=%d, CFO=%.3f KHz\n",
          ret, q->N_id_2, peak_pos, peak_unnormalized*1000,energy*1000,q->peak_value, q->mean_energy*1000, 
          q->threshold, q->sf_idx, 15*q->cfo);
 
@@ -328,4 +329,5 @@ int sync_find(sync_t *q, cf_t *input, uint32_t find_offset, uint32_t *peak_posit
 
 void sync_reset(sync_t *q) {
   q->frame_cnt = 0;
+  pss_synch_reset(&q->pss);
 }

lte/phy/lib/sync/test/CMakeLists.txt

@@ -24,12 +24,21 @@
 ########################################################################
 
 LIST(FIND OPTIONAL_LIBS cuhd CUHD_FIND)
+LIST(FIND OPTIONAL_LIBS graphics GRAPHICS_FIND)
   IF(${CUHD_FIND} GREATER -1)
   ADD_EXECUTABLE(pss_usrp pss_usrp.c)
   TARGET_LINK_LIBRARIES(pss_usrp lte_phy cuhd)
 ENDIF(${CUHD_FIND} GREATER -1)
 
+IF(${GRAPHICS_FIND} EQUAL -1)
+  SET_TARGET_PROPERTIES(pss_usrp PROPERTIES COMPILE_DEFINITIONS "DISABLE_GRAPHICS")
+ELSE(${GRAPHICS_FIND} EQUAL -1)
+  target_link_libraries(pss_usrp graphics)
+ENDIF(${GRAPHICS_FIND} EQUAL -1)
+
+
 BuildMex(MEXNAME pss SOURCES pss_mex.c LIBRARIES lte_phy liblte_mex)
+BuildMex(MEXNAME sss SOURCES sss_mex.c LIBRARIES lte_phy liblte_mex)
 
 ########################################################################
 # SYNC TEST  

lte/phy/lib/sync/test/pss_mex.c

@@ -38,7 +38,6 @@ void help()
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
 
-  int i;
   lte_cell_t cell; 
   pss_synch_t pss; 
   cf_t *input_symbols;

lte/phy/lib/sync/test/pss_usrp.c

@@ -38,7 +38,16 @@
 #include "liblte/phy/phy.h"
 #include "liblte/cuhd/cuhd.h"
 
-uint32_t N_id_2 = 100;
+
+#ifndef DISABLE_GRAPHICS
+void init_plots();
+void do_plots(float *corr, float energy, uint32_t size, cf_t ce[PSS_LEN]);
+void do_plots_sss(float *corr_m0, float *corr_m1);
+#endif
+
+
+bool disable_plots = false; 
+int cell_id = -1;
 char *uhd_args="";
 float uhd_gain=40.0, uhd_freq=-1.0;
 int nof_frames = -1;
@@ -46,18 +55,23 @@ uint32_t fft_size=64;
 float threshold = 0.4; 
 
 void usage(char *prog) {
-  printf("Usage: %s [agtvnp] -f rx_frequency_hz -i N_id_2\n", prog);
+  printf("Usage: %s [adgtvnp] -f rx_frequency_hz -i cell_id\n", prog);
   printf("\t-a UHD args [Default %s]\n", uhd_args);
   printf("\t-g UHD Gain [Default %.2f dB]\n", uhd_gain);
   printf("\t-n nof_frames [Default %d]\n", nof_frames);
   printf("\t-s symbol_sz [Default %d]\n", fft_size);
   printf("\t-t threshold [Default %.2f]\n", threshold);
+#ifndef DISABLE_GRAPHICS
+  printf("\t-d disable plots [Default enabled]\n");
+#else
+  printf("\t plots are disabled. Graphics library not available\n");
+#endif
   printf("\t-v verbose\n");
 }
 
 void parse_args(int argc, char **argv) {
   int opt;
-  while ((opt = getopt(argc, argv, "agtvsfi")) != -1) {
+  while ((opt = getopt(argc, argv, "adgtvsfi")) != -1) {
     switch (opt) {
     case 'a':
       uhd_args = argv[optind];
@@ -72,7 +86,7 @@ void parse_args(int argc, char **argv) {
       threshold = atof(argv[optind]);
       break;
     case 'i':
-      N_id_2 = atoi(argv[optind]);
+      cell_id = atoi(argv[optind]);
       break;
     case 's':
       fft_size = atoi(argv[optind]);
@@ -80,6 +94,9 @@ void parse_args(int argc, char **argv) {
     case 'n':
       nof_frames = atoi(argv[optind]);
       break;
+    case 'd':
+      disable_plots = true;
+      break;
     case 'v':
       verbose++;
       break;
@@ -88,28 +105,38 @@ void parse_args(int argc, char **argv) {
       exit(-1);
     }
   }
-  if (N_id_2 > 2 || uhd_freq < 0) {
+  if (cell_id < 0 || uhd_freq < 0) {
     usage(argv[0]);
     exit(-1);
   }
 }
+  float m0_value, m1_value; 
 
 int main(int argc, char **argv) {
   cf_t *buffer; 
   int frame_cnt, n; 
   void *uhd;
   pss_synch_t pss; 
+  cfo_t cfocorr; 
+  sss_synch_t sss; 
   int32_t flen; 
   int peak_idx, last_peak;
   float peak_value; 
   float mean_peak; 
   uint32_t nof_det, nof_nodet, nof_nopeak, nof_nopeakdet;
+  cf_t ce[PSS_LEN]; 
   
   parse_args(argc, argv);
 
+  uint32_t N_id_2 = cell_id%3;
+  uint32_t N_id_1 = cell_id/3;
+  
+  if (!disable_plots)
+    init_plots();
+  
   flen = 4800*(fft_size/64);
   
-  buffer = malloc(sizeof(cf_t) * flen);
+  buffer = malloc(sizeof(cf_t) * flen * 2);
   if (!buffer) {
     perror("malloc");
     exit(-1);
@@ -124,11 +151,21 @@ int main(int argc, char **argv) {
     exit(-1);
   }
   
+  cfo_init(&cfocorr, flen); 
+ 
+  if (sss_synch_init(&sss, fft_size)) {
+    fprintf(stderr, "Error initializing SSS object\n");
+    return LIBLTE_ERROR;
+  }
+
+  sss_synch_set_N_id_2(&sss, N_id_2);
+
   printf("Opening UHD device...\n");
   if (cuhd_open(uhd_args, &uhd)) {
     fprintf(stderr, "Error opening uhd\n");
     exit(-1);
   }
+  printf("N_id_2: %d\n", N_id_2);  
   printf("Set RX rate: %.2f MHz\n", cuhd_set_rx_srate(uhd, flen*2*100) / 1000000);
   printf("Set RX gain: %.1f dB\n", cuhd_set_rx_gain(uhd, uhd_gain));
   printf("Set RX freq: %.2f MHz\n", cuhd_set_rx_freq(uhd, uhd_freq) / 1000000);
@@ -142,8 +179,14 @@ int main(int argc, char **argv) {
   frame_cnt = 0;
   last_peak = 0; 
   mean_peak = 0;
+  int peak_offset = 0;
+  float cfo; 
+  float mean_cfo = 0; 
+  uint32_t m0, m1; 
+  uint32_t sss_error1 = 0, sss_error2 = 0, sss_error3 = 0; 
+  
   while(frame_cnt < nof_frames || nof_frames == -1) {
-    n = cuhd_recv(uhd, buffer, flen, 1);
+    n = cuhd_recv(uhd, buffer, flen - peak_offset, 1);
     if (n < 0) {
       fprintf(stderr, "Error receiving samples\n");
       exit(-1);
@@ -155,28 +198,59 @@ int main(int argc, char **argv) {
       exit(-1);
     }
     
-    float y = sqrtf(crealf(vec_dot_prod_conj_ccc(&buffer[peak_idx-fft_size], 
+    mean_peak = VEC_CMA(peak_value, mean_peak, frame_cnt);
-                                                 &buffer[peak_idx-fft_size], 
-                                                 fft_size)) /
-                                                 fft_size);
-    float x = peak_value/y;
-    
-    mean_peak = VEC_CMA(x, mean_peak, frame_cnt);
     
-    if (x >= threshold) {
+    if (peak_value >= threshold) {
       nof_det++;
+        
+      if (peak_idx >= fft_size) {
+
+        // Estimate CFO 
+        cfo = pss_synch_cfo_compute(&pss, &buffer[peak_idx-fft_size]);
+        mean_cfo = VEC_CMA(cfo, mean_cfo, frame_cnt);        
+
+        // Correct CFO
+        cfo_correct(&cfocorr, buffer, buffer, -mean_cfo / fft_size);               
+
+        // Estimate channel
+        if (pss_synch_chest(&pss, &buffer[peak_idx-fft_size], ce)) {
+          fprintf(stderr, "Error computing channel estimation\n");
+          exit(-1);
+        }
+        
+        // Find SSS 
+        int sss_idx = peak_idx-flen/10+SLOT_IDX_CPNORM(5,fft_size);             
+        if (sss_idx >= 0 && sss_idx < flen-fft_size) {
+          sss_synch_m0m1_diff(&sss, &buffer[sss_idx], &m0, &m0_value, &m1, &m1_value);
+          if (sss_synch_N_id_1(&sss, m0, m1) != N_id_1) {
+            sss_error3++;            
+          }
+          sss_synch_m0m1_partial(&sss, &buffer[sss_idx], 3, NULL, &m0, &m0_value, &m1, &m1_value);
+          if (sss_synch_N_id_1(&sss, m0, m1) != N_id_1) {
+            sss_error2++;            
+          }
+          sss_synch_m0m1_partial(&sss, &buffer[sss_idx], 1, NULL, &m0, &m0_value, &m1, &m1_value);
+          if (sss_synch_N_id_1(&sss, m0, m1) != N_id_1) {
+            sss_error1++;            
+          }
+        }
+        
+        
+      } else {
+        INFO("No space for CFO computation. Frame starts at \n",peak_idx);
+      }
+      
+      if (!disable_plots && sss_synch_subframe(m0,m1) == 0)
+        do_plots_sss(sss.corr_output_m0, sss.corr_output_m1);
+      
     } else {
       nof_nodet++;
     }
 
     if (frame_cnt > 100) {
       if (abs(last_peak-peak_idx) > 10) {
-        if (x >= threshold) {
+        if (peak_value >= threshold) {
           nof_nopeakdet++;
-        } else {
-          if (nof_nodet > 0) {
-            nof_nodet--;            
-          }
         } 
         nof_nopeak++;                  
       } 
@@ -184,16 +258,23 @@ int main(int argc, char **argv) {
     
     frame_cnt++;
    
-    printf("[%5d]: Pos: %5d, En: %.4f Val: %.3f MeanVal: %.3f, Det: %.3f, No-Det: %.3f, NoPeak: %.3f, NoPeakDet: %.3f\r",
+    printf("[%5d]: Pos: %5d, PSR: %6.3f MeanPSR: %6.3f, Pdet: %6.3f, Pmiss: %.3f, "
+           "FP: %d, FPThres: %d, CFO: %+4.1f KHz SSSmiss: %.3f/%.3f/%.3f\r",
            frame_cnt, 
-           peak_idx, y, x, mean_peak,
+           peak_idx, 
+           peak_value, mean_peak,
            (float) nof_det/frame_cnt, (float) nof_nodet/frame_cnt, 
-           (float) nof_nopeak/frame_cnt, (float) nof_nopeakdet/nof_nopeak);
+           nof_nopeak, nof_nopeakdet, mean_cfo*15, 
+           (float) sss_error1/nof_det,(float) sss_error2/nof_det,(float) sss_error3/nof_det
+          );
     
     if (VERBOSE_ISINFO()) {
       printf("\n");
     }
   
+    if (!disable_plots)
+      do_plots(pss.conv_output_avg, pss.conv_output_avg[peak_idx], pss.fft_size+pss.frame_size-1, ce);
+    
     last_peak = peak_idx;
 
   }
@@ -206,3 +287,75 @@ int main(int argc, char **argv) {
   exit(0);
 }
 
+extern cf_t *tmp2;
+
+
+/**********************************************************************
+ *  Plotting Functions
+ ***********************************************************************/
+#ifndef DISABLE_GRAPHICS
+
+
+#include "liblte/graphics/plot.h"
+plot_real_t pssout;
+//plot_complex_t pce; 
+
+plot_real_t psss1;//, psss2;
+
+float tmp[100000];
+cf_t tmpce[PSS_LEN];
+
+
+void init_plots() {
+  plot_init();
+  plot_real_init(&pssout);
+  plot_real_setTitle(&pssout, "PSS xCorr");
+  plot_real_setLabels(&pssout, "Index", "Absolute value");
+  plot_real_setYAxisScale(&pssout, 0, 1);
+  
+  /*
+  plot_complex_init(&pce);
+  plot_complex_setTitle(&pce, "Channel Estimates");
+  plot_complex_setYAxisScale(&pce, Ip, -2, 2);
+  plot_complex_setYAxisScale(&pce, Q, -2, 2);
+  plot_complex_setYAxisScale(&pce, Magnitude, 0, 2);
+  plot_complex_setYAxisScale(&pce, Phase, -M_PI, M_PI);
+  */
+  
+  plot_real_init(&psss1);
+  plot_real_setTitle(&psss1, "SSS xCorr m0");
+  plot_real_setLabels(&psss1, "Index", "Absolute value");
+  plot_real_setYAxisScale(&psss1, 0, 1);
+  
+  /*
+  plot_real_init(&psss2);
+  plot_real_setTitle(&psss2, "SSS xCorr m1");
+  plot_real_setLabels(&psss2, "Index", "Absolute value");
+  plot_real_setYAxisScale(&psss2, 0, 1);
+  */
+  
+
+}
+
+void do_plots(float *corr, float energy, uint32_t size, cf_t ce[PSS_LEN]) {  
+  vec_sc_prod_fff(corr,1./energy,tmp, size);
+  plot_real_setNewData(&pssout, tmp, size);        
+  
+//  float norm = vec_avg_power_cf(ce, PSS_LEN);
+ // vec_sc_prod_cfc(ce, 1.0/sqrt(norm), tmpce, PSS_LEN);
+  
+  //plot_complex_setNewData(&pce, tmpce, PSS_LEN);
+}
+
+void do_plots_sss(float *corr_m0, float *corr_m1) {  
+  if (m0_value > 0) 
+    vec_sc_prod_fff(corr_m0,1./m0_value,corr_m0, N_SSS);
+  plot_real_setNewData(&psss1, corr_m0, N_SSS);        
+  
+//  if (m1_value > 0) 
+//    vec_sc_prod_fff(corr_m1,1./m1_value,corr_m1, N_SSS);
+//  plot_real_setNewData(&psss2, corr_m1, N_SSS);        
+}
+
+#endif
+

lte/phy/lib/sync/test/sss_mex.c

@@ -0,0 +1,119 @@
+/* 
+ * Copyright (c) 2012, Ismael Gomez-Miguelez <ismael.gomez@tsc.upc.edu>.
+ * This file is part of ALOE++ (http://flexnets.upc.edu/)
+ * 
+ * ALOE++ is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * ALOE++ is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with ALOE++.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <string.h>
+#include "liblte/phy/phy.h"
+#include "liblte/mex/mexutils.h"
+
+/** MEX function to be called from MATLAB to test the channel estimator 
+ */
+
+#define ENBCFG  prhs[0]
+#define INPUT   prhs[1]
+#define ALGO    prhs[2]
+#define NOF_INPUTS 2
+
+
+void help()
+{
+  mexErrMsgTxt
+    ("[N_id_1,sf_idx,corr_output_m0,corr_output_m1] = liblte_sss(enbConfig, inputSignal, [Algorithm])\n"
+      "\tinputSignal must be aligned to the subframe. CP length is assumed Normal.\n"
+      "\tAlgorithm is an optional parameter: Can be 'partial','diff','full'\n\n");
+}
+
+/* the gateway function */
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+
+  lte_cell_t cell; 
+  sss_synch_t sss; 
+  cf_t *input_symbols;
+  int frame_len; 
+  uint32_t m0, m1;
+  float m0_value, m1_value; 
+  char alg[64];
+
+  if (nrhs < NOF_INPUTS) {
+    help();
+    return;
+  }
+    
+  if (mexutils_read_cell(ENBCFG, &cell)) {
+    help();
+    return;
+  }
+  
+  if (nrhs > NOF_INPUTS) {
+    mxGetString(ALGO, alg, (mwSize)sizeof(alg));       
+  } else {
+    strcpy(alg, "full");
+  }
+  
+  /** Allocate input buffers */
+  frame_len = mexutils_read_cf(INPUT, &input_symbols);
+  if (frame_len < 0) {
+    mexErrMsgTxt("Error reading input symbols\n");
+    return;
+  }
+  
+  if (sss_synch_init(&sss, lte_symbol_sz(cell.nof_prb))) {
+    mexErrMsgTxt("Error initializing SSS object\n");
+    return;
+  }
+
+  sss_synch_set_N_id_2(&sss, cell.id%3);
+      
+  // Find SSS 
+  uint32_t sss_idx = SLOT_IDX_CPNORM(5,lte_symbol_sz(cell.nof_prb));     
+  if (sss_idx > frame_len) {
+    mexErrMsgTxt("Error too few samples provided.\n");
+    return;
+  }
+  //mexPrintf("SSS begins at %d/%d. Running algorithm %s\n", sss_idx, frame_len, alg);
+  if (!strcmp(alg, "partial")) {
+    sss_synch_m0m1_partial(&sss, &input_symbols[sss_idx], 3, NULL, &m0, &m0_value, &m1, &m1_value);      
+  } else if (!strcmp(alg, "diff")) {
+    sss_synch_m0m1_diff(&sss, &input_symbols[sss_idx], &m0, &m0_value, &m1, &m1_value);      
+  } else if (!strcmp(alg, "full")) {
+    sss_synch_m0m1_partial(&sss, &input_symbols[sss_idx], 1, NULL, &m0, &m0_value, &m1, &m1_value);      
+  } else {
+    mexErrMsgTxt("Unsupported algorithm type\n");
+    return;
+  }
+  
+  //mexPrintf("m0: %d, m1: %d, N_id_1: %d\n", m0, m1, sss_synch_N_id_1(&sss, m0, m1));
+  
+  if (nlhs >= 1) { 
+    plhs[0] = mxCreateDoubleScalar(sss_synch_N_id_1(&sss, m0, m1));
+  }
+  if (nlhs >= 2) {
+    plhs[1] = mxCreateDoubleScalar(sss_synch_subframe(m0, m1));
+  }
+  if (nlhs >= 3) {
+    mexutils_write_f(sss.corr_output_m0, &plhs[2], N_SSS, 1);  
+  }
+  if (nlhs >= 4) {
+    mexutils_write_f(sss.corr_output_m1, &plhs[3], N_SSS, 1);  
+  }
+  sss_synch_free(&sss);
+  free(input_symbols);
+
+  return;
+}
+

matlab/tests/get_sc.m

@@ -0,0 +1,111 @@
+function [ s, c0, c1 ] = get_sc( N_id_2 )
+
+    if (N_id_2 == 0) 
+
+        s(1,:)=[1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, ];
+        s(2,:)=[1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, ];
+        s(3,:)=[1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, ];
+        s(4,:)=[1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, ];
+        s(5,:)=[-1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, ];
+        s(6,:)=[1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, ];
+        s(7,:)=[1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, ];
+        s(8,:)=[-1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, ];
+        s(9,:)=[1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, ];
+        s(10,:)=[-1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, ];
+        s(11,:)=[-1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, ];
+        s(12,:)=[1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, ];
+        s(13,:)=[1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, ];
+        s(14,:)=[-1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, ];
+        s(15,:)=[-1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, ];
+        s(16,:)=[-1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, ];
+        s(17,:)=[-1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, ];
+        s(18,:)=[-1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, ];
+        s(19,:)=[1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, ];
+        s(20,:)=[1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, ];
+        s(21,:)=[1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, ];
+        s(22,:)=[-1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, ];
+        s(23,:)=[-1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, ];
+        s(24,:)=[1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, ];
+        s(25,:)=[-1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, ];
+        s(26,:)=[-1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, ];
+        s(27,:)=[-1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, ];
+        s(28,:)=[1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, ];
+        s(29,:)=[-1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, ];
+        s(30,:)=[1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, ];
+        s(31,:)=[-1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, ];
+        c0=[1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, ];
+        c1=[1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, ];
+    elseif (N_id_2 == 1) 
+        s(1,:)=[1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, ];
+        s(2,:)=[1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, ];
+        s(3,:)=[1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, ];
+        s(4,:)=[1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, ];
+        s(5,:)=[-1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, ];
+        s(6,:)=[1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, ];
+        s(7,:)=[1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, ];
+        s(8,:)=[-1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, ];
+        s(9,:)=[1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, ];
+        s(10,:)=[-1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, ];
+        s(11,:)=[-1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, ];
+        s(12,:)=[1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, ];
+        s(13,:)=[1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, ];
+        s(14,:)=[-1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, ];
+        s(15,:)=[-1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, ];
+        s(16,:)=[-1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, ];
+        s(17,:)=[-1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, ];
+        s(18,:)=[-1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, ];
+        s(19,:)=[1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, ];
+        s(20,:)=[1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, ];
+        s(21,:)=[1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, ];
+        s(22,:)=[-1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, ];
+        s(23,:)=[-1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, ];
+        s(24,:)=[1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, ];
+        s(25,:)=[-1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, ];
+        s(26,:)=[-1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, ];
+        s(27,:)=[-1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, ];
+        s(28,:)=[1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, ];
+        s(29,:)=[-1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, ];
+        s(30,:)=[1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, ];
+        s(31,:)=[-1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, ];
+        c0=[1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, ];
+        c1=[-1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, ];
+
+    else
+        s(1,:)=[1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, ];
+        s(2,:)=[1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, ];
+        s(3,:)=[1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, ];
+        s(4,:)=[1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, ];
+        s(5,:)=[-1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, ];
+        s(6,:)=[1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, ];
+        s(7,:)=[1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, ];
+        s(8,:)=[-1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, ];
+        s(9,:)=[1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, ];
+        s(10,:)=[-1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, ];
+        s(11,:)=[-1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, ];
+        s(12,:)=[1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, ];
+        s(13,:)=[1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, ];
+        s(14,:)=[-1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, ];
+        s(15,:)=[-1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, ];
+        s(16,:)=[-1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, ];
+        s(17,:)=[-1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, ];
+        s(18,:)=[-1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, ];
+        s(19,:)=[1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, ];
+        s(20,:)=[1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, ];
+        s(21,:)=[1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, ];
+        s(22,:)=[-1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, ];
+        s(23,:)=[-1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, ];
+        s(24,:)=[1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, ];
+        s(25,:)=[-1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, ];
+        s(26,:)=[-1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, ];
+        s(27,:)=[-1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, ];
+        s(28,:)=[1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, ];
+        s(29,:)=[-1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, ];
+        s(30,:)=[1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, ];
+        s(31,:)=[-1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, ];
+        c0=[1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, ];
+        c1=[1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, -1, ];
+
+    end
+
+end
+

matlab/tests/test_sss.m

@@ -0,0 +1,143 @@
+
+SNR_values = linspace(-6,4,10);
+Npackets = 200; 
+CFO=4/15;
+m0=7;
+m1=10;
+%m0=26;
+%m1=21;
+
+recordedWaveform = signal;
+if (~isempty(recordedWaveform))
+    Npackets = floor(length(signal)/19200)-1;
+    SNR_values = 0;
+end
+
+error = zeros(6,length(SNR_values));
+
+enb = struct('NCellID',196,'NSubframe',0,'NDLRB',6,'CellRefP',1,'CyclicPrefix','Normal','DuplexMode','FDD');
+sss=lteSSS(enb);
+
+cfg.Seed = 2;                  % Random channel seed
+cfg.NRxAnts = 1;               % 1 receive antenna
+cfg.DelayProfile = 'ETU';      % EVA delay spread
+cfg.DopplerFreq = 144;         % 120Hz Doppler frequency
+cfg.MIMOCorrelation = 'Low';   % Low (no) MIMO correlation
+cfg.NTerms = 16;               % Oscillators used in fading model
+cfg.ModelType = 'GMEDS';       % Rayleigh fading model type
+cfg.InitPhase = 'Random';      % Random initial phases
+cfg.NormalizePathGains = 'On'; % Normalize delay profile power 
+cfg.NormalizeTxAnts = 'On';    % Normalize for transmit antennas       % Initialize at time zero
+
+[s, c0, c1] = get_sc(mod(enb.NCellID,3));
+
+subframe = lteDLResourceGrid(enb);
+sssSym = lteSSS(enb);
+sssInd = lteSSSIndices(enb);
+subframe(sssInd) = sssSym;
+N_id_1 = floor(enb.NCellID/3);
+
+[txWaveform,info] = lteOFDMModulate(enb,subframe); 
+cfg.SamplingRate = info.SamplingRate;
+fftSize = info.Nfft;
+
+
+addpath('../../debug/lte/phy/lib/sync/test')
+
+for snr_idx=1:length(SNR_values)
+    SNRdB = SNR_values(snr_idx);
+    for i=1:Npackets
+        %% Noise Addition
+        SNR = 10^(SNRdB/10);    % Linear SNR  
+
+        if (isempty(recordedWaveform))
+            cfg.InitTime = i*(10^-3);              
+            [rxWaveform, info]= lteFadingChannel(cfg,txWaveform);
+            rxWaveform = txWaveform; 
+
+            % Add CFO
+            freq = CFO/double(fftSize);
+            rxWaveform = rxWaveform.*exp(1i*2*pi*freq*(1:length(txWaveform))');
+
+            N0 = 1/(sqrt(2.0*enb.CellRefP*double(fftSize))*SNR);
+            noise = N0*complex(randn(size(rxWaveform)), randn(size(rxWaveform)));  % Generate noise
+
+            rxWaveform = rxWaveform + noise;
+        else
+            rxWaveform = recordedWaveform(i*19200+1:(i+1)*19200); 
+        end
+        
+        offset = lteDLFrameOffset(enb,rxWaveform);
+        rxWaveform = [rxWaveform(1+offset:end,:); zeros(offset,1)];
+
+        subframe_rx = lteOFDMDemodulate(enb,rxWaveform,1);
+        
+        sss_rx = subframe_rx(lteSSSIndices(enb));
+        sss0=sss_rx(1:2:end);
+        sss1=sss_rx(2:2:end);
+
+        beta0=sss0.*c0';
+        beta1=sss1.*c1';
+
+        corr0=zeros(31,1);
+        for m=1:31
+            corr0(m)=sum(beta0.*s(m,:)');    
+        end
+        corr0=abs(corr0).^2;
+        [m, idx]=max(corr0);
+        
+        error(1,snr_idx) = error(1,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+        M=2;
+        Nm=10;
+
+        corr2=zeros(31,1);
+        for m=1:31
+            for j=0:M
+                idx=1+j*Nm:(j+1)*Nm;
+                corr2(m)=corr2(m)+abs(sum(beta0(idx).*s(m,idx)')).^2;
+            end
+        end
+        [m, idx]=max(corr2);
+        
+        error(2,snr_idx) = error(2,snr_idx) + ((idx ~= m0 && idx ~= m1));
+
+        corr3=zeros(31,1);
+        for m=1:31
+            corr3(m)=abs(sum(beta0(2:end).*conj(beta0(1:end-1)).*transpose(s(m,2:end).*conj(s(m,1:end-1))))).^2;            
+        end
+        [m, idx]=max(corr3);
+        
+        error(3,snr_idx) = error(3,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+        % libLTE results
+        [n,sf_idx,lt_corr0,lt_corr1,lt_sss0,lt_sss1]=liblte_sss(enb,rxWaveform,'full');
+        [m, idx]=max(lt_corr0);
+        error(4,snr_idx) = error(4,snr_idx) + ((idx ~= m0 && idx ~= m1));
+
+        [n,sf_idx,lt_corr2]=liblte_sss(enb,rxWaveform,'partial');
+        [m, idx]=max(lt_corr2);
+        error(5,snr_idx) = error(5,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+        [n,sf_idx,lt_corr3]=liblte_sss(enb,rxWaveform,'diff');
+        [m, idx]=max(lt_corr3);
+        error(6,snr_idx) = error(6,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+    end
+end
+
+if (length(SNR_values) > 1) 
+    plot(SNR_values,1-error/Npackets)
+    legend('Full','Partial','Differential','Full-lt','Partial-lt','Differential-lt')
+    grid on
+else
+    e=error/Npackets;
+    fprintf('Full (mt/lt): \t%f/%f\n',e(1),e(4));
+    fprintf('Partial (mt/lt):%f/%f\n',e(2),e(5));
+    fprintf('Diff (mt/lt): \t%f/%f\n',e(3),e(6));
+end
+
+
+
+
+