Extend nr_dlsim to simulate performance over MIMO Rayleigh1 fading for MIMO test

openair1/SIMULATION/NR_PHY/dlsim.c

@@ -297,7 +297,7 @@ int main(int argc, char **argv)
 
   //uint8_t frame_mod4,num_pdcch_symbols = 0;
 
-  SCM_t channel_model=AWGN;//Rayleigh1_anticorr;
+  SCM_t channel_model=Rayleigh1;//AWGN Rayleigh1 Rayleigh1_anticorr;
 
   NB_UE_INST = 1;
   //double pbch_sinr;
@@ -753,7 +753,7 @@ int main(int argc, char **argv)
   gNB2UE = new_channel_desc_scm(n_tx,
                                 n_rx,
                                 channel_model,
- 				fs,
+                                fs/1e6,//sampling frequency in MHz
 				bw,
 				30e-9,
                                 0,
@@ -1030,7 +1030,6 @@ int main(int argc, char **argv)
           if (n_trials==1) printf("txlev[%d] = %d (%f dB) txlev_sum %d\n",aa,txlev[aa],10*log10((double)txlev[aa]),txlev_sum);
         }
         
-        
         for (i=(frame_parms->get_samples_slot_timestamp(slot,frame_parms,0)); 
              i<(frame_parms->get_samples_slot_timestamp(slot+1,frame_parms,0)); 
              i++) {
@@ -1040,27 +1039,43 @@ int main(int argc, char **argv)
             s_im[aa][i] = ((double)(((short *)txdata[aa]))[(i<<1)+1]);
           }
         }
+
         double ts = 1.0/(frame_parms->subcarrier_spacing * frame_parms->ofdm_symbol_size); 
-        //AWGN
+        //Compute AWGN variance
         sigma2_dB = 10 * log10((double)txlev_sum * ((double)UE->frame_parms.ofdm_symbol_size/(12*rel15->rbSize))) - SNR;
         sigma2    = pow(10, sigma2_dB/10);
         if (n_trials==1) printf("sigma2 %f (%f dB), txlev_sum %f (factor %f)\n",sigma2,sigma2_dB,10*log10((double)txlev_sum),(double)(double)UE->frame_parms.ofdm_symbol_size/(12*rel15->rbSize));
+
         for (aa=0; aa<n_rx; aa++) {
           bzero(r_re[aa],frame_length_complex_samples*sizeof(double));
           bzero(r_im[aa],frame_length_complex_samples*sizeof(double));
         }
         
