IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
channel.cpp
Go to the documentation of this file.
1 
30 #include <itpp/comm/channel.h>
31 #include <itpp/base/math/error.h>
33 #include <itpp/base/bessel.h>
34 #include <itpp/base/matfunc.h>
35 #include <itpp/base/specmat.h>
36 #include <itpp/signal/resampling.h>
37 #include <itpp/signal/transforms.h>
38 #include <itpp/signal/window.h>
39 #include <itpp/base/math/min_max.h>
40 #include <itpp/stat/misc_stat.h>
41 
42 
43 namespace itpp
44 {
45 
46 
47 // --------------------------------------------------------------------------
48 // Fading_Generator class
49 // --------------------------------------------------------------------------
50 
52 {
53  // no default LOS component
54  set_LOS_power(0.0);
55 }
56 
57 void Fading_Generator::set_LOS_power(double relative_power)
58 {
59  it_assert(relative_power >= 0.0,
60  "Fading_Generator::set_LOS_power(): Relative_power can not be negative");
61  los_power = relative_power;
62  los_diffuse = std::sqrt(1.0 / (1.0 + los_power));
64 }
65 
67 {
68  it_warning("Fading_Generator::set_LOS_doppler(): This function has no effect on this kind of generator");
69 }
70 
72 {
73  it_warning("Fading_Generator::set_time_offset(): This function has no effect on this kind of generator");
74 }
75 
77 {
78  it_warning("Fading_Generator::set_norm_doppler(): This function has no effect on this kind of generator");
79 }
80 
82 {
83  it_warning("Fading_Generator::set_filter_length(): This function has no effect on this kind of generator");
84 }
85 
87 {
88  it_warning("Fading_Generator::set_doppler_spectrum(): This function has no effect on this kind of generator");
89 }
90 
92 {
93  it_warning("Fading_Generator::set_no_frequencies(): This function has no effect on this kind of generator");
94 }
95 
97 {
98  it_warning("Fading_Generator::set_rice_method(): This function has no effect on this kind of generator");
99 }
100 
102 {
103  it_warning("Fading_Generator::get_LOS_doppler(): This function has no effect on this kind of generator");
104  return 0;
105 }
106 
108 {
109  it_warning("Fading_Generator::get_time_offset(): This function has no effect on this kind of generator");
110  return 0;
111 }
112 
114 {
115  it_warning("Fading_Generator::get_filter_length(): This function has no effect on this kind of generator");
116  return 0;
117 }
118 
120 {
121  it_warning("Fading_Generator::get_norm_doppler(): This function has no effect on this kind of generator");
122  return 0;
123 }
124 
126 {
127  it_warning("Fading_Generator::get_doppler_spectrum(): This function has no effect on this kind of generator");
128  return Jakes;
129 }
130 
132 {
133  it_warning("Fading_Generator::get_no_frequencies(): This function has no effect on this kind of generator");
134  return 0;
135 }
136 
138 {
139  it_warning("Fading_Generator::get_rice_method(): This function has no effect on this kind of generator");
140  return MEDS;
141 }
142 
144 {
145  it_warning("Fading_Generator::shift_time_offset(): This function has no effect on this kind of generator");
146 }
147 
148 cvec Fading_Generator::generate(int no_samples)
149 {
150  cvec output;
151  this->generate(no_samples, output);
152  return output;
153 }
154 
155 
156 // --------------------------------------------------------------------------
157 // Independent_Fading_Generator class
158 // --------------------------------------------------------------------------
159 
160 void Independent_Fading_Generator::generate(int no_samples, cvec& output)
161 {
162  output.set_size(no_samples, false);
163  if (los_power > 0.0) {
164  for (int i = 0; i < no_samples; ++i) {
165  output(i) = los_diffuse * randn_c() + los_direct;
166  }
167  }
168  else {
169  output = randn_c(no_samples);
170  }
171 }
172 
173 
174 // --------------------------------------------------------------------------
175 // Static_Fading_Generator class
176 // --------------------------------------------------------------------------
177 
179 {
181  if (los_power > 0.0) {
184  }
185  init_flag = true;
186 }
187 
188 void Static_Fading_Generator::generate(int no_samples, cvec& output)
189 {
190  if (init_flag == false)
191  init();
192 
193  output.set_size(no_samples, false);
194  output = static_sample;
195 }
196 
197 
198 // --------------------------------------------------------------------------
199 // Correlated_Fading_Generator class
200 // --------------------------------------------------------------------------
201 
203  Fading_Generator(), los_dopp(0.7), time_offset(0.0)
204 {
205  set_norm_doppler(norm_doppler);
206 }
207 
209 {
210  it_assert((norm_doppler > 0) && (norm_doppler <= 1.0),
211  "Correlated_Fading_Generator: Normalized Doppler out of range");
212  n_dopp = norm_doppler;
213  init_flag = false;
214 }
215 
216 void Correlated_Fading_Generator::set_LOS_doppler(double relative_doppler)
217 {
218  it_assert((relative_doppler >= 0) && (relative_doppler <= 1.0),
219  "Correlated_Fading_Generator::set_LOS_doppler(): Relative Doppler out of range");
220  los_dopp = relative_doppler;
221 }
222 
224 {
225  time_offset = static_cast<double>(offset);
226 }
227 
229 {
230  time_offset += static_cast<double>(no_samples);
231 }
232 
233 void Correlated_Fading_Generator::add_LOS(int idx, std::complex<double>& sample)
234 {
235  double tmp_arg = m_2pi * los_dopp * n_dopp * (idx + time_offset);
236  sample *= los_diffuse;
237  sample += los_direct * std::complex<double>(std::cos(tmp_arg),
238  std::sin(tmp_arg));
239 }
240 
241 
242 // --------------------------------------------------------------------------
243 // Rice_Fading_Generator class
244 // --------------------------------------------------------------------------
245 
248  int no_freq, RICE_METHOD method) :
249  Correlated_Fading_Generator(norm_doppler)
250 {
251  set_doppler_spectrum(spectrum);
252  set_no_frequencies(no_freq);
253  set_rice_method(method);
254 }
255 
257 {
259  init_flag = false;
260 }
261 
263 {
264  it_assert(no_freq >= 7,
265  "Rice_Fading_Generator::set_no_frequencies(): Too low number of Doppler frequencies");
266  Ni = no_freq;
267  init_flag = false;
268 }
269 
271 {
272  // check if this method works for the given spectrum
273  rice_method = method;
274  init_flag = false;
275 }
276 
278 {
279  switch (rice_method) {
280  case MEDS: // Method of Exact Doppler Spread (MEDS)
281  init_MEDS();
282  break;
283  default:
284  it_error("Rice_Fading_Generator::init(): Wrong Rice method for this fading generator");
285  };
286 
287  init_flag = true; // generator ready to use
288 }
289 
290 void Rice_Fading_Generator::generate(int no_samples, cvec &output)
291 {
292  if (init_flag == false)
293  init();
294 
295  output.set_size(no_samples, false);
296 
297  switch (dopp_spectrum) {
298  case Jakes: {
299  double tmp_re, tmp_im;
300  if (los_power > 0.0) { // LOS component exists
301  for (int i = 0; i < no_samples; i++) {
302  tmp_re = sum(elem_mult(c1, cos(m_2pi * f1 * n_dopp * (i + time_offset) + th1)));
303  tmp_im = sum(elem_mult(c2, cos(m_2pi * f2 * n_dopp * (i + time_offset) + th2)));
304  output(i) = std::complex<double>(tmp_re, tmp_im);
305  add_LOS(i, output(i));
306  }
307  }
308  else {
309  for (int i = 0; i < no_samples; i++) {
310  tmp_re = sum(elem_mult(c1, cos(m_2pi * f1 * n_dopp * (i + time_offset) + th1)));
311  tmp_im = sum(elem_mult(c2, cos(m_2pi * f2 * n_dopp * (i + time_offset) + th2)));
312  output(i) = std::complex<double>(tmp_re, tmp_im);
313  }
314  }
315  break;
316  }
317  case GaussI:
318  case GaussII: {
319  double tmp;
320  for (int i = 0; i < no_samples; i++) {
321  tmp = m_2pi * n_dopp * (i + time_offset);
322  output(i) = sum(elem_mult(c1, cos(f1 * tmp + th1)))
323  * std::complex<double>(std::cos(f01 * tmp), -std::sin(f01 * tmp))
324  + sum(elem_mult(c2, cos(f2 * tmp + th2)))
325  * std::complex<double>(std::cos(f02 * tmp), -std::sin(f02 * tmp));
326  }
327  break;
328  }
329  }
330 
331  time_offset += no_samples;
332 }
333 
335 {
336  vec n;
337  double sgm_0_2;
338 
339  switch (dopp_spectrum) {
340  case Jakes:
341  n = linspace(1, Ni, Ni);
342  f1 = sin(pi / (2 * Ni) * (n - 0.5));
343  c1 = std::sqrt(1.0 / Ni) * ones(Ni);
344  th1 = randu(Ni) * 2 * pi;
345  n = linspace(1, Ni + 1, Ni + 1);
346  f2 = sin(pi / (2 * (Ni + 1)) * (n - 0.5));
347  c2 = std::sqrt(1.0 / (Ni + 1)) * ones(Ni + 1);
348  th2 = randu(Ni + 1) * 2 * pi;
349  f01 = f02 = 0;
350  break;
351  case GaussI:
352  n = linspace(1, Ni, Ni);
353  sgm_0_2 = 5.0 / 6.0;
354  c1 = std::sqrt(sgm_0_2 * 2.0 / Ni) * ones(Ni);
355  f1 = std::sqrt(2.0) * 0.05 * erfinv((2 * n - 1) / (2 * Ni));
356  th1 = randu(Ni) * 2 * pi;
357  sgm_0_2 = 1.0 / 6.0;
358  c2 = std::sqrt(sgm_0_2 * 2.0 / Ni) * ones(Ni);
359  f2 = std::sqrt(2.0) * 0.1 * erfinv((2 * n - 1) / (2 * Ni));
360  th2 = randu(Ni) * 2 * pi;
361  f01 = 0.8;
362  f02 = -0.4;
363  break;
364  case GaussII:
365  n = linspace(1, Ni, Ni);
366  sgm_0_2 = std::sqrt(10.0) / (std::sqrt(10.0) + 0.15);
367  c1 = std::sqrt(sgm_0_2 * 2.0 / Ni) * ones(Ni);
368  f1 = std::sqrt(2.0) * 0.1 * erfinv((2 * n - 1) / (2 * Ni));
369  th1 = randu(Ni) * 2 * pi;
370  sgm_0_2 = 0.15 / (std::sqrt(10.0) + 0.15);
371  c2 = std::sqrt(sgm_0_2 * 2.0 / Ni) * ones(Ni);
372  f2 = std::sqrt(2.0) * 0.15 * erfinv((2 * n - 1) / (2 * Ni));
373  th2 = randu(Ni) * 2 * pi;
374  f01 = -0.7;
375  f02 = 0.4;
376  break;
377  default:
378  it_error("Rice_Fading_Generator::init_MEDS(): Wrong spectrum method for this fading generator");
379  };
380 }
381 
382 
383 // --------------------------------------------------------------------------
384 // FIR_Fading_Generator class methods
385 // --------------------------------------------------------------------------
386 
388  int filter_length) :
389  Correlated_Fading_Generator(norm_doppler)
390 {
391  set_filter_length(filter_length);
392 }
393 
395 {
396  it_assert(filter_length >= 50,
397  "FIR_Fading_Generator::set_filter_length(): Filter length should be at least 50");
398  fir_length = filter_length;
399  init_flag = false;
400 }
401 
403 {
404  // calculate a reasonable upsample rate so that normalized doppler is > 0.1
405  double norm_dopp = n_dopp;
406  upsample_rate = 1;
407  while (norm_dopp < 0.1) {
408  norm_dopp *= 2;
409  upsample_rate *= 2;
410  }
412 
413  // fill filter with dummy data
414  cvec dummy = fir_filter(randn_c(fir_length));
415 
416  left_overs.set_size(0, false);
417 
418  init_flag = true; // generator ready to use
419 }
420 
421 void FIR_Fading_Generator::generate(int no_samples, cvec &output)
422 {
423  if (init_flag == false)
424  init();
425 
426  // calculate number of samples before upsampling
427  int no_upsamples = ceil_i(static_cast<double>(no_samples - left_overs.size())
428  / upsample_rate) + 1;
429 
430  // should make a smarter interpolation here!!!
431  lininterp(fir_filter(randn_c(no_upsamples)), upsample_rate, output);
432  output = concat(left_overs, output); // add left-overs from previous filtering
433  left_overs = output.right(output.size() - no_samples); // save left-overs for next round of filtering
434  output.set_size(no_samples, true);
435 
436  if (los_power > 0.0) { // LOS component exist
437  for (int i = 0; i < no_samples; i++) {
438  add_LOS(i, output(i));
439  }
440  }
441 
442  time_offset += no_samples;
443 }
444 
445 vec FIR_Fading_Generator::Jakes_filter(double norm_dopp, int order)
446 {
447  int L = order / 2;
448  vec x_pos(L), x_neg(L), x(2*L + 1), h(2*L + 1);
449  for (int i = 1; i <= L; i++) {
450  x_pos(i - 1) = besselj(0.25, m_2pi * norm_dopp * i) / std::pow(i, 0.25);
451  // / std::sqrt(std::sqrt(static_cast<double>(i)));
452  }
453  double x0 = 1.468813 * std::pow(norm_dopp, 0.25); // std::sqrt(std::sqrt(norm_dopp));
454  x_neg = reverse(x_pos);
455  x = concat(concat(x_neg, x0), x_pos);
456  h = elem_mult(hamming(2 * L + 1), x);
457  h /= norm(h);
458  return h;
459 }
460 
461 
462 // --------------------------------------------------------------------------
463 // IFFT_Fading_Generator class methods
464 // --------------------------------------------------------------------------
465 
466 void IFFT_Fading_Generator::generate(int no_samples, cvec &output)
467 {
468  if (init_flag == false)
469  init();
470 
471  generate_Jakes(no_samples, output);
472 
473  if (los_power > 0.0) { // LOS component exist
474  for (int i = 0; i < no_samples; i++) {
475  add_LOS(i, output(i));
476  }
477  }
478 
479  time_offset += no_samples;
480 }
481 
482 void IFFT_Fading_Generator::generate_Jakes(int no_samples, cvec &output)
483 {
484  int Nfft = pow2i(levels2bits(no_samples));
485  double df = 1.0 / Nfft;
486  int noisesamp = ceil_i(n_dopp / df);
487  int no_upsample = 1;
488 
489  while (noisesamp <= 10) { // if too few samples, increase the FFT size
490  Nfft *= 2;
491  no_upsample *= 2;
492  df = 1.0 / Nfft;
493  noisesamp = ceil_i(n_dopp / df);
494  it_assert(no_upsample < 128,
495  "IFFT_Fading_Generator::generate_Jakes(): Too low normalized doppler or too small blocks of data. Results in an inefficient algorithm with lots of zero-padding");
496  }
497 
498  vec Fpos = linspace(0, 0.5, Nfft / 2 + 1);
499  vec F = concat(Fpos, reverse(-Fpos(1, Nfft / 2 - 1)));
500  vec S = zeros(Nfft);
501 
502  for (int i = 0; i < F.size(); i++) {
503  if (std::fabs(F(i)) < n_dopp)
504  S(i) = std::sqrt(1.5 / (pi * n_dopp * std::sqrt(1 - std::pow(F(i) / n_dopp, 2))));
505  else if (std::fabs(F(i)) == n_dopp)
506  S(i) = 1000000;
507  }
508 
509  S /= norm(S, 2);
510  S *= Nfft;
511 
512  cvec x = zeros_c(Nfft);
513 
514  for (int i = 0; i < noisesamp; ++i) {
515  x(i) = S(i) * randn_c();
516  x(Nfft - 1 - i) = S(Nfft - 1 - i) * randn_c();
517  }
518 
519  x = ifft(x);
520 
521  output = x.mid(0, no_samples);
522 }
523 
524 
525 // --------------------------------------------------------------------------
526 // Channel_Specification class methods
527 // --------------------------------------------------------------------------
528 
530  const vec &delay_prof)
531 {
532  set_channel_profile(avg_power_dB, delay_prof);
533 }
534 
536 {
537  set_channel_profile(profile);
538 }
539 
540 void Channel_Specification::set_channel_profile(const vec &avg_power_dB, const vec &delay_prof)
541 {
542  it_assert(min(delay_prof) == 0,
543  "Channel_Specification::set_channel_profile(): Minimum relative delay must be 0");
544  it_assert(avg_power_dB.size() == delay_prof.size(),
545  "Channel_Specification::set_channel_profile(): Power and delay vectors must be of equal length");
546  it_assert(delay_prof(0) == 0,
547  "Channel_Specification::set_channel_profile(): First tap must be at zero delay");
548  for (int i = 1; i < delay_prof.size(); i++) {
549  it_assert(delay_prof(i) > delay_prof(i - 1),
550  "Channel_Specification::set_channel_profile(): Delays should be sorted and unique");
551  }
552 
553  N_taps = delay_prof.size();
554  a_prof_dB = avg_power_dB;
555  d_prof = delay_prof;
556 
557  // set doppler spectrum to Jakes per default
559  tap_doppler_spectrum = Jakes;
560 
561  // set LOS parameters to zeros per default
562  set_LOS(zeros(N_taps));
563 }
564 
566 {
567  switch (profile) {
568  // -------------- ITU Channel models -----------------
569  case ITU_Vehicular_A:
570  set_channel_profile(vec("0 -1 -9 -10 -15 -20"),
571  vec("0 310 710 1090 1730 2510") * 1e-9);
572  break;
573 
574  case ITU_Vehicular_B:
575  set_channel_profile(vec("-2.5 0 -12.8 -10 -25.2 -16"),
576  vec("0 300 8900 12900 17100 20000") * 1e-9);
577  break;
578 
579  case ITU_Pedestrian_A:
580  set_channel_profile(vec("0 -9.7 -19.2 -22.8"),
581  vec("0 110 190 410") * 1e-9);
582  break;
583 
584  case ITU_Pedestrian_B:
585  set_channel_profile(vec("0 -0.9 -4.9 -8 -7.8 -23.9"),
586  vec("0 200 800 1200 2300 3700") * 1e-9);
587  break;
588 
589  // -------------- COST259 Channel models -----------------
590  case COST259_TUx:
591  set_channel_profile(vec("-5.7 -7.6 -10.1 -10.2 -10.2 -11.5 -13.4 -16.3 -16.9 -17.1 -17.4 -19 -19 -19.8 -21.5 -21.6 -22.1 -22.6 -23.5 -24.3"),
592  vec("0 217 512 514 517 674 882 1230 1287 1311 1349 1533 1535 1622 1818 1836 1884 1943 2048 2140") * 1e-9);
593  break;
594 
595  case COST259_RAx:
596  set_channel_profile(vec("-5.2 -6.4 -8.4 -9.3 -10 -13.1 -15.3 -18.5 -20.4 -22.4"),
597  vec("0 42 101 129 149 245 312 410 469 528") * 1e-9);
598  set_LOS(0, sqr(0.91 / 0.41), 0.7);
599  break;
600 
601  case COST259_HTx:
602  set_channel_profile(vec("-3.6 -8.9 -10.2 -11.5 -11.8 -12.7 -13.0 -16.2 -17.3 -17.7 -17.6 -22.7 -24.1 -25.8 -25.8 -26.2 -29 -29.9 -30 -30.7"),
603  vec("0 356 441 528 546 609 625 842 916 941 15000 16172 16492 16876 16882 16978 17615 17827 17849 18016") * 1e-9);
604  break;
605 
606  // -------------- COST207 Channel models -----------------
607  case COST207_RA:
608  set_channel_profile(vec("0 -2 -10 -20"),
609  vec("0 200 400 600") * 1e-9);
610  set_LOS(0, sqr(0.91 / 0.41), 0.7);
611  break;
612 
613  case COST207_RA6:
614  set_channel_profile(vec("0 -4 -8 -12 -16 -20"),
615  vec("0 100 200 300 400 500") * 1e-9);
616  set_LOS(0, sqr(0.91 / 0.41), 0.7);
617  break;
618 
619  case COST207_TU:
620  set_channel_profile(vec("-3 0 -2 -6 -8 -10"),
621  vec("0 200 600 1600 2400 5000") * 1e-9);
622  set_doppler_spectrum(2, GaussI);
623  set_doppler_spectrum(3, GaussI);
624  set_doppler_spectrum(4, GaussII);
625  set_doppler_spectrum(5, GaussII);
626  break;
627 
628  case COST207_TU6alt:
629  set_channel_profile(vec("-3 0 -2 -6 -8 -10"),
630  vec("0 200 500 1600 2300 5000") * 1e-9);
631  set_doppler_spectrum(3, GaussI);
632  set_doppler_spectrum(4, GaussII);
633  set_doppler_spectrum(5, GaussII);
634  break;
635 
636  case COST207_TU12:
637  set_channel_profile(vec("-4 -3 0 -2 -3 -5 -7 -5 -6 -9 -11 -10"),
638  vec("0 200 400 600 800 1200 1400 1800 2400 3000 3200 5000") * 1e-9);
639  set_doppler_spectrum(3, GaussI);
640  set_doppler_spectrum(4, GaussI);
641  set_doppler_spectrum(5, GaussI);
642  set_doppler_spectrum(6, GaussI);
643  set_doppler_spectrum(7, GaussI);
644  set_doppler_spectrum(8, GaussII);
645  set_doppler_spectrum(9, GaussII);
646  set_doppler_spectrum(10, GaussII);
647  set_doppler_spectrum(11, GaussII);
648  break;
649 
650  case COST207_TU12alt:
651  set_channel_profile(vec("-4 -3 0 -2.6 -3 -5 -7 -5 -6.5 -8.6 -11 -10"),
652  vec("0 200 400 600 800 1200 1400 1800 2400 3000 3200 5000") * 1e-9);
653  set_doppler_spectrum(4, GaussI);
654  set_doppler_spectrum(5, GaussI);
655  set_doppler_spectrum(6, GaussI);
656  set_doppler_spectrum(7, GaussI);
657  set_doppler_spectrum(8, GaussII);
658  set_doppler_spectrum(9, GaussII);
659  set_doppler_spectrum(10, GaussII);
660  set_doppler_spectrum(11, GaussII);
661  break;
662 
663  case COST207_BU:
664  set_channel_profile(vec("-3 0 -3 -5 -2 -4"),
665  vec("0 400 1000 1600 5000 6600") * 1e-9);
666  set_doppler_spectrum(2, GaussI);
667  set_doppler_spectrum(3, GaussI);
668  set_doppler_spectrum(4, GaussII);
669  set_doppler_spectrum(5, GaussII);
670  break;
671 
672  case COST207_BU6alt:
673  set_channel_profile(vec("-2.5 0 -3 -5 -2 -4"),
674  vec("0 300 1000 1600 5000 6600") * 1e-9);
675  set_doppler_spectrum(2, GaussI);
676  set_doppler_spectrum(3, GaussI);
677  set_doppler_spectrum(4, GaussII);
678  set_doppler_spectrum(5, GaussII);
679  break;
680 
681  case COST207_BU12:
682  set_channel_profile(vec("-7 -3 -1 0 -2 -6 -7 -1 -2 -7 -10 -15"),
683  vec("0 200 400 800 1600 2200 3200 5000 6000 7200 8200 10000") * 1e-9);
684  set_doppler_spectrum(3, GaussI);
685  set_doppler_spectrum(4, GaussI);
686  set_doppler_spectrum(5, GaussII);
687  set_doppler_spectrum(6, GaussII);
688  set_doppler_spectrum(7, GaussII);
689  set_doppler_spectrum(8, GaussII);
690  set_doppler_spectrum(9, GaussII);
691  set_doppler_spectrum(10, GaussII);
692  set_doppler_spectrum(11, GaussII);
693  break;
694 
695  case COST207_BU12alt:
696  set_channel_profile(vec("-7.7 -3.4 -1.3 0 -2.3 -5.6 -7.4 -1.4 -1.6 -6.7 -9.8 -15.1"),
697  vec("0 100 300 700 1600 2200 3100 5000 6000 7200 8100 10000") * 1e-9);
698  set_doppler_spectrum(3, GaussI);
699  set_doppler_spectrum(4, GaussI);
700  set_doppler_spectrum(5, GaussII);
701  set_doppler_spectrum(6, GaussII);
702  set_doppler_spectrum(7, GaussII);
703  set_doppler_spectrum(8, GaussII);
704  set_doppler_spectrum(9, GaussII);
705  set_doppler_spectrum(10, GaussII);
706  set_doppler_spectrum(11, GaussII);
707  break;
708 
709 
710  case COST207_HT:
711  set_channel_profile(vec("0 -2 -4 -7 -6 -12"),
712  vec("0 200 400 600 15000 17200") * 1e-9);
713  set_doppler_spectrum(4, GaussII);
714  set_doppler_spectrum(5, GaussII);
715  break;
716 
717  case COST207_HT6alt:
718  set_channel_profile(vec("0 -1.5 -4.5 -7.5 -8 -17.7"),
719  vec("0 100 300 500 15000 17200") * 1e-9);
720  set_doppler_spectrum(4, GaussII);
721  set_doppler_spectrum(5, GaussII);
722  break;
723 
724  case COST207_HT12:
725  set_channel_profile(vec("-10 -8 -6 -4 0 0 -4 -8 -9 -10 -12 -14"),
726  vec("0 200 400 600 800 2000 2400 15000 15200 15800 17200 20000") * 1e-9);
727  set_doppler_spectrum(3, GaussI);
728  set_doppler_spectrum(4, GaussI);
729  set_doppler_spectrum(5, GaussI);
730  set_doppler_spectrum(6, GaussII);
731  set_doppler_spectrum(7, GaussII);
732  set_doppler_spectrum(8, GaussII);
733  set_doppler_spectrum(9, GaussII);
734  set_doppler_spectrum(10, GaussII);
735  set_doppler_spectrum(11, GaussII);
736  break;
737 
738  case COST207_HT12alt:
739  set_channel_profile(vec("-10 -8 -6 -4 0 0 -4 -8 -9 -10 -12 -14"),
740  vec("0 100 300 500 700 1000 1300 15000 15200 15700 17200 20000") * 1e-9);
741  set_doppler_spectrum(4, GaussI);
742  set_doppler_spectrum(5, GaussI);
743  set_doppler_spectrum(6, GaussI);
744  set_doppler_spectrum(7, GaussII);
745  set_doppler_spectrum(8, GaussII);
746  set_doppler_spectrum(9, GaussII);
747  set_doppler_spectrum(10, GaussII);
748  set_doppler_spectrum(11, GaussII);
749  break;
750  };
751 }
752 
753 
755 {
756  for (int i = 0; i < N_taps; i++)
757  tap_doppler_spectrum(i) = tap_spectrum[i];
758 }
759 
761 {
762  tap_doppler_spectrum(tap_number) = tap_spectrum;
763 }
764 
765 void Channel_Specification::set_LOS(int tap_number, double relative_power,
766  double relative_doppler)
767 {
768  it_assert(N_taps >= 1,
769  "Channel_Specification::set_LOS(): Cannot set LOS component if not set channel profile");
770  it_assert((tap_number >= 0) && (tap_number < N_taps),
771  "Channel_Specification::set_LOS(): Tap number out of range");
772  it_assert((relative_doppler >= 0) && (relative_doppler <= 1.0),
773  "Channel_Specification::set_LOS(): Normalized Doppler out of range");
774  it_assert(relative_power >= 0.0,
775  "Channel_Specification::set_LOS(): Rice factor out of range");
776 
777  los_power.set_size(N_taps, true);
778  los_dopp.set_size(N_taps, true);
779  los_power(tap_number) = relative_power;
780  los_dopp(tap_number) = relative_doppler;
781 }
782 
783 void Channel_Specification::set_LOS(const vec& relative_power,
784  const vec& relative_doppler)
785 {
786  it_assert((relative_power.size() == N_taps),
787  "Channel_Specification::set_LOS(): Improper size of input vectors");
788 
789  if (relative_doppler.size() == 0) {
790  los_power.set_size(relative_power.size());
791  los_dopp.set_size(relative_power.size());
792  for (int i = 0; i < relative_power.size(); i++) {
793  it_assert(relative_power(i) >= 0.0,
794  "Channel_Specification::set_LOS(): Rice factor out of range");
795  los_power(i) = relative_power(i);
796  los_dopp(i) = 0.7;
797  }
798  }
799  else {
800  it_assert(relative_doppler.size() == N_taps,
801  "Channel_Specification::set_LOS(): Improper size of input vectors");
802  los_power.set_size(relative_power.size());
803  los_dopp.set_size(relative_power.size());
804  for (int i = 0; i < relative_power.size(); i++) {
805  it_assert((relative_doppler(i) >= 0) && (relative_doppler(i) <= 1.0),
806  "Channel_Specification::set_LOS(): Normalized Doppler out of range");
807  it_assert(relative_power(i) >= 0.0,
808  "Channel_Specification::set_LOS(): Rice factor out of range");
809  los_power(i) = relative_power(i);
810  los_dopp(i) = relative_doppler(i);
811  }
812  }
813 }
814 
816  vec &delay_prof) const
817 {
818  avg_power_dB = a_prof_dB;
819  delay_prof = d_prof;
820 }
821 
823 {
824  it_assert((index >= 0) && (index < N_taps),
825  "Channel_Specification::get_doppler_spectrum(): Index of of range");
826  return tap_doppler_spectrum(index);
827 }
828 
830 {
831  vec a_prof = inv_dB(a_prof_dB);
832  return (a_prof * d_prof / sum(a_prof));
833 }
834 
836 {
837  vec a_prof = inv_dB(a_prof_dB);
838  double a = a_prof * d_prof / sum(a_prof);
839  double b = a_prof * sqr(d_prof) / sum(a_prof);
840 
841  return std::sqrt(b - a*a);
842 }
843 
844 
845 // --------------------------------------------------------------------------
846 // TDL_Channel class methods
847 // --------------------------------------------------------------------------
848 
849 TDL_Channel::TDL_Channel(const vec &avg_power_dB, const ivec &delay_prof):
850  init_flag(false), n_dopp(0.0), fading_type(Independent), method(Rice_MEDS),
851  filter_length(0), nrof_freq(16), discrete_Ts(0.0)
852 {
853  set_channel_profile(avg_power_dB, delay_prof);
854 
855  // initialize LOS parameters to all zeros
856  set_LOS(zeros(delay_prof.size()));
857 
858  // initialize Doppler spectra
859  tap_doppler_spectrum.set_size(delay_prof.size());
860  tap_doppler_spectrum = Jakes;
861 }
862 
863 TDL_Channel::TDL_Channel(const Channel_Specification &channel_spec, double sampling_time):
864  init_flag(false), n_dopp(0.0), fading_type(Independent), method(Rice_MEDS),
865  filter_length(0), nrof_freq(16), discrete_Ts(sampling_time)
866 {
867  set_channel_profile(channel_spec, sampling_time);
868 
869  // set Doppler spectrum
871 }
872 
874 {
875  if (fading_gen.size() > 0) { // delete all old generators
876  for (int i = 0; i < fading_gen.size(); i++) {
877  if (fading_gen(i) != NULL) {
878  delete fading_gen(i);
879  fading_gen(i) = NULL;
880  }
881  }
882  }
883 }
884 
885 void TDL_Channel::set_channel_profile(const vec &avg_power_dB,
886  const ivec &delay_prof)
887 {
888  it_assert(min(delay_prof) == 0,
889  "TDL_Channel::set_channel_profile(): Minimum relative delay must be 0.");
890  it_assert(avg_power_dB.size() == delay_prof.size(),
891  "TDL_Channel::set_channel_profile(): Power and delay vectors must be of equal length!");
892  it_assert(delay_prof(0) == 0,
893  "TDL_Channel::set_channel_profile(): First tap must be at zero delay");
894  for (int i = 1; i < delay_prof.size(); i++) {
895  it_assert(delay_prof(i) > delay_prof(i - 1),
896  "TDL_Channel::set_channel_profile(): Delays should be sorted and unique");
897  }
898 
899  N_taps = delay_prof.size();
900  a_prof = pow(10.0, avg_power_dB / 20.0); // Convert power profile to amplitude profile
901  a_prof /= norm(a_prof); // Normalize amplitude profile
902  d_prof = delay_prof;
903 
904  // initialize Doppler spectra
906  tap_doppler_spectrum = Jakes;
907 
908  // set size of Rice parameters according to the new channel profile
909  set_LOS(zeros(N_taps));
910 
911  // changes in PDP require initialisation
912  init_flag = false;
913 }
914 
916 {
917  it_assert(no_taps >= 1, "TDL_Channel::set_channel_profile_uniform(): Minimum number of taps is 1.");
918 
919  vec avg_power_dB = zeros(no_taps);
920  ivec delay_prof(no_taps);
921  for (int i = 0; i < no_taps; i++)
922  delay_prof(i) = i;
923 
924  set_channel_profile(avg_power_dB, delay_prof);
925 }
926 
928 {
929  it_assert(no_taps >= 1, "TDL_Channel::set_channel_profile_exponential(): Minimum number of taps is 1.");
930 
931  vec avg_power_dB(no_taps);
932  ivec delay_prof(no_taps);
933  for (int i = 0; i < no_taps; i++) {
934  delay_prof(i) = i;
935  // p(i*ts) = exp(-i*ts), k = 0...no_taps-1
936  avg_power_dB(i) = dB(std::exp(static_cast<double>(-i)));
937  }
938 
939  set_channel_profile(avg_power_dB, delay_prof);
940 }
941 
942 void TDL_Channel::set_channel_profile(const Channel_Specification &channel_spec, double sampling_time)
943 {
944  vec avg_power_dB;
945  vec delay_profile;
946 
947  // set power profile and delays
948  channel_spec.get_channel_profile(avg_power_dB, delay_profile);
949  discrete_Ts = sampling_time;
950  N_taps = avg_power_dB.size();
951  a_prof = pow(10.0, avg_power_dB / 20.0); // Convert power profile to amplitude profile
952  a_prof /= norm(a_prof); // Normalize amplitude profile
953 
954  // set size of Rice parameters according to the new channel profile
955  set_LOS(channel_spec.get_LOS_power(), channel_spec.get_LOS_doppler());
956 
957  // set Doppler spectrum
959 
960  // sets discretized delay profile
961  discretize(delay_profile);
962 
963  init_flag = false;
964 }
965 
966 
968 {
969  fading_type = Correlated;
970  method = correlated_method;
971  init_flag = false;
972 }
973 
975 {
976  fading_type = fading_type_in;
977  init_flag = false;
978 }
979 
980 
981 void TDL_Channel::set_norm_doppler(double norm_doppler)
982 {
983  it_assert((norm_doppler > 0) && (norm_doppler <= 1.0),
984  "TDL_Channel::set_norm_doppler(): Normalized Doppler out of range");
985  n_dopp = norm_doppler;
986  // if normalized Doppler is set, we have correlated fading
987  fading_type = Correlated;
988  init_flag = false;
989 }
990 
991 
992 void TDL_Channel::set_LOS(const vec& relative_power, const vec& relative_doppler)
993 {
994  it_assert((relative_power.size() == N_taps),
995  "TDL_Channel::set_LOS(): Improper size of input vectors");
996 
997  if (relative_doppler.size() == 0) {
998  los_power.set_size(relative_power.size());
999  los_dopp.set_size(relative_power.size());
1000  for (int i = 0; i < relative_power.size(); i++) {
1001  it_assert(relative_power(i) >= 0.0,
1002  "TDL_Channel::set_LOS(): Rice factor out of range");
1003  los_power(i) = relative_power(i);
1004  los_dopp(i) = (relative_power(i) > 0) ? 0.7 : 0.0;
1005  }
1006  }
1007  else {
1008  it_assert(relative_doppler.size() == N_taps,
1009  "TDL_Channel::set_LOS(): Improper size of input vectors");
1010  los_power.set_size(relative_power.size());
1011  los_dopp.set_size(relative_power.size());
1012  for (int i = 0; i < relative_power.size(); i++) {
1013  it_assert((relative_doppler(i) >= 0) && (relative_doppler(i) <= 1.0),
1014  "TDL_Channel::set_LOS(): Normalized Doppler out of range");
1015  it_assert(relative_power(i) >= 0.0,
1016  "TDL_Channel::set_LOS(): Rice factor out of range");
1017  los_power(i) = relative_power(i);
1018  los_dopp(i) = relative_doppler(i);
1019  }
1020  }
1021 }
1022 
1023 void TDL_Channel::set_LOS_power(const vec& relative_power)
1024 {
1025  it_assert(relative_power.size() == N_taps,
1026  "TDL_Channel::set_LOS_power(): Improper size of input vector");
1027 
1028  los_power.set_size(relative_power.size());
1029  los_dopp.set_size(relative_power.size());
1030  for (int i = 0; i < los_power.size(); ++i) {
1031  los_power(i) = relative_power(i);
1032  los_dopp(i) = (relative_power(i) > 0) ? 0.7 : 0.0;
1033  }
1034  init_flag = false;
1035 }
1036 
1037 void TDL_Channel::set_LOS_doppler(const vec& relative_doppler)
1038 {
1039  it_assert(relative_doppler.size() == los_power.size(),
1040  "TDL_Channel::set_LOS_doppler(): Improper size of input vector");
1041 
1042  it_assert(n_dopp > 0, "TDL_Channel::set_LOS_doppler(): Normalized Doppler needs to be non zero to set the LOS Doppler in a Correlated fading generator");
1043 
1044  los_dopp.set_size(relative_doppler.size());
1045  for (int i = 0; i < relative_doppler.size(); ++i) {
1046  it_assert((relative_doppler(i) >= 0) && (relative_doppler(i) <= 1.0),
1047  "TDL_Channel::set_LOS_doppler(): Normalized Doppler out of range");
1048  los_dopp(i) = relative_doppler(i);
1049  }
1050 
1051  init_flag = false;
1052 }
1053 
1054 
1056 {
1057  it_assert(N_taps > 0, "TDL_Channel::set_doppler_spectrum(): Channel profile not defined yet");
1058 
1059  it_assert(n_dopp > 0, "TDL_Channel::set_doppler_spectrum(): Normalized Doppler needs to be non zero to set the Doppler spectrum in the Correlated Rice MEDS fading generator");
1060 
1061  if (method != Rice_MEDS)
1062  method = Rice_MEDS;
1063 
1065  for (int i = 0; i < N_taps; i++)
1066  tap_doppler_spectrum(i) = tap_spectrum[i];
1067 
1068  init_flag = false;
1069 }
1070 
1071 void TDL_Channel::set_doppler_spectrum(int tap_number, DOPPLER_SPECTRUM tap_spectrum)
1072 {
1073  it_assert((tap_number >= 0) && (tap_number < N_taps),
1074  "TDL_Channel::set_doppler_spectrum(): Improper tap number");
1075 
1076  it_assert(n_dopp > 0, "TDL_Channel::set_doppler_spectrum(): Normalized Doppler needs to be non zero to set the Doppler spectrum in the Correlated Rice MEDS fading generator");
1077 
1078  if (method != Rice_MEDS)
1079  method = Rice_MEDS;
1080 
1082  tap_doppler_spectrum(tap_number) = tap_spectrum;
1083 
1084  init_flag = false;
1085 }
1086 
1088 {
1089  it_assert(n_dopp > 0, "TDL_Channel::set_no_frequencies(): Normalized Doppler needs to be non zero to set the number of frequencies in the Correlated Rice MEDS fading generator");
1090  nrof_freq = no_freq;
1091  if (method != Rice_MEDS)
1092  method = Rice_MEDS;
1093 
1094  init_flag = false;
1095 }
1096 
1097 
1099 {
1100  it_assert(n_dopp > 0, "TDL_Channel::set_filter_length(): Normalized Doppler needs to be non zero to use the Correlated FIR fading generator");
1101 
1102  filter_length = fir_length;
1103  if (method != FIR)
1104  method = FIR;
1105 
1106  init_flag = false;
1107 }
1108 
1109 
1111 {
1112  it_assert(n_dopp > 0, "TDL_Channel::set_time_offset(): Normalized Doppler needs to be non zero to set time offset in a Correlated fading generator");
1113 
1114  if (init_flag == false)
1115  init();
1116 
1117  for (int i = 0; i < N_taps; i++) {
1118  fading_gen(i)->set_time_offset(offset);
1119  }
1120 }
1121 
1122 
1124 {
1125  it_assert(n_dopp > 0, "TDL_Channel::shift_time_offset(): Normalized Doppler needs to be non zero to shift time offset in a Correlated fading generator");
1126 
1127  if (init_flag == false)
1128  init();
1129 
1130  for (int i = 0; i < N_taps; i++) {
1131  fading_gen(i)->shift_time_offset(no_samples);
1132  }
1133 }
1134 
1135 
1136 void TDL_Channel::get_channel_profile(vec &avg_power_dB,
1137  ivec &delay_prof) const
1138 {
1139  avg_power_dB = 20 * log10(a_prof);
1140  delay_prof = d_prof;
1141 }
1142 
1144 {
1145  return (20 * log10(a_prof));
1146 }
1147 
1149 {
1150  if (fading_gen(0) != NULL)
1151  return fading_gen(0)->get_time_offset();
1152  else
1153  return -1.0;
1154 }
1155 
1157 {
1158  return (sqr(a_prof)*d_prof / sum_sqr(a_prof));
1159 }
1160 
1162 {
1163  double a = (sqr(a_prof) * d_prof / sum_sqr(a_prof));
1164  double b = (sqr(a_prof) * sqr(to_vec(d_prof)) / sum_sqr(a_prof));
1165 
1166  return (std::sqrt(b - a*a));
1167 }
1168 
1170 {
1171  it_assert(N_taps > 0, "TDL_Channel::init(): Channel profile not defined yet");
1172  it_assert(N_taps == los_power.size(),
1173  "TDL_Channel::init(): LOS profile does not mach the channel profile");
1174 
1175  if (fading_gen.size() > 0) { // delete all old generators
1176  for (int i = 0; i < fading_gen.size(); i++) {
1177  if (fading_gen(i) != NULL) {
1178  delete fading_gen(i);
1179  fading_gen(i) = NULL;
1180  }
1181  }
1182  }
1183 
1184  // create all generators and set the parameters
1185  fading_gen.set_size(N_taps, false);
1186  switch (fading_type) {
1187 
1188  case Independent:
1189  for (int i = 0; i < N_taps; ++i) {
1191  if (los_power(i) > 0)
1192  fading_gen(i)->set_LOS_power(los_power(i));
1193  fading_gen(i)->init();
1194  }
1195  break;
1196 
1197  case Static:
1198  for (int i = 0; i < N_taps; ++i) {
1200  if (los_power(i) > 0)
1201  fading_gen(i)->set_LOS_power(los_power(i));
1202  fading_gen(i)->init();
1203  }
1204  break;
1205 
1206  case Correlated:
1207  it_assert(n_dopp > 0,
1208  "TDL_Channel::init(): Correlated fading requires non zero normalized Doppler");
1209 
1210  switch (method) {
1211  case Rice_MEDS:
1212  // The third parameter is the number of sine waveforms that create the
1213  // fading process. It is increased by 2 for each tap to make taps
1214  // uncorrelated. Minimum number of waveforms set explicitly to 16.
1215  for (int i = 0; i < N_taps; ++i) {
1217  nrof_freq + 2*i, MEDS);
1218  if (los_power(i) > 0) {
1219  fading_gen(i)->set_LOS_power(los_power(i));
1220  fading_gen(i)->set_LOS_doppler(los_dopp(i));
1221  }
1222  fading_gen(i)->init();
1223  }
1224  break;
1225 
1226  case FIR:
1227  for (int i = 0; i < N_taps; ++i) {
1228  it_assert(tap_doppler_spectrum(i) == Jakes,
1229  "TDL_Channel::init(): FIR fading generator can be used with Jakes spectrum only");
1231  if (los_power(i) > 0) {
1232  fading_gen(i)->set_LOS_power(los_power(i));
1233  fading_gen(i)->set_LOS_doppler(los_dopp(i));
1234  }
1235  if (filter_length > 0)
1236  fading_gen(i)->set_filter_length(filter_length);
1237  fading_gen(i)->init();
1238  }
1239  break;
1240 
1241  case IFFT:
1242  for (int i = 0; i < N_taps; ++i) {
1243  it_assert(tap_doppler_spectrum(i) == Jakes,
1244  "TDL_Channel::init(): IFFT fading generator can be used with Jakes spectrum only");
1246  if (los_power(i) > 0) {
1247  fading_gen(i)->set_LOS_power(los_power(i));
1248  fading_gen(i)->set_LOS_doppler(los_dopp(i));
1249  }
1250  fading_gen(i)->init();
1251  }
1252  break;
1253 
1254  default:
1255  it_error("TDL_Channel::init(): No such fading generation method");
1256  }
1257  break;
1258 
1259  default:
1260  it_error("TDL_Channel::init(): No such fading type");
1261  };
1262 
1263  init_flag = true;
1264 }
1265 
1266 void TDL_Channel::generate(int no_samples, Array<cvec> &channel_coeff)
1267 {
1268  if (init_flag == false)
1269  init();
1270 
1271  channel_coeff.set_size(N_taps, false);
1272  for (int i = 0; i < N_taps; i++)
1273  channel_coeff(i) = a_prof(i) * fading_gen(i)->generate(no_samples);
1274 }
1275 
1276 void TDL_Channel::generate(int no_samples, cmat &channel_coeff)
1277 {
1278  if (init_flag == false)
1279  init();
1280 
1281  channel_coeff.set_size(no_samples, N_taps, false);
1282  for (int i = 0; i < N_taps; i++)
1283  channel_coeff.set_col(i, a_prof(i) * fading_gen(i)->generate(no_samples));
1284 }
1285 
1286 
1287 void TDL_Channel::filter_known_channel(const cvec &input, cvec &output, const Array<cvec> &channel_coeff)
1288 {
1289  int maxdelay = max(d_prof);
1290 
1291  output.set_size(input.size() + maxdelay, false);
1292  output.zeros();
1293 
1294  for (int i = 0; i < N_taps; i++)
1295  output += concat(zeros_c(d_prof(i)), elem_mult(input, channel_coeff(i)), zeros_c(maxdelay - d_prof(i)));
1296 }
1297 
1298 void TDL_Channel::filter_known_channel(const cvec &input, cvec &output, const cmat &channel_coeff)
1299 {
1300  int maxdelay = max(d_prof);
1301 
1302  output.set_size(input.size() + maxdelay, false);
1303  output.zeros();
1304 
1305  for (int i = 0; i < N_taps; i++)
1306  output += concat(zeros_c(d_prof(i)), elem_mult(input, channel_coeff.get_col(i)), zeros_c(maxdelay - d_prof(i)));
1307 }
1308 
1309 void TDL_Channel::filter(const cvec &input, cvec &output, Array<cvec> &channel_coeff)
1310 {
1311  generate(input.size(), channel_coeff);
1312  filter_known_channel(input, output, channel_coeff);
1313 }
1314 
1315 void TDL_Channel::filter(const cvec &input, cvec &output, cmat &channel_coeff)
1316 {
1317  generate(input.size(), channel_coeff);
1318  filter_known_channel(input, output, channel_coeff);
1319 }
1320 
1321 cvec TDL_Channel::filter(const cvec &input, Array<cvec> &channel_coeff)
1322 {
1323  cvec output;
1324  filter(input, output, channel_coeff);
1325  return output;
1326 }
1327 
1328 cvec TDL_Channel::filter(const cvec &input, cmat &channel_coeff)
1329 {
1330  cvec output;
1331  filter(input, output, channel_coeff);
1332  return output;
1333 }
1334 
1335 void TDL_Channel::filter(const cvec &input, cvec &output)
1336 {
1337  Array<cvec> channel_coeff;
1338  filter(input, output, channel_coeff);
1339 }
1340 
1341 cvec TDL_Channel::filter(const cvec &input)
1342 {
1343  cvec output;
1344  filter(input, output);
1345  return output;
1346 }
1347 
1348 
1349 void TDL_Channel::operator()(const cvec &input, cvec &output, Array<cvec> &channel_coeff)
1350 {
1351  filter(input, output, channel_coeff);
1352 }
1353 
1354 void TDL_Channel::operator()(const cvec &input, cvec &output, cmat &channel_coeff)
1355 {
1356  filter(input, output, channel_coeff);
1357 }
1358 
1359 
1360 cvec TDL_Channel::operator()(const cvec &input, Array<cvec> &channel_coeff)
1361 {
1362  return filter(input, channel_coeff);
1363 }
1364 
1365 cvec TDL_Channel::operator()(const cvec &input, cmat &channel_coeff)
1366 {
1367  return filter(input, channel_coeff);
1368 }
1369 
1370 cvec TDL_Channel::operator()(const cvec &input)
1371 {
1372  return filter(input);
1373 }
1374 
1375 
1376 void TDL_Channel::calc_impulse_response(const Array<cvec> &channel_coeff, Array<cvec> &impulse_response)
1377 {
1378  it_assert(init_flag == true, "calc_impulse_response: TDL_Channel is not initialized");
1379  it_assert(N_taps == channel_coeff.size(), "calc_impulse_response: number of channel taps do not match");
1380 
1381  int no_samples = channel_coeff(0).size();
1382  it_assert(no_samples > 0, "calc_impulse_response: channel_coeff must contain samples");
1383 
1384  impulse_response.set_size(no_samples);
1385 
1386  for (int i = 0; i < no_samples; i++) {
1387  impulse_response(i).set_size(d_prof(N_taps - 1) + 1, false);
1388  impulse_response(i).zeros();
1389 
1390  for (int j = 0; j < N_taps; j++)
1391  impulse_response(i)(d_prof(j)) = channel_coeff(j)(i);
1392 
1393  }
1394 }
1395 
1396 void TDL_Channel::calc_frequency_response(const Array<cvec> &channel_coeff, Array<cvec> &frequency_response, const int fft_size)
1397 {
1398  it_assert(init_flag == true, "calc_frequency_response: TDL_Channel is not initialized");
1399  it_assert(N_taps == channel_coeff.size(), "calc_frequency_response: number of channel taps do not match");
1400 
1401  int no_samples = channel_coeff(0).size();
1402  it_assert(no_samples > 0, "calc_frequency_response: channel_coeff must contain samples");
1403 
1404  frequency_response.set_size(no_samples);
1405 
1406  it_assert(fft_size > d_prof(N_taps - 1), "calc_frequency_response: fft_size must be larger than the maximum delay in samples");
1407  cvec impulse_response(fft_size);
1408 
1409  for (int i = 0; i < no_samples; i++) {
1410  impulse_response.zeros();
1411 
1412  for (int j = 0; j < N_taps; j++)
1413  impulse_response(d_prof(j)) = channel_coeff(j)(i);
1414 
1415  fft(impulse_response, frequency_response(i));
1416 
1417  }
1418 }
1419 
1420 void TDL_Channel::calc_frequency_response(const cmat &channel_coeff, cmat &frequency_response, const int fft_size)
1421 {
1422  it_assert(init_flag == true, "calc_frequency_response: TDL_Channel is not initialized");
1423  it_assert(N_taps == channel_coeff.cols(), "calc_frequency_response: number of channel taps do not match");
1424 
1425  int no_samples = channel_coeff.rows();
1426  it_assert(no_samples > 0, "calc_frequency_response: channel_coeff must contain samples");
1427 
1428  frequency_response.set_size(fft_size, no_samples, false);
1429 
1430  it_assert(fft_size > d_prof(N_taps - 1), "calc_frequency_response: fft_size must be larger than the maximum delay in samples");
1431  cvec impulse_response(fft_size);
1432  cvec freq;
1433 
1434  for (int i = 0; i < no_samples; i++) {
1435  impulse_response.zeros();
1436 
1437  for (int j = 0; j < N_taps; j++)
1438  impulse_response(d_prof(j)) = channel_coeff(i, j);
1439 
1440  fft(impulse_response, freq);
1441  frequency_response.set_col(i, freq);
1442  }
1443 }
1444 
1445 void TDL_Channel::discretize(const vec &delay_profile)
1446 {
1447  it_assert(N_taps > 0, "TDL_Channel::discretize(): No channel profile specified");
1448  it_assert(delay_profile(0) == 0, "TDL_Channel::discretize(): First tap should be at zero delay");
1449  it_assert(discrete_Ts > 0, "TDL_Channel::discretize(): Incorrect sampling time");
1450  it_assert((a_prof.size() == N_taps) && (delay_profile.size() == N_taps)
1451  && (los_power.size() == N_taps) && (tap_doppler_spectrum.size() == N_taps),
1452  "TDL_Channel:: discretize(): Channel profile lenghts must be equal to the number of taps");
1453 
1454  vec p_prof = sqr(a_prof); // Power profile
1455  ivec delay_prof(N_taps);
1456  vec power(N_taps);
1457  double spower;
1458  vec scattered(N_taps), direct(N_taps);
1459  vec los_doppler(N_taps);
1460  Array <DOPPLER_SPECTRUM> tap_spectrum(N_taps);
1461 
1462  delay_prof(0) = round_i(delay_profile(0) / discrete_Ts); // should be at zero delay anyway
1463  power(0) = p_prof(0);
1464  spower = p_prof(0) / (1 + los_power(0));
1465  scattered(0) = spower;
1466  direct(0) = los_power(0) * spower;
1467  los_doppler(0) = los_dopp(0);
1468  tap_spectrum(0) = tap_doppler_spectrum(0);
1469 
1470  // taps within ((j-0.5)Ts,(j+0.5)Ts] are included in the j-th tap
1471  int j = 0, j_delay = 0;
1472  for (int i = 1; i < N_taps; i++) {
1473  if (delay_profile(i) > (j_delay + 0.5)*discrete_Ts) {
1474  // first skip empty taps
1475  while (delay_profile(i) > (j_delay + 0.5)*discrete_Ts) { j_delay++; }
1476  // create a new tap at (j+1)Ts
1477  j++;
1478  delay_prof(j) = j_delay;
1479  power(j) = p_prof(i);
1480  spower = p_prof(i) / (1 + los_power(i));
1481  scattered(j) = spower;
1482  direct(j) = los_power(i) * spower;
1483  los_doppler(j) = los_dopp(i);
1484  tap_spectrum(j) = tap_doppler_spectrum(i);
1485  }
1486  else {
1487  // add to the previously created tap
1488  power(j) += p_prof(i);
1489  spower = p_prof(i) / (1 + los_power(i));
1490  scattered(j) += spower;
1491  direct(j) += los_power(i) * spower;
1492  it_assert(tap_spectrum(j) == tap_doppler_spectrum(i),
1493  "TDL_Channel::discretize(): Sampling frequency too low. Can not discretize the channel with different Doppler spectra on merged taps.");
1494  it_warning("TDL_Channel::discretize(): Sampling frequency too low. Merging original tap " << i << " with new tap " << j << ".");
1495  if (los_doppler(j) != los_dopp(i)) {
1496  los_doppler(j) = 0.7;
1497  it_warning("TDL_Channel::discretize(): LOS Doppler value reset to 0.7 for tap " << j << " due to the merging process.");
1498  }
1499  }
1500  }
1501 
1502  int no_taps = j + 1; // number of taps found
1503  if (no_taps < N_taps) {
1504  delay_prof.set_size(no_taps, true);
1505  power.set_size(no_taps, true);
1506  direct.set_size(no_taps, true);
1507  scattered.set_size(no_taps, true);
1508  los_doppler.set_size(no_taps, true);
1509  tap_spectrum.set_size(no_taps, true);
1510 
1511  // write over the existing channel profile with its new version
1512  N_taps = no_taps;
1513  a_prof = sqrt(power);
1514  los_power = elem_div(direct, scattered);
1515  los_dopp = los_doppler;
1516  tap_doppler_spectrum = tap_spectrum;
1517  }
1518  // new discretized path's delays
1519  d_prof = delay_prof; // in samples
1520 }
1521 
1522 
1523 // --------------------------------------------------------------------------
1524 // Binary Symetric Channel class methods
1525 // --------------------------------------------------------------------------
1526 
1527 bvec BSC::operator()(const bvec &input)
1528 {
1529  int i, length = input.length();
1530  bvec output(length);
1531 
1532  for (i = 0; i < length; i++) {
1533  if (u() <= p) {
1534  output(i) = input(i) + bin(1);
1535  }
1536  else {
1537  output(i) = input(i);
1538  }
1539  }
1540  return output;
1541 }
1542 
1543 
1544 // --------------------------------------------------------------------------
1545 // AWGN_Channel class methods
1546 // --------------------------------------------------------------------------
1547 
1548 cvec AWGN_Channel::operator()(const cvec &input)
1549 {
1550  int n = input.size();
1551  cvec noise(n);
1552  rng_cn.sample_vector(n, noise);
1553  noise *= sigma;
1554  noise += input;
1555  return noise;
1556 }
1557 
1558 vec AWGN_Channel::operator()(const vec &input)
1559 {
1560  int n = input.size();
1561  vec noise(n);
1562  rng_n.sample_vector(n, noise);
1563  noise *= sigma;
1564  noise += input;
1565  return noise;
1566 }
1567 
1568 
1569 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2