Added matlab models for prach Pdet a Pfa tests. Exposed PRACH detection scaling factor to the public API

matlab/tests/prach_detect_test.m

@@ -1,8 +1,10 @@
 %% PRACH Detection Conformance Test
 clear
 
-numSubframes = 50;  % Number of subframes frames to simulate at each SNR
-SNRdB = linspace(-12,-4,7);  % SNR points to simulate
+detect_factor=5;
+
+numSubframes = 100;  % Number of subframes frames to simulate at each SNR
+SNRdB = linspace(-12,-6,8);  % SNR points to simulate
 foffset = 270.0;                        % Frequency offset in Hertz
 
 %% UE Configuration
@@ -75,30 +77,25 @@ for nSNR = 1:length(SNRdB)
 
         % Channel modeling
         chcfg.InitTime = (nsf-1)/1000;
-        %[rxwave, fadinginfo] = lteFadingChannel(chcfg, ...
-        %                        [txwave; zeros(25, 1)]);
+        [rxwave, fadinginfo] = lteFadingChannel(chcfg, ...
+                                [txwave; zeros(25, 1)]);
 
-        rxwave=txwave;
-        
         % Add noise
         noise = N*complex(randn(size(rxwave)), randn(size(rxwave)));            
         rxwave = rxwave + noise;            
 
         % Remove the implementation delay of the channel modeling
-        %rxwave = rxwave((fadinginfo.ChannelFilterDelay + 1):1920, :);  
+        rxwave = rxwave((fadinginfo.ChannelFilterDelay + 1):1920, :);  
 
         % Apply frequency offset
