ieee802154a.c

一个UWB仿真程序包

/** \file */

#include "ieee802154a.h"

#ifdef MATLAB
#define free mxFree
#define malloc mxMalloc
#define realloc mxRealloc
#endif

#define CHUNK_SIZE 1000

void
malloc_channel_cluster_L(ieee802154a_channel_cluster* channel)
{
  int L = channel->L;

  channel->K_l = (int *)malloc(L*sizeof(int));
  channel->T_l = (double *)malloc(L*sizeof(double));
  channel->gamma_l = (double *)malloc(L*sizeof(double));
  channel->Omega_l = (double *)malloc(L*sizeof(double));
}

void
malloc_channel_cluster_K_l(ieee802154a_channel_cluster *channel,
			   int alloc_chunk_size)
{
  /* mexPrintf("malloc_cluster_KL size: %d\n",alloc_chunk_size); */
  channel->alpha_kl = (double *)malloc(alloc_chunk_size*sizeof(double));
  channel->phi_kl = (double *)malloc(alloc_chunk_size*sizeof(double));
  channel->tau_kl = (double *)malloc(alloc_chunk_size*sizeof(double));
  channel->Ealpha2_kl = (double *)malloc(alloc_chunk_size*sizeof(double));
  channel->m_kl = (double *)malloc(alloc_chunk_size*sizeof(double));

}

void
realloc_channel_cluster_K_l(ieee802154a_channel_cluster *channel,
			   int alloc_chunk_size)
{
  /* Do the reallocations */
  channel->alpha_kl = (double *)realloc(channel->alpha_kl,alloc_chunk_size*sizeof(double));
  channel->phi_kl = (double *)realloc(channel->phi_kl,alloc_chunk_size*sizeof(double));
  channel->tau_kl = (double *)realloc(channel->tau_kl,alloc_chunk_size*sizeof(double));
  channel->Ealpha2_kl = (double *)realloc(channel->Ealpha2_kl,alloc_chunk_size*sizeof(double));
  channel->m_kl = (double *)realloc(channel->m_kl,alloc_chunk_size*sizeof(double));

}

void
malloc_channel_cont_K_l(ieee802154a_channel_cont *channel_cont,
			int alloc_chunk_size)
{

  /* mexPrintf("malloc_cont_KL\n"); */
  channel_cont->alpha = (double *)malloc(alloc_chunk_size*sizeof(double));
  channel_cont->phi = (double *)malloc(alloc_chunk_size*sizeof(double));
  channel_cont->tau = (double *)malloc(alloc_chunk_size*sizeof(double));

}


void
malloc_channel_discr(ieee802154a_channel_discr* channel, int L)
{
  channel->tau = (double *)malloc(L*sizeof(double));
  channel->alpha = (double *)malloc(L*sizeof(double));
  channel->phi = (double *)malloc(L*sizeof(double));
  channel->complex_h = (gsl_complex *)malloc(L*sizeof(gsl_complex));
  channel->real_h = (double *)malloc(L*sizeof(double));
}

void
free_channel_cluster(ieee802154a_channel_cluster* channel)
{
  free(channel->K_l);
  free(channel->alpha_kl);
  free(channel->phi_kl);  
  free(channel->T_l);  
  free(channel->tau_kl);  
  free(channel->gamma_l);  
  free(channel->Omega_l);  
  free(channel->Ealpha2_kl);  
  free(channel->m_kl);  
}

void
free_channel_cont(ieee802154a_channel_cont* channel)
{
  free(channel->tau);
  free(channel->alpha);
  free(channel->phi);  
}

void
free_channel_discr(ieee802154a_channel_discr* channel)
{
  free(channel->tau);
  free(channel->alpha);
  free(channel->phi);  
  free(channel->complex_h);  
  free(channel->real_h);  
}

void
get_L(ieee802154a_channel_cluster* channel, int barL, const gsl_rng * r)
{
  double mu;

  mu = barL;
  /* We add 1 to avoid having L = 0 */
  channel->L = gsl_ran_poisson(r, mu) + 1;

#ifdef DEBUG
  fprintf(stderr,"L = %d\n",channel->L);
#endif
}