+        // Apply MIMO Channel
+        if (channel_model != AWGN) multipath_tv_channel(gNB2UE,
+                             s_re,
+                             s_im,
+                             r_re,
+                             r_im,
+                             frame_length_complex_samples,
+                             0);
+
         for (i=frame_parms->get_samples_slot_timestamp(slot,frame_parms,0); 
              i<frame_parms->get_samples_slot_timestamp(slot+1,frame_parms,0);
              i++) {
 
           for (int aa_rx=0; aa_rx<n_rx; aa_rx++) {
-            // Apply MIMO Channel and then add noise
-            // sum up signals from different Tx antennas
-            for (aa=0; aa<n_tx; aa++) {
-              r_re[aa_rx][i] += s_re[aa][i];
-              r_im[aa_rx][i] += s_im[aa][i];
+
+            if (channel_model == AWGN) {
+              // sum up signals from different Tx antennas
+              r_re[aa_rx][i] = 0;
+              r_im[aa_rx][i] = 0;
+              for (aa=0; aa<n_tx; aa++) {
+                r_re[aa_rx][i] += s_re[aa][i];
+                r_im[aa_rx][i] += s_im[aa][i];
+              }
             }
+            // Add Gaussian noise
             ((short*) UE->common_vars.rxdata[aa_rx])[2*i]   = (short) ((r_re[aa_rx][i] + sqrt(sigma2/2)*gaussdouble(0.0,1.0)));
             ((short*) UE->common_vars.rxdata[aa_rx])[2*i+1] = (short) ((r_im[aa_rx][i] + sqrt(sigma2/2)*gaussdouble(0.0,1.0)));
             /* Add phase noise if enabled */
@@ -1246,6 +1261,8 @@ int main(int argc, char **argv)
     }
   }*/
 
+  free_channel_desc_scm(gNB2UE);
+
   for (i = 0; i < n_tx; i++) {
     free(s_re[i]);
     free(s_im[i]);

openair1/SIMULATION/TOOLS/multipath_tv_channel.c

@@ -48,9 +48,9 @@ void multipath_tv_channel(channel_desc_t *desc,
   printf("[TV CHANNEL] keep = %d : path_loss = %g (%f), nb_rx %d, nb_tx %d, dd %d, len %d max_doppler %g\n",keep_channel,path_loss,desc->path_loss_dB,desc->nb_rx,desc->nb_tx,dd,desc->channel_length,
          desc->max_Doppler);
 #endif
-  tx = (double complex **)malloc(desc->nb_tx*sizeof(double complex));
-  rx = (double complex **)malloc(desc->nb_rx*sizeof(double complex));
-  H_t= (double complex ** *) malloc(desc->nb_tx*desc->nb_rx*sizeof(double complex **));
+  tx = (double complex **)malloc(desc->nb_tx*sizeof(double complex *));
+  rx = (double complex **)malloc(desc->nb_rx*sizeof(double complex *));
+  H_t= (double complex ***) malloc(desc->nb_tx*desc->nb_rx*sizeof(double complex **));
   //  tv_H_t = (double complex *) malloc(length*sizeof(double complex));
   rx_temp= (double complex *) calloc(length,sizeof(double complex));
 
@@ -63,10 +63,10 @@ void multipath_tv_channel(channel_desc_t *desc,
   }
 
   for(i=0; i<desc->nb_tx*desc->nb_rx; i++) {
-    H_t[i] = (double complex **) malloc(length*sizeof(double complex *));
+    H_t[i] = (double complex **)malloc(desc->nb_taps*sizeof(double complex *));
 
-    for(j=0; j<length; j++) {
-      H_t[i][j] = (double complex *) calloc (desc->nb_taps,sizeof(double complex));
+    for(j=0; j<desc->nb_taps; j++) {
+      H_t[i][j] = (double complex *)calloc(length,sizeof(double complex));
     }
   }
 
@@ -116,7 +116,7 @@ void multipath_tv_channel(channel_desc_t *desc,
   free(rx);
 
   for(i=0; i<desc->nb_rx*desc->nb_tx; i++) {
-    for(j=0; j<length; j++) {
+    for(j=0; j<desc->nb_taps; j++) {
       free(H_t[i][j]);
     }
 
@@ -131,14 +131,29 @@ void multipath_tv_channel(channel_desc_t *desc,
 void tv_channel(channel_desc_t *desc,double complex ***H,uint32_t length){
 
   int i,j,p,l,k;
-  double *alpha,*phi_rad,pi=acos(-1),*w_Hz;
+  double *alpha,***phi_rad,pi=acos(-1),***w_Hz;
   alpha = (double *)calloc(desc->nb_paths,sizeof(double));
-  phi_rad = (double *)calloc(desc->nb_paths,sizeof(double));
-  w_Hz = (double *)calloc(desc->nb_paths,sizeof(double));
+  phi_rad = (double ***)malloc(desc->nb_rx*desc->nb_tx*sizeof(double **));
+  w_Hz = (double ***)malloc(desc->nb_rx*desc->nb_tx*sizeof(double **));
 
-  for(i=0; i<desc->nb_paths; i++) {
-    w_Hz[i]=desc->max_Doppler*cos(frand_a_b(0,2*pi));
-    phi_rad[i]=frand_a_b(0,2*pi);
+  for(i=0; i<desc->nb_tx*desc->nb_rx; i++) {
+    phi_rad[i]   = (double **) malloc(desc->nb_taps*sizeof(double *));
+    w_Hz[i]      = (double **) malloc(desc->nb_taps*sizeof(double *));
+    for(j=0; j<desc->nb_taps; j++) {
+      phi_rad[i][j]   = (double *) malloc(desc->nb_paths*sizeof(double));
+      w_Hz[i][j]      = (double *) malloc(desc->nb_paths*sizeof(double));
+    }
+  }
+
+  for(i=0; i<desc->nb_tx*desc->nb_rx; i++) {
+    for (j = 0; j<desc->nb_taps; j++) {
+      for(k=0; k<desc->nb_paths; k++) {
+        w_Hz[i][j][k] = desc->max_Doppler*cos(frand_a_b(0,2*M_PI));
+        phi_rad[i][j][k] = frand_a_b(0,2*M_PI);
+        //printf("w_hz[%d][%d][%d]=f_d*cos(theta) = %f\n",i,j,k,w_Hz[i][j][k]);
+        //printf("phi_rad[%d][%d][%d] = %f\n",i,j,k,phi_rad[i][j][k]);
+        }
+    }
   }
 
   if(desc->ricean_factor == 1) {
@@ -153,44 +168,34 @@ void tv_channel(channel_desc_t *desc,double complex ***H,uint32_t length){
     }
   }
 
-  /*
-  // This is the code when we only consider a SISO case
-  for(i=0;i<length;i++)
-  {
-  for(j=0;j<desc->nb_taps;j++)
-     {
-    for(p=0;p<desc->nb_paths;p++)
-       {
-         H[i][j] += sqrt(desc->amps[j]/2)*alpha[p]*cexp(-I*(2*pi*w_Hz[p]*i*(1/(desc->sampling_rate*1e6))+phi_rad[p]));
-       }
-       }
-   }
-  for(j=0;j<desc->nb_paths;j++)
-   {
-  phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->sampling_rate)+phi_rad[j],2*pi);
-   }
-  */
-
-  // if MIMO
+  // SISO or MIMO
   for (i=0; i<desc->nb_rx; i++) {
     for(j=0; j<desc->nb_tx; j++) {
-      for(k=0; k<length; k++) {
-        for(l=0; l<desc->nb_taps; l++) {
-          H[i+(j*desc->nb_rx)][k][l] = 0;
-
+      for(k=0; k<desc->nb_taps; k++) {
+        for(l=0; l<length; l++) {
           for(p=0; p<desc->nb_paths; p++) {
-            H[i+(j*desc->nb_rx)][k][l] += sqrt(desc->amps[l]/2)*alpha[p]*cexp(I*(2*pi*w_Hz[p]*k*(1/(desc->sampling_rate*1e6))+phi_rad[p]));
+            H[i+(j*desc->nb_rx)][k][l] += sqrt(desc->amps[k])*alpha[p]*cexp(I*(2*pi*w_Hz[i+(j*desc->nb_rx)][k][p]*l*(1/(desc->sampling_rate*1e6))+phi_rad[i+(j*desc->nb_rx)][k][p]));
           }
         }
-      }
-
-      for(j=0; j<desc->nb_paths; j++) {
-        phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->sampling_rate)+phi_rad[j],2*pi);
+        //printf("H[tx%d][rx%d][k%d] = %f+j%f \n",j,i,k,creal(H[i+(j*desc->nb_rx)][k][0]),cimag(H[i+(j*desc->nb_rx)][k][0]));
       }
     }
   }
-
+  //accumlate the phase
+  /*for(k=0; k<desc->nb_taps; k++) {
+   * for(j=0; j<desc->nb_paths; j++) {
+   * desc->random_phase[k][j] = fmod(2*pi*w_Hz[k][j]*(length-1)*(1/(desc->sampling_rate*1e6))+phi_rad[k][j],2*pi);
+   * }
+   * }*/
   free(alpha);
+  for(i=0; i<desc->nb_rx*desc->nb_tx; i++) {
+    for (j=0; j<desc->nb_taps; j++) {
+      free(w_Hz[i][j]);
+      free(phi_rad[i][j]);
+    }
+    free(w_Hz[i]);
+    free(phi_rad[i]);
+  }
   free(w_Hz);
   free(phi_rad);
 }
@@ -203,7 +208,7 @@ void tv_conv(double complex **h, double complex *x, double complex *y, uint32_t
   for(i=0; i<((int)nb_samples-dd); i++) {
     for(j=0; j<nb_taps; j++) {
       if(i>j)
-        y[i+dd] += creal(h[i][j])*creal(x[i-j])-cimag(h[i][j])*cimag(x[i-j]) + I*(creal(h[i][j])*cimag(x[i-j])+cimag(h[i][j])*creal(x[i-j]));
+        y[i+dd] += creal(h[j][i])*creal(x[i-j])-cimag(h[j][i])*cimag(x[i-j]) + I*(creal(h[j][i])*cimag(x[i-j])+cimag(h[j][i])*creal(x[i-j]));
     }
   }
 }

openair1/SIMULATION/TOOLS/random_channel.c

@@ -61,7 +61,7 @@ static telnetshell_vardef_t channelmod_vardef[] = {
 static double snr_dB=25;
 static double sinr_dB=0;
 static unsigned int max_chan;
-static channel_desc_t**  defined_channels;
+static channel_desc_t **defined_channels;
 void fill_channel_desc(channel_desc_t *chan_desc,
                        uint8_t nb_tx,
                        uint8_t nb_rx,
@@ -1185,11 +1185,11 @@ channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
                         1);
       break;
 
-    case Rayleigh1:
+    case Rayleigh1://MIMO Test uses Rayleigh1
       nb_taps = 1;
       Td = 0;
       channel_length = 1;
-      ricean_factor = 1;
+      ricean_factor = 0.0;
       aoa = .03;
       maxDoppler = 0;
       fill_channel_desc(chan_desc,nb_tx,
@@ -1632,7 +1632,7 @@ channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
      
 void free_channel_desc_scm(channel_desc_t *ch) {
   // Must be made cleanly, a lot of leaks...
-  defined_channels[ch->chan_idx]=NULL; 
+  if (max_chan != 0) defined_channels[ch->chan_idx]=NULL;
   if(ch->free_flags&CHANMODEL_FREE_AMPS)
     free(ch->amps);
   for (int i = 0; i<ch->nb_tx*ch->nb_rx; i++) {