Fixed issue in PRACH generation. Added Matlab test

lte/phy/include/liblte/phy/phy.h

@@ -96,6 +96,7 @@
 #include "liblte/phy/phch/pdsch.h"
 #include "liblte/phy/phch/phich.h"
 #include "liblte/phy/phch/pusch.h"
+#include "liblte/phy/phch/prach.h"
 #include "liblte/phy/phch/ra.h"
 #include "liblte/phy/phch/regs.h"
 #include "liblte/phy/phch/sch.h"

lte/phy/lib/phch/src/prach.c

@@ -29,6 +29,7 @@
 #include <string.h>
 #include "liblte/phy/phch/prach.h"
 #include "liblte/phy/utils/debug.h"
+#include "liblte/phy/utils/vector.h"
 
 #define   N_SEQS        64    // Number of prach sequences available
 #define   N_RB_SC       12    // Number of subcarriers per resource block
@@ -201,6 +202,7 @@ int prach_gen_seqs(prach_t *p)
       }else{
         u = prach_zc_roots[p->rsi + p->N_roots];
       }
+
       for(int j=0;j<p->N_zc;j++){
         double phase = -M_PI*u*j*(j+1)/p->N_zc;
         root[j] = cexp(phase*I);
@@ -305,23 +307,23 @@ int prach_init(prach_t *p,
     }
 
     // Set up containers
-    p->prach_bins = malloc(sizeof(cf_t)*p->N_zc);
-    p->corr_spec = malloc(sizeof(cf_t)*p->N_zc);
-    p->corr = malloc(sizeof(float)*p->N_zc);
+    p->prach_bins = vec_malloc(sizeof(cf_t)*p->N_zc);
+    p->corr_spec = vec_malloc(sizeof(cf_t)*p->N_zc);
+    p->corr = vec_malloc(sizeof(float)*p->N_zc);
 
     // Set up ZC FFTS
-    p->zc_fft = (dft_plan_t*)malloc(sizeof(dft_plan_t));
+    p->zc_fft = (dft_plan_t*)vec_malloc(sizeof(dft_plan_t));
     if(dft_plan(p->zc_fft, p->N_zc, FORWARD, COMPLEX)){
       return LIBLTE_ERROR;
     }
-    dft_plan_set_mirror(p->zc_fft, true);
+    dft_plan_set_mirror(p->zc_fft, false);
     dft_plan_set_norm(p->zc_fft, true);
 
-    p->zc_ifft = (dft_plan_t*)malloc(sizeof(dft_plan_t));
+    p->zc_ifft = (dft_plan_t*)vec_malloc(sizeof(dft_plan_t));
     if(dft_plan(p->zc_ifft, p->N_zc, BACKWARD, COMPLEX)){
       return LIBLTE_ERROR;
     }
-    dft_plan_set_mirror(p->zc_ifft, true);
+    dft_plan_set_mirror(p->zc_ifft, false);
     dft_plan_set_norm(p->zc_ifft, true);
 
     // Generate our 64 sequences
@@ -341,16 +343,16 @@ int prach_init(prach_t *p,
       p->N_ifft_prach = p->N_ifft_ul * DELTA_F/DELTA_F_RA;
     }
 
-    p->ifft_in = (cf_t*)malloc(p->N_ifft_prach*sizeof(cf_t));
-    p->ifft_out = (cf_t*)malloc(p->N_ifft_prach*sizeof(cf_t));
-    p->ifft = (dft_plan_t*)malloc(sizeof(dft_plan_t));
+    p->ifft_in = (cf_t*)vec_malloc(p->N_ifft_prach*sizeof(cf_t));
+    p->ifft_out = (cf_t*)vec_malloc(p->N_ifft_prach*sizeof(cf_t));
+    p->ifft = (dft_plan_t*)vec_malloc(sizeof(dft_plan_t));
     if(dft_plan(p->ifft, p->N_ifft_prach, BACKWARD, COMPLEX)){
        return -1;
     }
     dft_plan_set_mirror(p->ifft, true);
     dft_plan_set_norm(p->ifft, true);
 
-    p->fft = (dft_plan_t*)malloc(sizeof(dft_plan_t));
+    p->fft = (dft_plan_t*)vec_malloc(sizeof(dft_plan_t));
     if(dft_plan(p->fft, p->N_ifft_prach, FORWARD, COMPLEX)){
        return -1;
     }
@@ -382,22 +384,24 @@ int prach_gen(prach_t *p,
     uint32_t K = DELTA_F/DELTA_F_RA;
     uint32_t begin = PHI + (K*k_0) + (K/2);
 
+    printf("N_zc: %d, N_cp: %d, N_seq: %d, N_ifft_prach=%d begin: %d\n", p->N_zc, p->N_cp, p->N_seq, p->N_ifft_prach, begin);
     // Map dft-precoded sequence to ifft bins
-    memset(p->ifft_in, 0, p->N_ifft_prach*sizeof(cf_t));
-    for(int i=0;i<p->N_zc;i++){
-      p->ifft_in[begin+i] = p->dft_seqs[seq_index][i];
-    }
+    memset(p->ifft_in, 0, begin*sizeof(cf_t));
+    memcpy(&p->ifft_in[begin], p->dft_seqs[seq_index], p->N_zc * sizeof(cf_t));
+    memset(&p->ifft_in[begin+p->N_zc], 0, (p->N_ifft_prach - begin - p->N_zc) * sizeof(cf_t));
+
     dft_run(p->ifft, p->ifft_in, p->ifft_out);
 
     // Copy CP into buffer
-    for(int i=0;i<p->N_cp;i++){
-      signal[i] = p->ifft_out[p->N_ifft_prach-p->N_cp+i];
-    }
+    memcpy(signal, &p->ifft_out[p->N_ifft_prach-p->N_cp], p->N_cp*sizeof(cf_t));
 
     // Copy preamble sequence into buffer
     for(int i=0;i<p->N_seq;i++){
       signal[p->N_cp+i] = p->ifft_out[i%p->N_ifft_prach];
     }
+    
+    // Normalize 
+    vec_sc_prod_cfc(signal, 1.0/sqrtf(p->N_ifft_ul), signal, (p->N_cp + p->N_seq));
 
     ret = LIBLTE_SUCCESS;
   }

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

@@ -177,3 +177,6 @@ ADD_TEST(prach_test_multi_n16 prach_test_multi -n 16)
 ADD_TEST(prach_test_multi_n8 prach_test_multi -n 8)
 ADD_TEST(prach_test_multi_n4 prach_test_multi -n 4)
  
+ 
+ BuildMex(MEXNAME prach SOURCES prach_test_mex.c LIBRARIES lte_phy liblte_mex)
+

lte/phy/lib/phch/test/prach_test.c

@@ -91,6 +91,6 @@ int main(int argc, char **argv) {
   prach_free(p);
   free(p);
 
-	printf("Done\n");
-	exit(0);
+  printf("Done\n");
+  exit(0);
 }

lte/phy/lib/phch/test/prach_test_mex.c

@@ -0,0 +1,112 @@
+/**
+ *
+ * \section COPYRIGHT
+ *
+ * Copyright 2013-2014 The libLTE Developers. See the
+ * COPYRIGHT file at the top-level directory of this distribution.
+ *
+ * \section LICENSE
+ *
+ * This file is part of the libLTE library.
+ *
+ * libLTE 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.
+ *
+ * libLTE 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.
+ *
+ * A copy of the GNU Lesser General Public License can be found in
+ * the LICENSE file in the top-level directory of this distribution
+ * and at 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 UECFG     prhs[0]
+#define PRACHCFG  prhs[1]
+#define NOF_INPUTS 2
+
+void help()
+{
+  mexErrMsgTxt
+    ("waveform = liblte_prach(ueConfig, prachConfig)\n\n");
+}
+
+/* the gateway function */
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+ 
+  if (nrhs != NOF_INPUTS) {
+    help();
+    return;
+  }
+   
+  uint32_t n_ul_rb = 0; 
+  if (mexutils_read_uint32_struct(UECFG, "NULRB", &n_ul_rb)) {
+    mexErrMsgTxt("Field NULRB not found in UE config\n");
+    return;
+  }  
+  int r = lte_symbol_sz(n_ul_rb);
+  if (r < 0) {
+    mexErrMsgTxt("Invalid NULRB\n");
+    return;
+  }
+  uint32_t N_ifft_ul = (uint32_t) r; 
+  
+  uint32_t sf_idx = 0; 
+  mexutils_read_uint32_struct(UECFG, "NSubframe", &sf_idx);
+  uint32_t nframe = 0; 
+  mexutils_read_uint32_struct(UECFG, "NFrame", &nframe);
+
+  uint32_t preamble_format = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "Format", &preamble_format);
+  uint32_t root_seq_idx = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "SeqIdx", &root_seq_idx);
+  uint32_t seq_idx = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "PreambleIdx", &seq_idx);
+  uint32_t zero_corr_zone = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "CyclicShiftIdx", &zero_corr_zone);
+  uint32_t high_speed_flag = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "HighSpeed", &high_speed_flag);
+  uint32_t timing_offset = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "TimingOffset", &timing_offset);
+  uint32_t frequency_offset = 0; 
+  mexutils_read_uint32_struct(PRACHCFG, "FreqOffset", &frequency_offset);
+
+  prach_t prach; 
+  if (prach_init(&prach, N_ifft_ul, preamble_format, root_seq_idx, high_speed_flag, zero_corr_zone)) {
+    mexErrMsgTxt("Error initiating PRACH\n");
+    return;
+  }
+
+  uint32_t nof_samples = lte_sampling_freq_hz(n_ul_rb) * 0.001;
+  
+  cf_t *signal = vec_malloc(sizeof(cf_t) * nof_samples);
+  if (!signal) {
+    mexErrMsgTxt("malloc");
+    return;
+  }
+  bzero(signal, sizeof(cf_t) * nof_samples);
+  if (prach_gen(&prach, seq_idx, frequency_offset, signal)) {
+    mexErrMsgTxt("Error generating PRACH\n");
+    return; 
+  }
+
+  
+  if (nlhs >= 0) {
+    mexutils_write_cf(signal, &plhs[0], nof_samples, 1);  
+  }
+  
+  return;
+}
+