-        %t = ((0:size(rxwave, 1)-1)/chcfg.SamplingRate).';
-        %rxwave = rxwave .* repmat(exp(1i*2*pi*foffset*t), ...
-        %    1, size(rxwave, 2));
+        t = ((0:size(rxwave, 1)-1)/chcfg.SamplingRate).';
+        rxwave = rxwave .* repmat(exp(1i*2*pi*foffset*t), ...
+            1, size(rxwave, 2));
 
         % PRACH detection for all cell preamble indices
         [detected, offsets] = ltePRACHDetect(ue, prach, rxwave, (0:63).');
         
-        txwave_srs = srslte_prach(ue, prach);
-        plot(abs(txwave-txwave_srs))
-        
-        [detected_srs] = srslte_prach_detect(ue, prach, rxwave);
+        [detected_srs] = srslte_prach_detect(ue, prach, rxwave, detect_factor);
         
         
         % Test for preamble detection

matlab/tests/prach_falsealarm_test.m

@@ -0,0 +1,83 @@
+%% PRACH False Alarm Probability Conformance Test
+clear
+
+detect_factor=5;
+
+ue.NULRB = 6;                   % Number of resource blocks
+ue.DuplexMode = 'FDD';          % FDD duplexing mode
+
+%% PRACH Configuration
+prach.Format = 0;               % Preamble format
+prach.SeqIdx = 2;               % Logical root sequence index
+prach.CyclicShiftIdx = 1;       % Cyclic shift configuration index
+prach.HighSpeed = 0;            % High speed flag
+prach.FreqOffset = 0;           % Use default frequency resource index
+prach.PreambleIdx = [];         % Empty since no preamble is transmitted
+
+%% Establish PRACH Generator Output Length for this Configuration
+info = ltePRACHInfo(ue, prach);
+nSamples = info.SamplingRate*info.TotSubframes*0.001;
+
+%% Loop for Detection in Each Subframe 
+numTrials = 2000;
+falseCount = 0;         % Initialize false detection counter
+falseCount_srs = 0;         % Initialize false detection counter
+rng('default');         % Random number generator to default state
+
+runningP=zeros(1, numTrials);
+runningP_srs=zeros(1, numTrials);
+for nt = 1:numTrials 
+    
+    % Create noise
+    noise = complex(randn(nSamples, 1), randn(nSamples, 1));
+
+    % Attempt detection for all cell preamble indices (0...63)
+    [detected,offset] = ltePRACHDetect(ue, prach, noise, 0:63);
+    [detected_srs] = srslte_prach_detect(ue, prach, noise, detect_factor);
+
+    % Record false alarm
+    if (~isempty(detected))
+        falseCount = falseCount+1;        
+    end
+
+    if (~isempty(detected_srs))
+        falseCount_srs = falseCount_srs+1;        
+    end
+
+    % Calculate running false alarm probability
+    runningP(nt) = falseCount/nt*100;
+    runningP_srs(nt) = falseCount_srs/nt*100;
+    
+    % Plot information about false alarm (if applicable)
+    if (~isempty(detected))
+        plot(nt,runningP(nt),'ro','LineWidth',2,'MarkerSize',7);
+        hold on;
+        text(nt,runningP(nt), sprintf(['Preamble index = %d' ...
+            '  \nTiming offset = %0.2f samples  '],detected,offset), ...
+            'HorizontalAlignment','right');
+    end
+
+    if (~isempty(detected_srs))
+        plot(nt,runningP(nt),'mx','LineWidth',2,'MarkerSize',7);
+        hold on;
+        text(nt,runningP(nt), sprintf(['SRS index = %d' ...
+            '  \nTiming offset = %0.2f samples  '],detected,offset), ...
+            'HorizontalAlignment','right');
+    end
+
+end
+
+%% Compute Final False Alarm Probability
+
+P = falseCount / numTrials;
+P_srs = falseCount_srs / numTrials;
+plot(1:numTrials,runningP,'b','LineWidth',2);
+plot(1:numTrials,runningP_srs,'k','LineWidth',2);
+hold off
+axis([0 numTrials+1 -0.1 0.2]);
+xlabel('Trials');
+ylabel('Running false alarm probability (%)');
+title('PRACH False Alarm Detection Probability');
+
+fprintf('\nFalse alarm probability = %0.4f%% - %.04f%%\n',P*100, P_srs*100);
+

srslte/include/srslte/phch/prach.h

@@ -85,6 +85,9 @@ typedef struct SRSLTE_API {
   // ZC-sequence FFT and IFFT
   srslte_dft_plan_t *zc_fft;
   srslte_dft_plan_t *zc_ifft;
+  
+  cf_t *signal_fft; 
+  float detect_factor; 
 
 } srslte_prach_t;
 
@@ -128,6 +131,9 @@ SRSLTE_API int srslte_prach_detect(srslte_prach_t *p,
                                    uint32_t *indices,
                                    uint32_t *ind_len);
 
+SRSLTE_API void srslte_prach_set_detect_factor(srslte_prach_t *p, 
+                                               float factor); 
+
 SRSLTE_API int srslte_prach_free(srslte_prach_t *p);
 
 SRSLTE_API int srslte_prach_print_seqs(srslte_prach_t *p);

srslte/lib/phch/prach.c

@@ -337,6 +337,7 @@ int srslte_prach_init(srslte_prach_t *p,
     p->rsi = root_seq_index;
     p->hs = high_speed_flag;
     p->zczc = zero_corr_zone_config;
+    p->detect_factor = PRACH_DETECT_FACTOR; 
     
     // Determine N_zc and N_cs
     if(4 == preamble_format){
@@ -403,6 +404,13 @@ int srslte_prach_init(srslte_prach_t *p,
       fprintf(stderr, "Error creating DFT plan\n");
       return -1;
     }
+    
+    p->signal_fft = srslte_vec_malloc(sizeof(cf_t)*p->N_ifft_prach);
+    if (!p->signal_fft) {
+      fprintf(stderr, "Error allocating memory\n");
+      return -1; 
+    }
+    
     srslte_dft_plan_set_mirror(p->fft, true);
     srslte_dft_plan_set_norm(p->fft, true);
 
@@ -455,6 +463,9 @@ int srslte_prach_gen(srslte_prach_t *p,
   return ret;
 }
 
+void srslte_prach_set_detect_factor(srslte_prach_t *p, float ratio) {
+  p->detect_factor = ratio; 
+}
 
 int srslte_prach_detect(srslte_prach_t *p,
                  uint32_t freq_offset,
@@ -475,7 +486,7 @@ int srslte_prach_detect(srslte_prach_t *p,
     }
 
     // FFT incoming signal
-    srslte_dft_run(p->fft, signal, signal);
+    srslte_dft_run(p->fft, signal, p->signal_fft);
 
     memset(p->prach_bins, 0, sizeof(cf_t)*p->N_zc);
     *n_indices = 0;
@@ -487,7 +498,7 @@ int srslte_prach_detect(srslte_prach_t *p,
     uint32_t begin = PHI + (K*k_0) + (K/2);
 
     for(int i=0;i<p->N_zc;i++){
-      p->prach_bins[i] = signal[begin+i];
+      p->prach_bins[i] = p->signal_fft[begin+i];
     }
 
     for(int i=0;i<p->N_roots;i++){
@@ -497,20 +508,18 @@ int srslte_prach_detect(srslte_prach_t *p,
       float corr_ave = 0;
 
       cf_t *root_spec = p->dft_seqs[p->root_seqs_idx[i]];
-
-      for(int j=0;j<p->N_zc;j++){
-        p->corr_spec[j] = p->prach_bins[j]*conjf(root_spec[j]);
-      }
+     
+      srslte_vec_prod_conj_ccc(p->prach_bins, root_spec, p->corr_spec, p->N_zc);
 
       srslte_dft_run(p->zc_ifft, p->corr_spec, p->corr_spec);
-
+      
+      srslte_vec_abs_cf(p->corr_spec, p->corr, p->N_zc);
+      
       float norm = sqrtf(p->N_zc);
-      for(int j=0;j<p->N_zc;j++){
-        p->corr[j] = cabsf(p->corr_spec[j])/norm;
-        corr_ave += p->corr[j];
-      }
-      corr_ave /= p->N_zc;
-
+      srslte_vec_sc_prod_fff(p->corr, 1.0/norm, p->corr, p->N_zc);
+      
+      corr_ave = srslte_vec_acc_ff(p->corr, p->N_zc)/p->N_zc;
+      
       uint32_t winsize = 0;
       if(p->N_cs != 0){
         winsize = p->N_cs;
@@ -527,7 +536,7 @@ int srslte_prach_detect(srslte_prach_t *p,
             corr_max = p->corr[k];
           }
         }
-        if(corr_max > PRACH_DETECT_FACTOR*corr_ave){
+        if(corr_max > p->detect_factor*corr_ave){
           indices[*n_indices] = (i*n_wins)+j;
           (*n_indices)++;
         }
@@ -554,6 +563,10 @@ int srslte_prach_free(srslte_prach_t *p) {
   srslte_dft_plan_free(p->zc_ifft);
   free(p->zc_ifft);
 
+  if (p->signal_fft) {
+    free(p->signal_fft); 
+  }
+  
   bzero(p, sizeof(srslte_prach_t));
 
   return 0;

srslte/lib/phch/test/prach_detect_test_mex.c

@@ -47,7 +47,7 @@ void help()
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
  
-  if (nrhs != NOF_INPUTS) {
+  if (nrhs < NOF_INPUTS) {
     help();
     return;
   }
@@ -93,6 +93,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   
   uint32_t preambles[64]; 
   uint32_t nof_detected = 0; 
+  
+  if (nrhs > NOF_INPUTS) {
+    float factor = mxGetScalar(prhs[NOF_INPUTS]);
+    srslte_prach_set_detect_factor(&prach, factor);
+  }
 
   if (srslte_prach_detect(&prach, frequency_offset, &input_signal[prach.N_cp], nof_samples, preambles, &nof_detected)) {    
     mexErrMsgTxt("Error detecting PRACH\n");

srslte/lib/phch/test/prach_test.c

@@ -31,14 +31,15 @@
 #include <unistd.h>
 #include <math.h>
 #include <time.h>
+#include <sys/time.h>
 #include <complex.h>
 
-#include "srslte/phch/prach.h"
+#include "srslte/srslte.h"
 
 #define MAX_LEN  70176
 
 
-uint32_t N_ifft_ul        = 128;
+uint32_t N_ifft_ul        = 1536;
 uint32_t preamble_format  = 0;
 uint32_t root_seq_idx     = 0;
 uint32_t zero_corr_zone   = 1;
@@ -109,8 +110,9 @@ int main(int argc, char **argv) {
     uint32_t prach_len = p->N_seq;
     if(preamble_format == 2 || preamble_format == 3)
       prach_len /= 2;
+    
     srslte_prach_detect(p, 0, &preamble[p->N_cp], prach_len, indices, &n_indices);
-
+    
     if(n_indices != 1 || indices[0] != seq_index)
       return -1;
   }