void
get_T_l(ieee802154a_channel_cluster* channel,
	ieee802154a_parameters* param, const gsl_rng * r)
{
  int i, L;
  double mu;
  double M_cluster;
  double tmp;
  
  L = channel->L;
  mu = 1/param->Lambda;

  /* Allocate space in the struct */
  malloc_channel_cluster_L(channel);

  /* Manuel: in the matlab version this is only true for los */
  channel->T_l[0] = 0;
#ifdef DEBUG
  fprintf(stderr,"T_0 = 0\n");
#endif

  for (i = 1; i < L; i++) {
      channel->T_l[i] = channel->T_l[i-1] + gsl_ran_exponential (r,mu);
#ifdef DEBUG
      fprintf(stderr,"T_%d = %.12f\n",i,channel->T_l[i]);
#endif
    }

  for (i = 0; i < L; i++) {
      channel->gamma_l[i] = param->k_gamma*channel->T_l[i] + param->gamma_0;

      M_cluster = gsl_ran_gaussian(r, param->sigma_cluster);
      tmp = 10*log10(exp(-1*(channel->T_l[i]/param->Gamma))) + M_cluster;
      channel->Omega_l[i] = pow(10,tmp*0.1);

#ifdef DEBUG
      fprintf(stderr,"gamma_%d = %.12f\n",i,channel->gamma_l[i]);
      fprintf(stderr,"Omega_%d = %.12f\n",i,channel->Omega_l[i]);
#endif
    }
}

void
get_phi_kl(ieee802154a_channel_cluster* channel, const gsl_rng * r)
{
  int k,l,last;

  /* fprintf(stderr,"get_phi_kl: L = %d\n",channel->L); */

  last = 0;
  for (l=0; l<channel->L; l++) {
    /* fprintf(stderr,"get_phi_kl: K_l = %d\n",channel->K_l[l]); */
    for (k=0; k<channel->K_l[l]; k++) {
      channel->phi_kl[last+k] = gsl_ran_flat (r,0,2*M_PI);
    }
    last = last + channel->K_l[l];
  }
}

void
get_alpha_kl(ieee802154a_channel_cluster* channel,
	     ieee802154a_parameters* param, const gsl_rng * r)
{
  int i,k,l,last;
  double mu_m_tau,sigma_m_tau;

  last = 0;
  for (l=0; l<channel->L; l++) {
    /* fprintf(stderr,"get_alpha_kl: K_l = %d\n",channel->K_l[l]); */

    for (k=0; k<channel->K_l[l]; k++) {
      i = last + k;
      /*
       * The strongest path m-factor does not seem to be implemented,
       * so comment out everything related to it for now.
       */
      if (k == 0 && param->tildem_0 != NA) {
	channel->m_kl[i] = param->tildem_0;
      } else {
	mu_m_tau = param->m_0 - param->k_m * channel->tau_kl[i];
	sigma_m_tau = param->hatm_0 - param->hatk_m * channel->tau_kl[i];
	channel->m_kl[i] = gsl_ran_lognormal (r, mu_m_tau, sigma_m_tau);
      }

      /*
       * We are now ready to compute the amplitude which is
       * nakagami(m,Ealpha2) -> sqrt(gamma(m,Ealphai2/m))
       */

      channel->alpha_kl[i] =
	sqrt(gsl_ran_gamma(r,channel->m_kl[i],channel->Ealpha2_kl[i] / 
			   channel->m_kl[i]));
    }
    last = last + channel->K_l[l];
  }
}

double
get_mean_power(ieee802154a_channel_cluster* channel,
	       ieee802154a_parameters* param, int l, int i)
{

  double mean_power;

  mean_power =
    channel->Omega_l[l] *
    exp(-channel->tau_kl[i]/channel->gamma_l[l]) /
    (channel->gamma_l[l]*
     ((1-param->beta)*param->lambda_1 + param->beta * param->lambda_2 + 1));

  return mean_power;
}