matlab/tests/prach_test.m

@@ -0,0 +1,48 @@
+clear 
+ueConfig=struct('NULRB',6,'DuplexMode','FDD','CyclicPrefix','Normal');
+prachConfig=struct('Format',0,'SeqIdx',0,'PreambleIdx',0,'CyclicShiftIdx',0,'HighSpeed',0,'TimingOffset',0,'FreqIdx',0,'FreqOffset',0);
+
+addpath('../../debug/lte/phy/lib/phch/test')
+
+NULRB=[6 15 25 50 100];
+
+% FreqIdx, FreqOffset and TimeOffset need to be tested
+
+for n_rb=1:length(NULRB)
+    for format=0:3
+        for seqIdx=0:837
+            fprintf('format %d, seqIdx: %d\n',format,seqIdx);
+            for preambleIdx=0:63
+                for CyclicShift=0:15
+                    %for hs=0:1
+                    hs=0;
+                        ueConfig.NULRB=NULRB(n_rb);
+                        prachConfig.Format=format;
+                        prachConfig.SeqIdx=seqIdx;
+                        prachConfig.PreambleIdx=preambleIdx;
+                        prachConfig.CyclicShiftIdx=CyclicShift;
+                        prachConfig.HighSpeed=hs;
+
+                        [mat, info]=ltePRACH(ueConfig,prachConfig);
+
+                        lib=liblte_prach(ueConfig,prachConfig);
+                         err=mean(abs(mat(:)-lib(:)));
+                         if (err > 10^-3)
+                             disp(err)    
+                             error('Error!');
+                         end
+                  %  end
+                end
+            end
+        end
+    end
+end
+
+% 
+% disp(info)
+ n=1:length(mat);
+% plot(abs(double(mat)-double(lib)))
+flib=fft(lib(199:end),1536);
+fmat=fft(mat(199:end),1536);
+n=1:1536;
+plot(n,real(flib(n)),n,real(fmat(n)))

matlab/tests/pusch_test.m

@@ -9,7 +9,7 @@ addpath('../../debug/lte/phy/lib/phch/test')
  mods={'QPSK','16QAM','64QAM'};
  rvs=0;
  betas=0:3:11;
-
+ 
 for i=1:length(TBs)
     for m=1:length(mods)
         for r=1:length(rvs)