void
get_tau_kl(ieee802154a_channel_cluster* channel,
	   ieee802154a_channel_cont* channel_cont,
	   ieee802154a_parameters* param,
	   const gsl_rng * r)
{

  int L, remaining_chunks, alloc_chunk_size;
  int l, k, i;
  double mu1, mu2;
  double thld_lin, peakMeanPow;
  
  
  /* Don't consider components not within thld_db dB of peak power */
  thld_lin = pow(10,-channel->thld_db*0.1);
  /* Grow table by this many indices at each reallocation */
  alloc_chunk_size = CHUNK_SIZE;
  
  L = channel->L;
  mu1 = 1/param->lambda_1;
  mu2 = 1/param->lambda_2;
  
  /* Initial allocation */
  malloc_channel_cluster_K_l(channel, alloc_chunk_size);
  channel->K_l[0] = alloc_chunk_size;
  
  /*initialize others to 0*/
  for (l=1;l<L;l++){
    channel->K_l[l] = 0;
  }
  
  i = 0; /* Counter to keep track of the overall index into the tables */
  for (l = 0; l < L; l++) {
    k = 0;

    while (1) {
      /* Check whether we have to make more space in the tables */
      if (k >= channel->K_l[l]) {
	int new_size = alloc_chunk_size + i; /* Grow the array */
	/* mexPrintf("realloc i=%d, k=%d, l=%d   ",i,k,l); */
        realloc_channel_cluster_K_l(channel, new_size);  
	channel->K_l[l] += alloc_chunk_size; /* Update K_l array */
      }
          
      /* Draw tau */
      if (k == 0) {
	/* First one is equal to cluster arrival time */
	channel->tau_kl[i] = channel->T_l[l];
      } else {
	/*
	 * Draw it from mixed exponential
	 * unifRnd = gsl_ran_flat (r, 0.0, 1.0); 
	 * unifRnd = gsl_rng_uniform (r); 
	 */
	if (param->beta > gsl_rng_uniform (r)) {
	  channel->tau_kl[i] = channel->tau_kl[i-1] +
	    gsl_ran_exponential (r,mu1);   
	} else {
	  channel->tau_kl[i] = channel->tau_kl[i-1] +
	    gsl_ran_exponential (r,mu2);
	}
      }
          
      /*
	fprintf(stderr,"Omega = %.12f, beta = %.12f, lambda_1 = %.12f,lambda_2 = %.12f, channel->gamma_l[l] = %.12f, tau = %.12f\n",
	channel->Omega_l[l],param->beta,param->lambda_1,param->lambda_2,
	channel->gamma_l[l],channel->tau_kl[i]);
      */

      /* Calculate mean power */
      channel->Ealpha2_kl[i] = get_mean_power(channel,param,l,i);

      /* fprintf(stderr,"meanPow = %.12f\n",channel->Ealpha2_kl[i]); */
      if (k==0) {
	peakMeanPow = channel->Ealpha2_kl[i];
      }

      /* Check whether we are below threshold */
      if (channel->Ealpha2_kl[i] < thld_lin * peakMeanPow) {
	/* Transfer remaining space in table to next cluster */
	remaining_chunks = channel->K_l[l] - k;
	/* mexPrintf("rem_chunks = %d\n",remaining_chunks); */

	if (remaining_chunks > 0) {
	  channel->K_l[l] -= remaining_chunks;

	  if (l != L-1) {
	    /* We are not at the last iteration (cluster) */
	    channel->K_l[l+1] += remaining_chunks;
	  } else {
	    /*
	     * We are at the last iteration and there is space left
	     * in the table -> shrink the whole thing
	     *
	     * i stores the total length at this moment
	     */
	    realloc_channel_cluster_K_l(channel, i); 
	  }

	  /* mexPrintf("k_l %d, %d  ", channel->K_l[l],channel->K_l[l+1]); */
	}
	break;
      }
      /* Increase counters */
      k++; i++;

    } /* End loop over k */
    /* fprintf(stderr,"no chunks = %d\n",channel->K_l[l]); */

  } /* End loop over l */

  /* Draw channel path amplitudes */
  get_alpha_kl(channel,param,r);

  /* Draw phi random uniform */
  get_phi_kl(channel, r);

  /* Save the total length of the arrays */
  malloc_channel_cont_K_l(channel_cont, i);
  channel_cont->length = i;



#ifdef DEBUG
  /* Check the whole thing */
  i = 0;
  for (l=0; l<L; l++) {

    fprintf(stderr,"T_l[%d] = %.3f\n",l,channel->T_l[l]);

    for (k=0; k<channel->K_l[l]; k++) {
      fprintf(stderr,
	      "tau_kl[%d] = %.3f, phi_kl[%d] = %.3f, alpha_kl[%d] = %.3f\n",
	      i,channel->tau_kl[i],i,channel->phi_kl[i],i,channel->alpha_kl[i]);
      i++;    
    }
  }

  fprintf(stderr,"--------------------------------------------\n");
#endif

    
  /* sort the whole thing according to tau  */
  sort_channel(channel, channel_cont);
    
}

void
sort_channel(ieee802154a_channel_cluster* channel,
	     ieee802154a_channel_cont* channel_cont)
{
  int i;
  unsigned int p[channel_cont->length]; /* this guy holds the permutation */

  gsl_sort_index (p, channel->tau_kl, 1, channel_cont->length);
    
  /* check the whole thing again, this time sorted */
  for (i=0; i<channel_cont->length; i++) {

    channel_cont->tau[i] = channel->tau_kl[p[i]];
    channel_cont->phi[i] = channel->phi_kl[p[i]];
    channel_cont->alpha[i] = channel->alpha_kl[p[i]];

#ifdef DEBUG
    fprintf(stderr,"tau[%d] = %.3f, phi[%d] = %.3f, alpha[%d] = %.3f\n",
	    i,channel_cont->tau[i],i,channel_cont->phi[i],i,
	    channel_cont->alpha[i]);
#endif

  }

}

void
convert_discrete(ieee802154a_channel_cont* channel_cont,
		 ieee802154a_channel_discr* channel_discr, double fs)
{

  int ch_len, nb_bins;
  int i,j;
  double ts, max_tau;
  double temp_a;

  /* tau is in ns, so get sampling period in ns */
  ts = 1.0 / fs;
  ch_len = channel_cont->length;
  
  /* Determine length of arrays */
  max_tau = channel_cont->tau[ch_len-1];
  nb_bins = ceil(max_tau / ts) + 1;

  malloc_channel_discr(channel_discr, nb_bins);
  channel_discr->length = nb_bins; 
  
  /* Create time axis plus complex channel coefs from sinc */
  for (i=0; i<nb_bins; i++) {
    channel_discr->tau[i] = i*ts;
    channel_discr->complex_h[i] = gsl_complex_rect(0.0,0.0);
    for (j=0; j<ch_len; j++) {
      temp_a = 
	channel_cont->alpha[j] *
	gsl_sf_sinc((channel_discr->tau[i] - channel_cont->tau[j]) / ts);
      channel_discr->complex_h[i] = 
	gsl_complex_add(
			channel_discr->complex_h[i],
			gsl_complex_polar(temp_a,channel_cont->phi[j]));
    }
    channel_discr->alpha[i] = gsl_complex_abs(channel_discr->complex_h[i]);
    channel_discr->phi[i] = gsl_complex_arg(channel_discr->complex_h[i]);
    channel_discr->real_h[i] = channel_discr->complex_h[i].dat[0]; 

#ifdef DEBUG
    fprintf(stderr,"%d: tau=%.5f, alpha=%.5f, phi=%.5f, real=%.9f, check=%.9f\n",
	    i,channel_discr->tau[i],channel_discr->alpha[i],
	    channel_discr->phi[i],channel_discr->real_h[i],
	    channel_discr->alpha[i]*cos(channel_discr->phi[i]));
#endif

  }
}


void
print_channel(ieee802154a_channel_cont* channel_cont,
	      ieee802154a_channel_discr* channel_discr, char * file)
{
  int i;
  FILE* fp;

  if (strcmp(file,"") == 0) {
    fp = stdout;
  } else {
    if ((fp=fopen(file,"w")) == NULL) {
      fprintf(stderr,"Cannot open file %s",file);
      exit(-1);
    }
  }

  for (i=0; i<channel_cont->length; i++) {
    fprintf(fp,"%d: tau=%.5f, alpha=%.5f, phi=%.5f\n",
	    i,channel_cont->tau[i],channel_cont->alpha[i],
	    channel_cont->phi[i]);
  }

  for (i=0; i<channel_discr->length; i++) {
    fprintf(fp,"%d: tau=%.5f, alpha=%.5f, phi=%.5f, real=%.9f\n",
	    i,channel_discr->tau[i],channel_discr->alpha[i],
	    channel_discr->phi[i],channel_discr->real_h[i]);
  }

  if (fp != stdout)
    fclose(fp);
}