text:
/* vector/gsl_vector_uint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_VECTOR_UINT_H__
#define __GSL_VECTOR_UINT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_uint.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned int *data;
  gsl_block_uint *block;
  int owner;
}
gsl_vector_uint;

typedef struct
{
  gsl_vector_uint vector;
} _gsl_vector_uint_view;

typedef _gsl_vector_uint_view gsl_vector_uint_view;

typedef struct
{
  gsl_vector_uint vector;
} _gsl_vector_uint_const_view;

typedef const _gsl_vector_uint_const_view gsl_vector_uint_const_view;

/* Allocation */

GSL_EXPORT gsl_vector_uint *gsl_vector_uint_alloc (const size_t n);
GSL_EXPORT gsl_vector_uint *gsl_vector_uint_calloc (const size_t n);

GSL_EXPORT gsl_vector_uint *gsl_vector_uint_alloc_from_block (gsl_block_uint * b,
                                                                const size_t offset,
                                                                const size_t n,
                                                                const size_t stride);

GSL_EXPORT gsl_vector_uint *gsl_vector_uint_alloc_from_vector (gsl_vector_uint * v,
                                                                 const size_t offset,
                                                                 const size_t n,
                                                                 const size_t stride);

GSL_EXPORT void gsl_vector_uint_free (gsl_vector_uint * v);

/* Views */

GSL_EXPORT
_gsl_vector_uint_view
gsl_vector_uint_view_array (unsigned int *v, size_t n);

GSL_EXPORT
_gsl_vector_uint_view
gsl_vector_uint_view_array_with_stride (unsigned int *base,
                                         size_t stride,
                                         size_t n);

GSL_EXPORT
_gsl_vector_uint_const_view
gsl_vector_uint_const_view_array (const unsigned int *v, size_t n);

GSL_EXPORT
_gsl_vector_uint_const_view
gsl_vector_uint_const_view_array_with_stride (const unsigned int *base,
                                               size_t stride,
                                               size_t n);

GSL_EXPORT
_gsl_vector_uint_view
gsl_vector_uint_subvector (gsl_vector_uint *v,
                            size_t i,
                            size_t n);

GSL_EXPORT
_gsl_vector_uint_view
gsl_vector_uint_subvector_with_stride (gsl_vector_uint *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_EXPORT
_gsl_vector_uint_const_view
gsl_vector_uint_const_subvector (const gsl_vector_uint *v,
                                  size_t i,
                                  size_t n);

GSL_EXPORT
_gsl_vector_uint_const_view
gsl_vector_uint_const_subvector_with_stride (const gsl_vector_uint *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_EXPORT unsigned int gsl_vector_uint_get (const gsl_vector_uint * v, const size_t i);
GSL_EXPORT void gsl_vector_uint_set (gsl_vector_uint * v, const size_t i, unsigned int x);

GSL_EXPORT unsigned int *gsl_vector_uint_ptr (gsl_vector_uint * v, const size_t i);
GSL_EXPORT const unsigned int *gsl_vector_uint_const_ptr (const gsl_vector_uint * v, const size_t i);

GSL_EXPORT void gsl_vector_uint_set_zero (gsl_vector_uint * v);
GSL_EXPORT void gsl_vector_uint_set_all (gsl_vector_uint * v, unsigned int x);
GSL_EXPORT int gsl_vector_uint_set_basis (gsl_vector_uint * v, size_t i);

GSL_EXPORT int gsl_vector_uint_fread (FILE * stream, gsl_vector_uint * v);
GSL_EXPORT int gsl_vector_uint_fwrite (FILE * stream, const gsl_vector_uint * v);
GSL_EXPORT int gsl_vector_uint_fscanf (FILE * stream, gsl_vector_uint * v);
GSL_EXPORT int gsl_vector_uint_fprintf (FILE * stream, const gsl_vector_uint * v,
                                         const char *format);

GSL_EXPORT int gsl_vector_uint_memcpy (gsl_vector_uint * dest, const gsl_vector_uint * src);

GSL_EXPORT int gsl_vector_uint_reverse (gsl_vector_uint * v);

GSL_EXPORT int gsl_vector_uint_swap (gsl_vector_uint * v, gsl_vector_uint * w);
GSL_EXPORT int gsl_vector_uint_swap_elements (gsl_vector_uint * v, const size_t i, const size_t j);

GSL_EXPORT unsigned int gsl_vector_uint_max (const gsl_vector_uint * v);
GSL_EXPORT unsigned int gsl_vector_uint_min (const gsl_vector_uint * v);
GSL_EXPORT void gsl_vector_uint_minmax (const gsl_vector_uint * v, unsigned int * min_out, unsigned int * max_out);

GSL_EXPORT size_t gsl_vector_uint_max_index (const gsl_vector_uint * v);
GSL_EXPORT size_t gsl_vector_uint_min_index (const gsl_vector_uint * v);
GSL_EXPORT void gsl_vector_uint_minmax_index (const gsl_vector_uint * v, size_t * imin, size_t * imax);

GSL_EXPORT int gsl_vector_uint_add (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_EXPORT int gsl_vector_uint_sub (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_EXPORT int gsl_vector_uint_mul (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_EXPORT int gsl_vector_uint_div (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_EXPORT int gsl_vector_uint_scale (gsl_vector_uint * a, const double x);
GSL_EXPORT int gsl_vector_uint_add_constant (gsl_vector_uint * a, const double x);

GSL_EXPORT int gsl_vector_uint_isnull (const gsl_vector_uint * v);

#ifdef HAVE_INLINE

extern inline
unsigned int
gsl_vector_uint_get (const gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_uint_set (gsl_vector_uint * v, const size_t i, unsigned int x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
unsigned int *
gsl_vector_uint_ptr (gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned int *) (v->data + i * v->stride);
}

extern inline
const unsigned int *
gsl_vector_uint_const_ptr (const gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned int *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_UINT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#include "asf.h"
#include "ceos.h"

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

struct pp_erfin_params {
    int npixels,nlines;
    double satellite_height; /* earth center to satellite height (PP "rsc") */
    double slant_first; /* slant range to first pixel */
    double slant_last; /* slant range to last pixel */
    double nominal_pixsize_range;
};

static double getObjective(double R, void *params)
{
    struct pp_erfin_params *p =
        (struct pp_erfin_params *)params;

    double H=p->satellite_height;
    double sf=p->slant_first;
    double HHRR=H*H + R*R;
    /* Convert slant_first to ground range */
    double g=R*acos( (HHRR - sf*sf) / (2*H*R));
    /* Add width of CEOS image in ground range */
    g += (p->npixels-1) * p->nominal_pixsize_range;
    /* Compute slant range to last pixel */
    double slant_last_R=sqrt(HHRR - 2*H*R*cos(g/R));
    /* Compare with slant_last from CEOS */
    return slant_last_R-p->slant_last;
}

void pp_get_corrected_vals(char *sarName, double *corrected_earth_radius,
                           double *corrected_azimuth_time_per_pixel)
{
    int status;
    int iter = 0, max_iter = 100;
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *s;
    gsl_function F;
    gsl_error_handler_t *prev;
    struct pp_erfin_params params;

    double nominal_pixsize_azimuth;
    double nadir_radius; /* earth radius at nadir */
    double pp_earth_radius; /* f'd up PP Earth radius from start of swath */
    double seconds_per_azimuth_line; /* PP's real azimuth resolution */
    double lo, hi;  /* starting points for the bisection algorithm */
    struct VFDRECV facdr;

    get_asf_facdr(sarName, &facdr);

    /* Find the PP's earth radius with an iterative search */
    F.function = &getObjective;
    F.params = &params;

    params.npixels=facdr.npixels;
    params.nlines=facdr.nlines;
    params.nominal_pixsize_range=facdr.rapixspc;
    nominal_pixsize_azimuth=facdr.azpixspc;
    nadir_radius=1.0e3*facdr.eradnadr;
    params.satellite_height=1.0e3*(facdr.scalt+facdr.eradnadr);
    params.slant_first=1.0e3*facdr.sltrngfp;
    params.slant_last=1.0e3*facdr.sltrnglp;

    prev = gsl_set_error_handler_off();

    lo = nadir_radius - 2000;
    hi = nadir_radius + 2000;

    T = gsl_root_fsolver_brent;
    s = gsl_root_fsolver_alloc (T);
    gsl_root_fsolver_set (s, &F, lo, hi);

    do {
        ++iter;
        status = gsl_root_fsolver_iterate(s);
        pp_earth_radius = gsl_root_fsolver_root(s);
        status = gsl_root_test_residual(
            getObjective(pp_earth_radius, (void*)&params), 1.0e-4);
    } while (status == GSL_CONTINUE && iter < max_iter);

    if (status == GSL_SUCCESS) {
        //printf("Converged after %d iterations.\n", iter);
        //printf("PP Earth Radius: %.3f m\n",pp_earth_radius);
        //printf("   (for comparison) Nadir Earth Radius: %.3f m\n",nadir_radius);
        *corrected_earth_radius = pp_earth_radius;
    } else {
        asfPrintWarning("Failed to determine PP earth radius!\n"
                        "iter: %d, pp_earth_radius=%.3f, res=%.5f\n"
                        "Proceeding using the nadir radius: %.3f m\n"
                        "Starting points were: lo: %.3f -> %.4f\n"
                        "                      hi: %.3f -> %.4f\n",
                        iter, pp_earth_radius,
                        getObjective(pp_earth_radius, (void*)&params),
                        nadir_radius,
                        lo, getObjective(lo, (void*)&params),
                        hi, getObjective(hi, (void*)&params));
        *corrected_earth_radius = nadir_radius;
    }

    gsl_set_error_handler(prev);

    // Find the PP's per-second azimuth pixel spacing
    seconds_per_azimuth_line=nominal_pixsize_azimuth/facdr.swathvel;
    //printf("PP seconds per azimuth line: %.9f s/line\n",seconds_per_azimuth_line);
    //printf("   (for comparison) PP interpolated lines per second: %.3f lines/s\n",1.0/seconds_per_azimuth_line);
    //printf("   (for comparison) FACDR swath velocity: %.3f m/s\n",facdr.swathvel);

    //double R=pp_earth_radius;
    //double H=params.satellite_height;
    //double HHRR=H*H + R*R;
    //double slant=0.5*(params.slant_first+params.slant_last);
    //double rg_center=R*acos( (HHRR - slant*slant) / (2*H*R));
    //double vs=sqrt(facdr.scxvel*facdr.scxvel + facdr.scyvel*facdr.scyvel + facdr.sczvel*facdr.sczvel);
    //double vsg = vs * pp_earth_radius/params.satellite_height*cos(rg_center/pp_earth_radius);

    //printf("   (for comparison) PP-style recalc velocity: %.3f m/s\n",vsg);

    *corrected_azimuth_time_per_pixel = seconds_per_azimuth_line;

    // Free the solver
    gsl_root_fsolver_free(s);
}

ext: c

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_dot (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int N = 1;
   float alpha = 0.0f;
   float X[] = { 0.733f };
   float Y[] = { 0.825f };
   int incX = 1;
   int incY = -1;
   float expected = 0.604725f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 1)");
  };

  {
   int N = 1;
   float alpha = 0.1f;
   float X[] = { 0.733f };
   float Y[] = { 0.825f };
   int incX = 1;
   int incY = -1;
   float expected = 0.704725f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 2)");
  };

  {
   int N = 1;
   float alpha = 1.0f;
   float X[] = { 0.733f };
   float Y[] = { 0.825f };
   int incX = 1;
   int incY = -1;
   float expected = 1.604725f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 3)");
  };

  {
   int N = 1;
   float alpha = 0.0f;
   float X[] = { -0.812f };
   float Y[] = { -0.667f };
   int incX = -1;
   int incY = 1;
   float expected = 0.541604f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 4)");
  };

  {
   int N = 1;
   float alpha = 0.1f;
   float X[] = { -0.812f };
   float Y[] = { -0.667f };
   int incX = -1;
   int incY = 1;
   float expected = 0.641604f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 5)");
  };

  {
   int N = 1;
   float alpha = 1.0f;
   float X[] = { -0.812f };
   float Y[] = { -0.667f };
   int incX = -1;
   int incY = 1;
   float expected = 1.541604f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 6)");
  };

  {
   int N = 1;
   float alpha = 0.0f;
   float X[] = { 0.481f };
   float Y[] = { 0.523f };
   int incX = -1;
   int incY = -1;
   float expected = 0.251563f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 7)");
  };

  {
   int N = 1;
   float alpha = 0.1f;
   float X[] = { 0.481f };
   float Y[] = { 0.523f };
   int incX = -1;
   int incY = -1;
   float expected = 0.351563f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 8)");
  };

  {
   int N = 1;
   float alpha = 1.0f;
   float X[] = { 0.481f };
   float Y[] = { 0.523f };
   int incX = -1;
   int incY = -1;
   float expected = 1.251563f;
   float f;
   f = cblas_sdsdot (N, alpha, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdsdot(case 9)");
  };

  {
   int N = 1;
   float X[] = { 0.785f };
   float Y[] = { -0.7f };
   int incX = 1;
   int incY = -1;
   float expected = -0.5495f;
   float f;
   f = cblas_sdot(N, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdot(case 10)");
  };

  {
   int N = 1;
   double X[] = { 0.79 };
   double Y[] = { -0.679 };
   int incX = 1;
   int incY = -1;
   double expected = -0.53641;
   double f;
   f = cblas_ddot(N, X, incX, Y, incY);
   gsl_test_rel(f, expected, dbleps, "ddot(case 11)");
  };

  {
   int N = 1;
   float X[] = { 0.474f, -0.27f };
   float Y[] = { -0.144f, -0.392f };
   int incX = 1;
   int incY = -1;
   float expected[2] = {-0.174096f, -0.146928f};
   float f[2];
   cblas_cdotu_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], flteps, "cdotu(case 12) real");
   gsl_test_rel(f[1], expected[1], flteps, "cdotu(case 12) imag");
  };

  {
   int N = 1;
   float X[] = { 0.474f, -0.27f };
   float Y[] = { -0.144f, -0.392f };
   int incX = 1;
   int incY = -1;
   float expected[2] = {0.037584f, -0.224688f};
   float f[2];
   cblas_cdotc_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], flteps, "cdotc(case 13) real");
   gsl_test_rel(f[1], expected[1], flteps, "cdotc(case 13) imag");
  };

  {
   int N = 1;
   double X[] = { -0.87, -0.631 };
   double Y[] = { -0.7, -0.224 };
   int incX = 1;
   int incY = -1;
   double expected[2] = {0.467656, 0.63658};
   double f[2];
   cblas_zdotu_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], dbleps, "zdotu(case 14) real");
   gsl_test_rel(f[1], expected[1], dbleps, "zdotu(case 14) imag");
  };

  {
   int N = 1;
   double X[] = { -0.87, -0.631 };
   double Y[] = { -0.7, -0.224 };
   int incX = 1;
   int incY = -1;
   double expected[2] = {0.750344, -0.24682};
   double f[2];
   cblas_zdotc_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], dbleps, "zdotc(case 15) real");
   gsl_test_rel(f[1], expected[1], dbleps, "zdotc(case 15) imag");
  };

  {
   int N = 1;
   float X[] = { -0.457f };
   float Y[] = { 0.839f };
   int incX = -1;
   int incY = 1;
   float expected = -0.383423f;
   float f;
   f = cblas_sdot(N, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdot(case 16)");
  };

  {
   int N = 1;
   double X[] = { 0.949 };
   double Y[] = { -0.873 };
   int incX = -1;
   int incY = 1;
   double expected = -0.828477;
   double f;
   f = cblas_ddot(N, X, incX, Y, incY);
   gsl_test_rel(f, expected, dbleps, "ddot(case 17)");
  };

  {
   int N = 1;
   float X[] = { 0.852f, -0.045f };
   float Y[] = { 0.626f, -0.164f };
   int incX = -1;
   int incY = 1;
   float expected[2] = {0.525972f, -0.167898f};
   float f[2];
   cblas_cdotu_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], flteps, "cdotu(case 18) real");
   gsl_test_rel(f[1], expected[1], flteps, "cdotu(case 18) imag");
  };

  {
   int N = 1;
   float X[] = { 0.852f, -0.045f };
   float Y[] = { 0.626f, -0.164f };
   int incX = -1;
   int incY = 1;
   float expected[2] = {0.540732f, -0.111558f};
   float f[2];
   cblas_cdotc_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], flteps, "cdotc(case 19) real");
   gsl_test_rel(f[1], expected[1], flteps, "cdotc(case 19) imag");
  };

  {
   int N = 1;
   double X[] = { -0.786, -0.341 };
   double Y[] = { -0.271, -0.896 };
   int incX = -1;
   int incY = 1;
   double expected[2] = {-0.09253, 0.796667};
   double f[2];
   cblas_zdotu_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], dbleps, "zdotu(case 20) real");
   gsl_test_rel(f[1], expected[1], dbleps, "zdotu(case 20) imag");
  };

  {
   int N = 1;
   double X[] = { -0.786, -0.341 };
   double Y[] = { -0.271, -0.896 };
   int incX = -1;
   int incY = 1;
   double expected[2] = {0.518542, 0.611845};
   double f[2];
   cblas_zdotc_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], dbleps, "zdotc(case 21) real");
   gsl_test_rel(f[1], expected[1], dbleps, "zdotc(case 21) imag");
  };

  {
   int N = 1;
   float X[] = { -0.088f };
   float Y[] = { -0.165f };
   int incX = -1;
   int incY = -1;
   float expected = 0.01452f;
   float f;
   f = cblas_sdot(N, X, incX, Y, incY);
   gsl_test_rel(f, expected, flteps, "sdot(case 22)");
  };

  {
   int N = 1;
   double X[] = { -0.434 };
   double Y[] = { -0.402 };
   int incX = -1;
   int incY = -1;
   double expected = 0.174468;
   double f;
   f = cblas_ddot(N, X, incX, Y, incY);
   gsl_test_rel(f, expected, dbleps, "ddot(case 23)");
  };

  {
   int N = 1;
   float X[] = { -0.347f, 0.899f };
   float Y[] = { -0.113f, -0.858f };
   int incX = -1;
   int incY = -1;
   float expected[2] = {0.810553f, 0.196139f};
   float f[2];
   cblas_cdotu_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], flteps, "cdotu(case 24) real");
   gsl_test_rel(f[1], expected[1], flteps, "cdotu(case 24) imag");
  };

  {
   int N = 1;
   float X[] = { -0.347f, 0.899f };
   float Y[] = { -0.113f, -0.858f };
   int incX = -1;
   int incY = -1;
   float expected[2] = {-0.732131f, 0.399313f};
   float f[2];
   cblas_cdotc_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], flteps, "cdotc(case 25) real");
   gsl_test_rel(f[1], expected[1], flteps, "cdotc(case 25) imag");
  };

  {
   int N = 1;
   double X[] = { -0.897, -0.204 };
   double Y[] = { -0.759, 0.557 };
   int incX = -1;
   int incY = -1;
   double expected[2] = {0.794451, -0.344793};
   double f[2];
   cblas_zdotu_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], dbleps, "zdotu(case 26) real");
   gsl_test_rel(f[1], expected[1], dbleps, "zdotu(case 26) imag");
  };

  {
   int N = 1;
   double X[] = { -0.897, -0.204 };
   double Y[] = { -0.759, 0.557 };
   int incX = -1;
   int incY = -1;
   double expected[2] = {0.567195, -0.654465};
   double f[2];
   cblas_zdotc_sub(N, X, incX, Y, incY, &f);
   gsl_test_rel(f[0], expected[0], dbleps, "zdotc(case 27) real");
   gsl_test_rel(f[1], expected[1], dbleps, "zdotc(case 27) imag");
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#ifndef BIO_MATH_H_
#define BIO_MATH_H_

#include "bio/defs.h"
#include "bio/gsl.h"

#include <boost/io/ios_state.hpp>

#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_exp.h>

#include <numeric>

BIO_NS_START

/** Calculates the log of n choose the observed counts.  */
template < typename ObsIt >
double
calc_ln_n_choose(
	ObsIt obs_begin,
	ObsIt obs_end)
{
	double result = 0.0;

	unsigned n = 0;
	while (obs_end != obs_begin)
	{
		result -= gsl_sf_lnfact(*obs_begin);
		n += *obs_begin;

		++obs_begin;
	}

	result += gsl_sf_lnfact(n);

	return result;
}

/** Calculates a factor involved in the likelihood of a multinomial distribution with a dirichlet prior. */
template < typename ObsIt, typename AlphaIt >
double
calc_ln_gamma_factor(
	ObsIt obs_begin,
	ObsIt obs_end,
	AlphaIt alpha_begin)
{
	double ln_pi_gamma_alpha_obs = 0.0; //work in logarithms
	double ln_pi_gamma_alpha = 0.0; //work in logarithms
	double sum_alpha = 0.0;
	double sum_obs = 0.0;
	while (obs_end != obs_begin)
	{
		sum_alpha += *alpha_begin;
		sum_obs += *obs_begin;
		ln_pi_gamma_alpha += gsl_sf_lngamma(*alpha_begin);
		ln_pi_gamma_alpha_obs += gsl_sf_lngamma(*alpha_begin + *obs_begin);

		++obs_begin;
		++alpha_begin;
	}

	double ln_result =
		gsl_sf_lngamma(sum_alpha)
		- ln_pi_gamma_alpha
		+ ln_pi_gamma_alpha_obs
		- gsl_sf_lngamma(sum_alpha + sum_obs);

	return ln_result;
}

/** Calculates the log of the likelihood of seeing the observed counts under a dirichlet prior. */
template < typename ObsIt, typename AlphaIt >
double
calc_multinomial_ln_likelihood_dirichlet_prior(
	ObsIt obs_begin,
	ObsIt obs_end,
	AlphaIt alpha_begin)
{
	return calc_ln_gamma_factor(obs_begin, obs_end, alpha_begin) + calc_ln_n_choose(obs_begin, obs_end);
}

/** Calculates the log of the likelihood of seeing the observed counts under a uniform multinomial distribution. */
template < typename ObsIt >
double
calc_multinomial_ln_likelihood_uniform_dist(
	ObsIt obs_begin,
	ObsIt obs_end)
{
	const unsigned n = std::accumulate(obs_begin, obs_end, unsigned(0));
	if (0 == n)
	{
		return 0.0;
	}

	const double ln_n_choose = calc_ln_n_choose(obs_begin, obs_end);
	return ln_n_choose - n * gsl_sf_log(double(obs_end - obs_begin));
}

/** Calculates the log of the likelihood of seeing the observed counts under a multinomial distribution with given p's. */
template < typename ObsIt, typename PIt >
double
calc_multinomial_ln_likelihood(
	ObsIt obs_begin,
	ObsIt obs_end,
	PIt p_begin)
{
	const unsigned n = std::accumulate(obs_begin, obs_end, unsigned(0));
	if (0 == n)
	{
		return 0.0;
	}

	double result = calc_ln_n_choose(obs_begin, obs_end);
	while (obs_end != obs_begin)
	{
		result += *obs_begin * gsl_sf_log(*p_begin);

		++obs_begin;
		++p_begin;
	}

	return result;
}

/** Calculates the ratio of evidence in favour of a uniform distribution over one with a dirichlet prior
defined by the alphas. */
template < typename ObsIt, typename AlphaIt >
double
calc_ln_multinomial_uniform_evidence(
	ObsIt obs_begin,
	ObsIt obs_end,
	AlphaIt alpha_begin)
{
	const double ln_likelihood_uniform_dist = calc_multinomial_ln_likelihood_uniform_dist(obs_begin, obs_end);
	const double ln_likelihood_dirichlet_prior = calc_multinomial_ln_likelihood_dirichlet_prior(obs_begin, obs_end, alpha_begin);
	return ln_likelihood_uniform_dist - ln_likelihood_dirichlet_prior;
}

/** Calculates the ratio of evidence in favour of a particular multinomial distribution (defined by the p's)
over one with a dirichlet prior defined by the alpha's. */
template < typename ObsIt, typename AlphaIt, typename PIt >
double
calc_ln_multinomial_evidence(
	ObsIt obs_begin,
	ObsIt obs_end,
	AlphaIt alpha_begin,
	PIt p_begin)
{
	return
		calc_multinomial_ln_likelihood(obs_begin, obs_end, p_begin)
		- calc_multinomial_ln_likelihood_dirichlet_prior(obs_begin, obs_end, alpha_begin);
}

BIO_NS_END

#endif //BIO_MATH_H_

ext: h

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
#ifndef NETWORK_UDP_ASYNCRECEIVER_H
#define NETWORK_UDP_ASYNCRECEIVER_H

/* Base class for implementing an asynchronous UDP receiver. A derived class must implement
   handle_receive which gets passed a view of the received data. The function start() will
   block and should be called in its own thread. */

#include <cstdint>
#include <array>
#include <string>

#include <gsl/gsl>
#include <asio.hpp>

namespace network::udp {

class Async_receiver
{
public:
  Async_receiver() : socket_{io_service_} {}
  virtual ~Async_receiver();
  Async_receiver(const Async_receiver&) = delete;
  Async_receiver& operator=(const Async_receiver&) = delete;

  void start(const std::string& ip, std::uint16_t port);
  void stop();
  bool is_running() const { return socket_.is_open(); }

protected:
  virtual void handle_receive(gsl::span<std::uint8_t> buffer) = 0;

private:
  void start_async_receive();

  asio::io_service io_service_;
  asio::ip::udp::socket socket_;
  std::array<std::uint8_t, 128> buffer_;
};

}  // namespace network::udp

#endif  // NETWORK_UDP_ASYNCRECEIVER_H

max_stars_repo_name: mwkpe/tincan

max_issues_repo_name: jwkpeter/tincan

max_forks_repo_name: jwkpeter/tincan


text:
/* ode-initval/rk4.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Runge-Kutta 4th order, Classical */

/* Author:  G. Jungman
 */

/* Reference: Abramowitz & Stegun, section 25.5. equation 25.5.10

   Error estimation by step doubling, see eg. Ascher, U.M., Petzold,
   L.R., Computer methods for ordinary differential and
   differential-algebraic equations, SIAM, Philadelphia, 1998.
*/

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k;
  double *k1;
  double *y0;
  double *ytmp;
  double *y_onestep;
}
rk4_state_t;

static void *
rk4_alloc (size_t dim)
{
  rk4_state_t *state = (rk4_state_t *) malloc (sizeof (rk4_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk4_state", GSL_ENOMEM);
    }

  state->k = (double *) malloc (dim * sizeof (double));

  if (state->k == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
    }

  state->k1 = (double *) malloc (dim * sizeof (double));

  if (state->k1 == 0)
    {
      free (state->k);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->y0);
      free (state->k);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      free (state->ytmp);
      free (state->y0);
      free (state->k);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  return state;
}

static int
rk4_step (double *y, const rk4_state_t *state,
	  const double h, const double t, const size_t dim,
	  const gsl_odeiv_system *sys)
{
  /* Makes a Runge-Kutta 4th order advance with step size h. */

  /* initial values of variables y. */
  const double *y0 = state->y0;

  /* work space */
  double *ytmp = state->ytmp;

  /* Runge-Kutta coefficients. Contains values of coefficient k1
     in the beginning
  */
  double *k = state->k;

  size_t i;

  /* k1 step */

  for (i = 0; i < dim; i++)
    {
      y[i] += h / 6.0 * k[i];
      ytmp[i] = y0[i] + 0.5 * h * k[i];
    }

  /* k2 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, k);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    {
      y[i] += h / 3.0 * k[i];
      ytmp[i] = y0[i] + 0.5 * h * k[i];
    }

  /* k3 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, k);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    {
      y[i] += h / 3.0 * k[i];
      ytmp[i] = y0[i] + h * k[i];
    }

  /* k4 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    {
      y[i] += h / 6.0 * k[i];
    }

  return GSL_SUCCESS;
}

static int
rk4_apply (void *vstate,
           size_t dim,
           double t,
           double h,
           double y[],
           double yerr[],
           const double dydt_in[],
           double dydt_out[],
           const gsl_odeiv_system * sys)
{
  rk4_state_t *state = (rk4_state_t *) vstate;

  size_t i;

  double *const k = state->k;
  double *const k1 = state->k1;
  double *const y0 = state->y0;
  double *const y_onestep = state->y_onestep;

  DBL_MEMCPY (y0, y, dim);

  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y0, k);

      if (s != GSL_SUCCESS)
	{
	  return s;
	}
    }

  /* Error estimation is done by step doubling procedure */

  /* Save first point derivatives*/

  DBL_MEMCPY (k1, k, dim);

  /* First traverse h with one step (save to y_onestep) */

  DBL_MEMCPY (y_onestep, y, dim);

  {
    int s = rk4_step (y_onestep, state, h, t, dim, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Then with two steps with half step length (save to y) */

  DBL_MEMCPY (k, k1, dim);

  {
    int s = rk4_step (y, state, h/2.0, t, dim, sys);

    if (s != GSL_SUCCESS)
      {
	/* Restore original values */
	DBL_MEMCPY (y, y0, dim);
	return s;
    }
  }

  /* Update before second step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h/2.0, y, k);

    if (s != GSL_SUCCESS)
      {
	/* Restore original values */
	DBL_MEMCPY (y, y0, dim);
	return s;
      }
  }

  /* Save original y0 to k1 for possible failures */
  DBL_MEMCPY (k1, y0, dim);

  /* Update y0 for second step */
  DBL_MEMCPY (y0, y, dim);

  {
    int s = rk4_step (y, state, h/2.0, t + h/2.0, dim, sys);

    if (s != GSL_SUCCESS)
      {
	/* Restore original values */
	DBL_MEMCPY (y, k1, dim);
	return s;
      }
  }

  /* Derivatives at output */

  if (dydt_out != NULL) {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

    if (s != GSL_SUCCESS)
      {
	/* Restore original values */
	DBL_MEMCPY (y, k1, dim);
	return s;
      }
  }

  /* Error estimation

     yerr = C * 0.5 * | y(onestep) - y(twosteps) | / (2^order - 1)

     constant C is approximately 8.0 to ensure 90% of samples lie within
     the error (assuming a gaussian distribution with prior p(sigma)=1/sigma.)

  */

  for (i = 0; i < dim; i++)
    {
      yerr[i] = 4.0 * (y[i] - y_onestep[i]) / 15.0;
    }

  return GSL_SUCCESS;
}

static int
rk4_reset (void *vstate, size_t dim)
{
  rk4_state_t *state = (rk4_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k, dim);
  DBL_ZERO_MEMSET (state->k1, dim);
  DBL_ZERO_MEMSET (state->y0, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);
  DBL_ZERO_MEMSET (state->y_onestep, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk4_order (void *vstate)
{
  rk4_state_t *state = (rk4_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 4;
}

static void
rk4_free (void *vstate)
{
  rk4_state_t *state = (rk4_state_t *) vstate;
  free (state->k);
  free (state->k1);
  free (state->y0);
  free (state->ytmp);
  free (state->y_onestep);
  free (state);
}

static const gsl_odeiv_step_type rk4_type = { "rk4",    /* name */
  1,                            /* can use dydt_in */
  1,                            /* gives exact dydt_out */
  &rk4_alloc,
  &rk4_apply,
  &rk4_reset,
  &rk4_order,
  &rk4_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rk4 = &rk4_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
// Copyright 2016-2018 Lauri Juvela and Manu Airaksinen
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef DEFINITIONS_H_
#define DEFINITIONS_H_

#include <gslwrap/vector_int.h>
#include <gslwrap/vector_double.h>
#include <gslwrap/matrix_double.h>

/* Enums */
enum DataType {ASCII, DOUBLE, FLOAT};
enum SignalPolarity {POLARITY_DEFAULT, POLARITY_INVERT, POLARITY_DETECT};
enum LpWeightingFunction {NONE, AME, STE};
enum WindowingFunctionType {HANN, HAMMING, BLACKMAN, COSINE, HANNING, RECT, NUTTALL};
enum ExcitationMethod {SINGLE_PULSE_EXCITATION, DNN_GENERATED_EXCITATION,
   PULSES_AS_FEATURES_EXCITATION, EXTERNAL_EXCITATION, IMPULSE_EXCITATION};

/* Structures */
struct Param
{
	Param();
	~Param();

  public:
	int fs;
	int frame_length;
	int frame_length_long;
	int frame_length_unvoiced;
	int frame_shift;
	int number_of_frames;
	int signal_length;
	int lpc_order_vt;
	int lpc_order_glot;
	int hnr_order;
	bool use_external_f0;
	std::string external_f0_filename;
	bool use_external_gci;
	std::string external_gci_filename;
	bool use_external_lsf_vt; // inverse filtering with external vocal tract filter
	std::string external_lsf_vt_filename;
	bool use_external_excitation;
	std::string external_excitation_filename;

	std::string dnn_path_basename;
	std::string data_directory;
	std::string file_basename;
	std::string file_path;
	bool save_to_datadir_root;
	DataType data_type;
	bool qmf_subband_analysis;
	SignalPolarity signal_polarity;
	double gif_pre_emphasis_coefficient;
	double unvoiced_pre_emphasis_coefficient;
	bool use_iterative_gif;
	LpWeightingFunction lp_weighting_function;
	int lpc_order_glot_iaif;
	double warping_lambda_vt;
	double ame_duration_quotient;
	double ame_position_quotient;
	WindowingFunctionType default_windowing_function;
	WindowingFunctionType psola_windowing_function;
	WindowingFunctionType paf_analysis_window;

	double max_pulse_len_diff;
	int paf_pulse_length;
	bool use_pulse_interpolation;
	bool use_highpass_filtering;
	bool use_waveforms_directly;
	bool use_paf_unvoiced_synthesis;
	bool use_velvet_unvoiced_paf;
	bool extract_f0;
	bool extract_gain;
	bool extract_lsf_vt;
	bool extract_lsf_glot;
	bool extract_hnr;
	bool extract_infofile;
	bool extract_glottal_excitation;
	bool extract_original_signal;
	bool extract_gci_signal;
	bool extract_pulses_as_features;
	bool use_paf_energy_normalization;
	int lpc_order_vt_qmf1;
	int lpc_order_vt_qmf2;
	double f0_max;
	double f0_min;
	double voicing_threshold;
	double zcr_threshold;
	double relative_f0_threshold;
	double speed_scale;
	double pitch_scale;
	bool use_postfiltering;
	bool use_spectral_matching;
	bool use_wsola;
	bool use_wsola_pitch_shift;
	bool noise_gated_synthesis;
	double noise_reduction_db;
	double noise_gate_limit_db;
	double postfilter_coefficient;
	double postfilter_coefficient_glot;
	bool use_trajectory_smoothing;
	int lsf_vt_smooth_len;
	int lsf_glot_smooth_len;
	int gain_smooth_len;
	int hnr_smooth_len;
	int filter_update_interval_vt;
	int filter_update_interval_specmatch;
	double noise_gain_unvoiced;
	double noise_gain_voiced;
	double noise_low_freq_limit_voiced;
	//double f0_check_range;
	ExcitationMethod excitation_method;
	bool use_pitch_synchronous_analysis;
	bool use_generic_envelope;

	/* directory paths for storing parameters */
	std::string dir_gain;
	std::string dir_lsf;
	std::string dir_lsfg;
	std::string dir_hnr;
	std::string dir_paf;
	std::string dir_f0;
	std::string dir_exc;
	std::string dir_syn;
	std::string dir_sp;

	/* extensions for parameter types */
	std::string extension_gain;
	std::string extension_lsf;
	std::string extension_lsfg;
	std::string extension_hnr;
	std::string extension_paf;
	std::string extension_f0;
	std::string extension_exc;
	std::string extension_src = ".src.wav";
	std::string extension_syn;
	std::string extension_wav;

	std::string wav_filename;
	std::string default_config_filename;
	std::string user_config_filename;
};

/* Define analysis data variable struct*/
struct AnalysisData {
	AnalysisData();
	~AnalysisData();
	int AllocateData(const Param &params);
	int SaveData(const Param &params);
public:
	gsl::vector signal;
	gsl::vector fundf;
	gsl::vector frame_energy;
	gsl::vector_int gci_inds;
	gsl::vector source_signal;
   gsl::vector source_signal_iaif;

	gsl::matrix poly_vocal_tract;
	gsl::matrix lsf_vocal_tract;
	gsl::matrix poly_glot;
	gsl::matrix lsf_glot;
	gsl::matrix excitation_pulses;
   gsl::matrix hnr_glot;

	/* QMF analysis specific */
	//gsl::matrix lsf_vt_qmf1;
	//gsl::matrix lsf_vt_qmf2;
	//gsl::vector gain_qmf;
};

/* Define analysis data variable struct*/
struct SynthesisData {
   SynthesisData();
   ~SynthesisData();
public:
   gsl::vector signal;
   gsl::vector fundf;
   gsl::vector frame_energy;
   gsl::vector excitation_signal;

   gsl::matrix poly_vocal_tract;
   gsl::matrix lsf_vocal_tract;
   gsl::matrix poly_glot;
   gsl::matrix lsf_glot;
   gsl::matrix excitation_pulses;
   gsl::matrix hnr_glot;

   gsl::matrix spectrum;

   /* QMF analysis specific */
   //gsl::matrix lsf_vt_qmf1;
   //gsl::matrix lsf_vt_qmf2;
   //gsl::vector gain_qmf;
};

#endif

max_stars_repo_name: entn-at/GlottDNN

max_issues_repo_name: entn-at/GlottDNN

max_forks_repo_name: entn-at/GlottDNN


text:
/*
 * Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The OpenAirInterface Software Alliance licenses this file to You under
 * the OAI Public License, Version 1.1  (the "License"); you may not use this file
 * except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.openairinterface.org/?page_id=698
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *-------------------------------------------------------------------------------
 * For more information about the OpenAirInterface (OAI) Software Alliance:
 *      contact@openairinterface.org
 */

#include <math.h>
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "PHY/TOOLS/tools_defs.h"
#include "sim.h"
#include "scm_corrmat.h"
#include "common/utils/LOG/log.h"
//#define DEBUG_CH

#include "assertions.h"

extern void print_shorts(char *s,__m128i *x);

void fill_channel_desc(channel_desc_t *chan_desc,
                       uint8_t nb_tx,
                       uint8_t nb_rx,
                       uint8_t nb_taps,
                       uint8_t channel_length,
                       double *amps,
                       double *delays,
                       struct complex** R_sqrt,
                       double Td,
                       double sampling_rate,
                       double channel_bandwidth,
                       double ricean_factor,
                       double aoa,
                       double forgetting_factor,
                       double max_Doppler,
                       int32_t channel_offset,
                       double path_loss_dB,
                       uint8_t random_aoa)
{

  uint16_t i,j;
  double delta_tau;

  LOG_I(OCM,"[CHANNEL] Getting new channel descriptor, nb_tx %d, nb_rx %d, nb_taps %d, channel_length %d\n",
      nb_tx,nb_rx,nb_taps,channel_length);

  chan_desc->nb_tx          = nb_tx;
  chan_desc->nb_rx          = nb_rx;
  chan_desc->nb_taps        = nb_taps;
  chan_desc->channel_length = channel_length;
  chan_desc->amps           = amps;
  LOG_D(OCM,"[CHANNEL] Doing delays ...\n");
  if (delays==NULL) {
    chan_desc->delays = (double*) malloc(nb_taps*sizeof(double));
    delta_tau = Td/nb_taps;
    for (i=0; i<nb_taps; i++)
      chan_desc->delays[i] = ((double)i)*delta_tau;
  }
  else
  chan_desc->delays              = delays;
  chan_desc->Td                         = Td;
  chan_desc->sampling_rate              = sampling_rate;
  chan_desc->channel_bandwidth          = channel_bandwidth;
  chan_desc->ricean_factor              = ricean_factor;
  chan_desc->aoa                        = aoa;
  chan_desc->random_aoa                 = random_aoa;
  chan_desc->forgetting_factor          = forgetting_factor;
  chan_desc->channel_offset             = channel_offset;
  chan_desc->path_loss_dB               = path_loss_dB;
  chan_desc->first_run                  = 1;
  chan_desc->ip                         = 0.0;
  chan_desc->max_Doppler                = max_Doppler;
  chan_desc->ch                         = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
  chan_desc->chF                        = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
  chan_desc->a                          = (struct complex**) malloc(nb_taps*sizeof(struct complex*));

  LOG_D(OCM,"[CHANNEL] Filling ch \n");

  for (i = 0; i<nb_tx*nb_rx; i++)
    chan_desc->ch[i] = (struct complex*) malloc(channel_length * sizeof(struct complex));

  for (i = 0; i<nb_tx*nb_rx; i++)
    chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));  // allocate for up to 100 RBs, 12 samples per RB

  LOG_D(OCM,"[CHANNEL] Filling a (nb_taps %d)\n",nb_taps);
  for (i = 0; i<nb_taps; i++) {
    LOG_D(OCM,"tap %d (%p,%zu)\n",i,&chan_desc->a[i],nb_tx*nb_rx * sizeof(struct complex));
    chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
  }

  LOG_D(OCM,"[CHANNEL] Doing R_sqrt ...\n");
  if (R_sqrt == NULL) {
      chan_desc->R_sqrt         = (struct complex**) calloc(nb_taps,sizeof(struct complex*));
    for (i = 0; i<nb_taps; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) calloc(nb_tx*nb_rx*nb_tx*nb_rx,sizeof(struct complex));
      for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
        chan_desc->R_sqrt[i][j].x = 1.0;
        chan_desc->R_sqrt[i][j].y = 0.0;
      }
    }
  }
  else {
    chan_desc->R_sqrt = (struct complex**) calloc(nb_taps,sizeof(struct complex*));
    for (i = 0; i<nb_taps; i++) {
        //chan_desc->R_sqrt[i]    = (struct complex*) calloc(nb_tx*nb_rx*nb_tx*nb_rx,sizeof(struct complex));
        //chan_desc->R_sqrt = (struct complex*)&R_sqrt[i][0];
        /* all chan_desc share the same R_sqrt, coming from caller */
        chan_desc->R_sqrt[i] = R_sqrt[0];
        }
  }

  for (i = 0; i<nb_taps; i++) {
    for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
      LOG_D(OCM,"Rsqrt[%d][%d] %f %f\n",i,j,chan_desc->R_sqrt[i][j].x,chan_desc->R_sqrt[i][j].y);
    }
  }

  LOG_D(OCM,"[CHANNEL] RF %f\n",chan_desc->ricean_factor);
  for (i=0;i<chan_desc->nb_taps;i++)
    LOG_D(OCM,"[CHANNEL] tap %d: amp %f, delay %f\n",i,chan_desc->amps[i],chan_desc->delays[i]);

  chan_desc->nb_paths=10;

  reset_meas(&chan_desc->random_channel);
  reset_meas(&chan_desc->interp_time);
  reset_meas(&chan_desc->interp_freq);
  reset_meas(&chan_desc->convolution);

}

double mbsfn_delays[] = {0,.03,.15,.31,.37,1.09,12.490,12.52,12.64,12.80,12.86,13.58,27.49,27.52,27.64,27.80,27.86,28.58};
double mbsfn_amps_dB[] = {0,-1.5,-1.4,-3.6,-0.6,-7.0,-10,-11.5,-11.4,-13.6,-10.6,-17.0,-20,-21.5,-21.4,-23.6,-20.6,-27};

double scm_c_delays[] = {0, 0.0125, 0.0250, 0.3625, 0.3750, 0.3875, 0.2500, 0.2625, 0.2750, 1.0375, 1.0500, 1.0625, 2.7250, 2.7375, 2.7500, 4.6000, 4.6125, 4.6250};
double scm_c_amps_dB[] = {0.00, -2.22, -3.98, -1.86, -4.08, -5.84, -1.08, -3.30, -5.06, -9.08, -11.30, -13.06, -15.14, -17.36, -19.12, -20.64, -22.85, -24.62};

double epa_delays[] = { 0,.03,.07,.09,.11,.19,.41};
double epa_amps_dB[] = {0.0,-1.0,-2.0,-3.0,-8.0,-17.2,-20.8};

double eva_delays[] = { 0,.03,.15,.31,.37,.71,1.09,1.73,2.51};
double eva_amps_dB[] = {0.0,-1.5,-1.4,-3.6,-0.6,-9.1,-7.0,-12.0,-16.9};

double etu_delays[] = { 0,.05,.12,.2,.23,.5,1.6,2.3,5.0};
double etu_amps_dB[] = {-1.0,-1.0,-1.0,0.0,0.0,0.0,-3.0,-5.0,-7.0};

double default_amps_lin[] = {0.3868472 , 0.3094778 , 0.1547389 , 0.0773694 , 0.0386847 , 0.0193424 , 0.0096712 , 0.0038685};
double default_amp_lin[] = {1};

double ts_shift_delays[] = {0, 1/7.68};
double ts_shift_amps[] = {0, 1};

//correlation matrix for a 2x2 channel with full Tx correlation
struct complex R_sqrt_22_corr_tap[16] = {{0.70711,0}, {0.0, 0.0}, {0.70711,0}, {0.0, 0.0},
                                        {0.0, 0.0}, {0.70711,0}, {0.0, 0.0}, {0.70711,0},
                                        {0.70711,0}, {0.0, 0.0}, {0.70711,0}, {0.0, 0.0},
                                        {0.0, 0.0}, {0.70711,0}, {0.0, 0.0}, {0.70711,0}};
struct complex *R_sqrt_22_corr[1]     = {R_sqrt_22_corr_tap};

//correlation matrix for a fully correlated 2x1 channel (h1==h2)
struct complex R_sqrt_21_corr_tap[4]  = {{0.70711,0}, {0.70711,0}, {0.70711,0}, {0.70711,0}};
struct complex *R_sqrt_21_corr[1]      = {R_sqrt_21_corr_tap};

//correlation matrix for a 2x2 channel with full Tx anti-correlation
struct complex R_sqrt_22_anticorr_tap[16] = {{0.70711,0}, {0.0, 0.0}, {-0.70711,0}, {0.0, 0.0},
                                             {0.0, 0.0}, {0.70711,0}, {0.0, 0.0}, {-0.70711,0},
                                             {-0.70711,0}, {0.0, 0.0}, {0.70711,0}, {0.0, 0.0},
                                             {0.0, 0.0}, {-0.70711,0}, {0.0, 0.0}, {0.70711,0}};
struct complex *R_sqrt_22_anticorr[1]     = {R_sqrt_22_anticorr_tap};

//correlation matrix for a fully anti-correlated 2x1 channel (h1==-h2)
struct complex R_sqrt_21_anticorr_tap[4]  = {{0.70711,0}, {-0.70711,0}, {-0.70711,0}, {0.70711,0}};
struct complex *R_sqrt_21_anticorr[1]     = {R_sqrt_21_anticorr_tap};

struct complex **R_sqrt_ptr2;

// full correlation matrix in vectorized form for 2x2 channel, where h1 is  perfectly orthogonal to h2

struct complex R_sqrt_22_orthogonal_tap[16] = {{0.70711,0.0}, {0.0, 0.0}, {0.0,0.0}, {0.0, 0.0},
                                                        {0.0, 0.0}, {0.0,0.0}, {0.0, 0.0}, {0.0,0.0},
                                                        {0.0,0.0}, {0.0, 0.0}, {0.0,0.0}, {0.0, 0.0},
                                                        {0.0, 0.0}, {0.0,0.0}, {0.0, 0.0}, {0.70711,0.0}};
struct complex *R_sqrt_22_orthogonal[1]     = {R_sqrt_22_orthogonal_tap};

// full correlation matrix for TM4 to make orthogonal effective channel

struct complex R_sqrt_22_orth_eff_ch_TM4_prec_real_tap[16] = {{0.70711,0.0}, {0.0, 0.0}, {0.70711,0.0}, {0.0, 0.0},
                                                        {0.0, 0.0}, {0.70711,0.0}, {0.0, 0.0}, {-0.70711,0.0},
                                                        {0.70711,0.0}, {0.0, 0.0}, {0.70711,0.0}, {0.0, 0.0},
                                                        {0.0, 0.0}, {-0.70711,0.0}, {0.0, 0.0}, {0.70711,0.0}};
struct complex *R_sqrt_22_orth_eff_ch_TM4_prec_real[1]     = {R_sqrt_22_orth_eff_ch_TM4_prec_real_tap};

struct complex R_sqrt_22_orth_eff_ch_TM4_prec_imag_tap[16] = {{0.70711,0.0}, {0.0,0.0}, {0.0, -0.70711}, {0.0,0.0},
                                                        {0.0, 0.0}, {0.70711,0.0}, {0.0, 0.0}, {0.0,0.70711},
                                                        {0.0,-0.70711}, {0.0, 0.0}, {-0.70711,0.0}, {0.0, 0.0},
                                                        {0.0, 0.0}, {0.0,0.70711}, {0.0, 0.0}, {-0.70711,0.0}};
struct complex *R_sqrt_22_orth_eff_ch_TM4_prec_imag[1]     = {R_sqrt_22_orth_eff_ch_TM4_prec_imag_tap};

//Correlation matrix for EPA channel
struct complex R_sqrt_22_EPA_low_tap[16] = {{1.0,0.0}, {0.0,0.0}, {0.0,0.0}, {0.0,0.0},
                                            {0.0,0.0}, {1.0,0.0}, {0.0,0.0}, {0.0,0.0},
                                            {0.0,0.0}, {0.0,0.0}, {1.0,0.0}, {0.0,0.0},
                                            {0.0,0.0}, {0.0,0.0}, {0.0,0.0}, {1.0,0.0}};
struct complex *R_sqrt_22_EPA_low[1]     = {R_sqrt_22_EPA_low_tap};

struct complex R_sqrt_22_EPA_high_tap[16] = {{0.7179,0.0}, {0.4500,0.0}, {0.4500,0.0}, {0.2821,0.0},
                                             {0.4500,0.0}, {0.7179,0.0}, {0.2821,0.0}, {0.4500,0.0},
                                             {0.4500,0.0}, {0.2821,0.0}, {0.7179,0.0}, {0.4500,0.0},
                                             {0.2821,0.0}, {0.4500,0.0}, {0.4500,0.0}, {0.7179,0.0}};
struct complex *R_sqrt_22_EPA_high[1]     = {R_sqrt_22_EPA_high_tap};

struct complex R_sqrt_22_EPA_medium_tap[16] = {{0.8375,0.0}, {0.5249,0.0}, {0.1286,0.0}, {0.0806,0.0},
                                               {0.5249,0.0}, {0.8375,0.0}, {0.0806,0.0}, {0.1286,0.0},
                                               {0.1286,0.0}, {0.0806,0.0}, {0.8375,0.0}, {0.5249,0.0},
                                               {0.0806,0.0}, {0.1286,0.0}, {0.5249,0.0}, {0.8375,0.0}};
struct complex *R_sqrt_22_EPA_medium[1]     = {R_sqrt_22_EPA_medium_tap};

//Rayleigh1_orth_eff_ch_TM4

channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
                                     uint8_t nb_rx,
                                     SCM_t channel_model,
                                     double sampling_rate,
                                     double channel_bandwidth,
                                     double forgetting_factor,
                                     int32_t channel_offset,
                                     double path_loss_dB)
{

  channel_desc_t *chan_desc = (channel_desc_t *)malloc(sizeof(channel_desc_t));
  uint16_t i,j;
  double sum_amps;
  double aoa,ricean_factor,Td,maxDoppler;
  int channel_length,nb_taps;

  chan_desc->nb_tx                      = nb_tx;
  chan_desc->nb_rx                      = nb_rx;
  chan_desc->sampling_rate              = sampling_rate;
  chan_desc->channel_bandwidth          = channel_bandwidth;
  chan_desc->forgetting_factor          = forgetting_factor;
  chan_desc->channel_offset             = channel_offset;
  chan_desc->path_loss_dB               = path_loss_dB;
  chan_desc->first_run                  = 1;
  chan_desc->ip                                 = 0.0;

  LOG_I(OCM,"Channel Model (inside of new_channel_desc_scm)=%d\n\n", channel_model);

  switch (channel_model) {
  case SCM_A:
    LOG_W(OCM,"channel model not yet supported\n");
    free(chan_desc);
    return(NULL);
  case SCM_B:
    LOG_W(OCM,"channel model not yet supported\n");
    free(chan_desc);
    return(NULL);
  case SCM_C:
    chan_desc->nb_taps        = 18;
    chan_desc->Td             = 4.625;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*scm_c_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = scm_c_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));

    chan_desc->R_sqrt  = (struct complex**) malloc(6*sizeof(struct complex**));
    if (nb_tx==2 && nb_rx==2) {
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
    }
    else if (nb_tx==2 && nb_rx==1) {
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R21_sqrt[i][0];
    }
    else if (nb_tx==1 && nb_rx==2) {
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R12_sqrt[i][0];
    }
    else {
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix not implemented for nb_tx==%d and nb_rx==%d, using identity\n", nb_tx, nb_rx);
      }
    }
    break;
  case SCM_D:
    LOG_W(OCM,"This is not the real SCM-D model! It is just SCM-C with an additional Rice factor!\n");
    chan_desc->nb_taps        = 18;
    chan_desc->Td             = 4.625;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*scm_c_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = scm_c_delays;
    chan_desc->ricean_factor  = 0.1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));

    chan_desc->R_sqrt  = (struct complex**) malloc(6*sizeof(struct complex**));
    if (nb_tx==2 && nb_rx==2) {
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
    }
    else if (nb_tx==2 && nb_rx==1) {
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R21_sqrt[i][0];
    }
    else if (nb_tx==1 && nb_rx==2) {
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R12_sqrt[i][0];
    }
    else {
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix not implemented for nb_tx==%d and nb_rx==%d, using identity\n", nb_tx, nb_rx);
      }
    }
    break;
  case EPA:
    chan_desc->nb_taps        = 7;
    chan_desc->Td             = .410;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*epa_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = epa_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
    if (nb_tx==2 && nb_rx==2) {
      chan_desc->R_sqrt  = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
    }
    else {
      chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
      }
    }
    break;
  case EPA_low:
    chan_desc->nb_taps        = 7;
    chan_desc->Td             = .410;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*epa_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = epa_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
    if (nb_tx==2 && nb_rx==2) {
      chan_desc->R_sqrt  = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex**));
      for (i = 0; i<chan_desc->nb_taps; i++)
        chan_desc->R_sqrt[i] = R_sqrt_22_EPA_low[0];
    }
      else {
      printf("Correlation matrices are implemented for 2 x 2 only");
    }
    /*else {
      chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
      }
    }*/
    break;
  case EPA_high:
    chan_desc->nb_taps        = 7;
    chan_desc->Td             = .410;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*epa_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = epa_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
    if (nb_tx==2 && nb_rx==2) {
      chan_desc->R_sqrt  = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex**));
      for (i = 0; i<chan_desc->nb_taps; i++)
        chan_desc->R_sqrt[i] = R_sqrt_22_EPA_high[0];
    }
    else {
      printf("Correlation matrices are implemented for 2 x 2 only");
    }
    /*else {
      chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
      }
    }*/
    break;
    case EPA_medium:
    chan_desc->nb_taps        = 7;
    chan_desc->Td             = .410;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*epa_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = epa_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
    if (nb_tx==2 && nb_rx==2) {
      chan_desc->R_sqrt  = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex**));
      for (i = 0; i<chan_desc->nb_taps; i++)
        chan_desc->R_sqrt[i] = R_sqrt_22_EPA_medium[0];
    } else {
      printf("Correlation matrices are implemented for 2 x 2 only");
    }
    /*else {
      chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
      }
    }*/
    break;
  case EVA:
    chan_desc->nb_taps        = 9;
    chan_desc->Td             = 2.51;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*eva_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = eva_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
    if (nb_tx==2 && nb_rx==2) {
      chan_desc->R_sqrt  = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
    }
    else {
      chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
      }
    }
    break;
  case ETU:
    chan_desc->nb_taps        = 9;
    chan_desc->Td             = 5.0;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*etu_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = etu_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
    if (nb_tx==2 && nb_rx==2) {
      chan_desc->R_sqrt  = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++)
        chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
    }
    else {
      chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
      for (i = 0; i<6; i++) {
        chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
        for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
          chan_desc->R_sqrt[i][j].x = 1.0;
          chan_desc->R_sqrt[i][j].y = 0.0;
        }
        LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
      }
    }
    break;
  case MBSFN:
    chan_desc->nb_taps        = 18;
    chan_desc->Td             = 28.58;
    chan_desc->channel_length = (int) (2*chan_desc->sampling_rate*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->sampling_rate*chan_desc->Td));
    sum_amps = 0;
    chan_desc->amps           = (double*) malloc(chan_desc->nb_taps*sizeof(double));
    for (i = 0; i<chan_desc->nb_taps; i++) {
      chan_desc->amps[i]      = pow(10,.1*mbsfn_amps_dB[i]);
      sum_amps += chan_desc->amps[i];
    }
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->amps[i] /= sum_amps;
    chan_desc->delays         = mbsfn_delays;
    chan_desc->ricean_factor  = 1;
    chan_desc->aoa            = 0;
    chan_desc->random_aoa     = 0;
    chan_desc->ch             = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->chF            = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
    chan_desc->a              = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
    for (i = 0; i<nb_tx*nb_rx; i++)
      chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
    for (i = 0; i<chan_desc->nb_taps; i++)
      chan_desc->a[i]         = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));

    chan_desc->R_sqrt  = (struct complex**) malloc(6*sizeof(struct complex*));
    for (i = 0; i<6; i++) {
      chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
      for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
        chan_desc->R_sqrt[i][j].x = 1.0;
        chan_desc->R_sqrt[i][j].y = 0.0;
      }
      LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
    }
    break;
  case Rayleigh8:

      nb_taps = 8;
      Td = 0.8;
      channel_length = (int)11+2*sampling_rate*Td;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,
                        nb_tx,
                        nb_rx,
                        nb_taps,
                        channel_length,
                        default_amps_lin,
                        NULL,
                        NULL,
                        Td,
                        sampling_rate,
                        channel_bandwidth,
                        ricean_factor,
                        aoa,
                        forgetting_factor,
                        maxDoppler,
                        channel_offset,
                        path_loss_dB,
                        0);
      break;

  case Rice8:
      nb_taps = 8;
      Td = 0.8;
      channel_length = (int)11+2*sampling_rate*Td;
      ricean_factor = 0.1;
      aoa = 0.7854;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amps_lin,
                                   NULL,
                                   NULL,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   1);
      break;

  case Rayleigh1:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   NULL,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

  case Rayleigh1_800:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 800;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   NULL,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

  case Rayleigh1_corr:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==1)) {
        R_sqrt_ptr2 = R_sqrt_21_corr;
      }
      else if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_corr;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

  case Rayleigh1_anticorr:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==1)) { //check this
        R_sqrt_ptr2 = R_sqrt_21_anticorr;
      }
      else if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_anticorr;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

  case Rice1:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 0.1;
      aoa = 0.7854;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   NULL,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

  case AWGN:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 0.0;
      aoa = 0.0;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   NULL,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      printf("AWGN: ricean_factor %f\n",chan_desc->ricean_factor);

      break;

  case TS_SHIFT:
      nb_taps = 2;
      Td = ts_shift_delays[1];
      channel_length = 10;
      ricean_factor = 0.0;
      aoa = 0.0;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   ts_shift_amps,
                                   ts_shift_delays,
                                   NULL,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      printf("TS_SHIFT: ricean_factor %f\n",chan_desc->ricean_factor);

      break;

  case Rice1_corr:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 0.1;
      aoa = .03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==1)) {
        R_sqrt_ptr2 = R_sqrt_21_corr;
      }
      else if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_corr;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   1);
      break;

  case Rice1_anticorr:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 0.1;
      aoa = .03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==1)) {
        R_sqrt_ptr2 = R_sqrt_21_anticorr;
      }
      else if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_anticorr;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   1);
      break;

  case Rayleigh1_orthogonal:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = 0.03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_orthogonal;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

  case Rayleigh1_orth_eff_ch_TM4_prec_real:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = 0.03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_real;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   1);
      break;

      case Rayleigh1_orth_eff_ch_TM4_prec_imag:
      nb_taps = 1;
      Td = 0;
      channel_length = 1;
      ricean_factor = 1;
      aoa = 0.03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_imag;
      }
      else
        R_sqrt_ptr2 = NULL;

      fill_channel_desc(chan_desc,nb_tx,
                                   nb_rx,
                                   nb_taps,
                                   channel_length,
                                   default_amp_lin,
                                   NULL,
                                   R_sqrt_ptr2,
                                   Td,
                                   sampling_rate,
                                   channel_bandwidth,
                                   ricean_factor,
                                   aoa,
                                   forgetting_factor,
                                   maxDoppler,
                                   channel_offset,
                                   path_loss_dB,
                                   0);
      break;

       case Rayleigh8_orth_eff_ch_TM4_prec_real:

    if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_real;
        //R_sqrt_ptr2 = NULL;
      }
      else
        R_sqrt_ptr2 = NULL;

      nb_taps = 8;
      Td = 0.8;
      channel_length = (int)11+2*sampling_rate*Td;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 0;

      fill_channel_desc(chan_desc,
                        nb_tx,
                        nb_rx,
                        nb_taps,
                        channel_length,
                        default_amps_lin,
                        NULL,
                        R_sqrt_ptr2,
                        Td,
                        sampling_rate,
                        channel_bandwidth,
                        ricean_factor,
                        aoa,
                        forgetting_factor,
                        maxDoppler,
                        channel_offset,
                        path_loss_dB,
                        0);

      break;

      case Rayleigh8_orth_eff_ch_TM4_prec_imag:
      nb_taps = 8;
      Td = 0.8;
      channel_length = (int)11+2*sampling_rate*Td;
      ricean_factor = 1;
      aoa = .03;
      maxDoppler = 0;

      if ((nb_tx==2) && (nb_rx==2)) {
        R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_imag;
      }
      else
        R_sqrt_ptr2 = NULL;

       fill_channel_desc(chan_desc,
                        nb_tx,
                        nb_rx,
                        nb_taps,
                        channel_length,
                        default_amps_lin,
                        NULL,
                        R_sqrt_ptr2,
                        Td,
                        sampling_rate,
                        channel_bandwidth,
                        ricean_factor,
                        aoa,
                        forgetting_factor,
                        maxDoppler,
                        channel_offset,
                        path_loss_dB,
                        0);
      break;

  default:
    LOG_W(OCM,"channel model not yet supported\n");
    free(chan_desc);
    return(NULL);
  }
  LOG_D(OCM,"[CHANNEL] RF %f\n",chan_desc->ricean_factor);
  for (i=0;i<chan_desc->nb_taps;i++)
    LOG_D(OCM,"[CHANNEL] tap %d: amp %f, delay %f\n",i,chan_desc->amps[i],chan_desc->delays[i]);

  chan_desc->nb_paths = 10;

  return(chan_desc);
}

int random_channel(channel_desc_t *desc, uint8_t abstraction_flag) {

  double s;
  int i,k,l,aarx,aatx;
  struct complex anew[NB_ANTENNAS_TX*NB_ANTENNAS_RX],acorr[NB_ANTENNAS_TX*NB_ANTENNAS_RX];
  struct complex phase, alpha, beta;

  AssertFatal(desc->nb_tx<=NB_ANTENNAS_TX && desc->nb_rx <= NB_ANTENNAS_RX,
	      "random_channel.c: Error: temporary buffer for channel not big enough (%d,%d)\n",desc->nb_tx,desc->nb_rx);

  start_meas(&desc->random_channel);
  for (i=0;i<(int)desc->nb_taps;i++) {
    for (aarx=0;aarx<desc->nb_rx;aarx++) {
      for (aatx=0;aatx<desc->nb_tx;aatx++) {

        anew[aarx+(aatx*desc->nb_rx)].x = sqrt(desc->ricean_factor*desc->amps[i]/2) * gaussdouble(0.0,1.0);
        anew[aarx+(aatx*desc->nb_rx)].y = sqrt(desc->ricean_factor*desc->amps[i]/2) * gaussdouble(0.0,1.0);

        if ((i==0) && (desc->ricean_factor != 1.0)) {
          if (desc->random_aoa==1) {
            desc->aoa = uniformrandom()*2*M_PI;
          }

          // this assumes that both RX and TX have linear antenna arrays with lambda/2 antenna spacing.
          // Furhter it is assumed that the arrays are parallel to each other and that they are far enough apart so
          // that we can safely assume plane wave propagation.
          phase.x = cos(M_PI*((aarx-aatx)*sin(desc->aoa)));
          phase.y = sin(M_PI*((aarx-aatx)*sin(desc->aoa)));

          anew[aarx+(aatx*desc->nb_rx)].x += phase.x * sqrt(1.0-desc->ricean_factor);
          anew[aarx+(aatx*desc->nb_rx)].y += phase.y * sqrt(1.0-desc->ricean_factor);
        }
#ifdef DEBUG_CH
        printf("(%d,%d,%d) %f->(%f,%f) (%f,%f) phase (%f,%f)\n",aarx,aatx,i,desc->amps[i],anew[aarx+(aatx*desc->nb_rx)].x,anew[aarx+(aatx*desc->nb_rx)].y,desc->aoa,desc->ricean_factor,phase.x,phase.y);
#endif
      } //aatx
    } //aarx

    /*
    // for debugging set a=anew;
    for (aarx=0;aarx<desc->nb_rx;aarx++) {
      for (aatx=0;aatx<desc->nb_tx;aatx++) {
        desc->a[i][aarx+(aatx*desc->nb_rx)].x = anew[aarx+(aatx*desc->nb_rx)].x;
        desc->a[i][aarx+(aatx*desc->nb_rx)].y = anew[aarx+(aatx*desc->nb_rx)].y;
        printf("anew(%d,%d) = %f+1j*%f\n",aatx,aarx,anew[aarx+(aatx*desc->nb_rx)].x, anew[aarx+(aatx*desc->nb_rx)].y);
     }
    }
    */
    //apply correlation matrix
    //compute acorr = R_sqrt[i] * anew
    alpha.x = 1.0;
    alpha.y = 0.0;
    beta.x = 0.0;
    beta.y = 0.0;

    cblas_zgemv(CblasRowMajor, CblasNoTrans, desc->nb_tx*desc->nb_rx, desc->nb_tx*desc->nb_rx,
                (void*) &alpha, (void*) desc->R_sqrt[i/3], desc->nb_rx*desc->nb_tx,
                (void*) anew, 1, (void*) &beta, (void*) acorr, 1);

    /*
    for (aarx=0;aarx<desc->nb_rx;aarx++) {
      for (aatx=0;aatx<desc->nb_tx;aatx++) {
        desc->a[i][aarx+(aatx*desc->nb_rx)].x = acorr[aarx+(aatx*desc->nb_rx)].x;
        desc->a[i][aarx+(aatx*desc->nb_rx)].y = acorr[aarx+(aatx*desc->nb_rx)].y;
        printf("tap %d, acorr1(%d,%d) = %f+1j*%f\n",i,aatx,aarx,acorr[aarx+(aatx*desc->nb_rx)].x, acorr[aarx+(aatx*desc->nb_rx)].y);
      }
    }
    */

    if (desc->first_run==1){
      cblas_zcopy(desc->nb_tx*desc->nb_rx, (void*) acorr, 1, (void*) desc->a[i], 1);
    }
    else {
      // a = alpha*acorr+beta*a
      // a = beta*a
      // a = a+alpha*acorr
      alpha.x = sqrt(1-desc->forgetting_factor);
      alpha.y = 0;
      beta.x = sqrt(desc->forgetting_factor);
      beta.y = 0;
      cblas_zscal(desc->nb_tx*desc->nb_rx, (void*) &beta, (void*) desc->a[i], 1);
      cblas_zaxpy(desc->nb_tx*desc->nb_rx, (void*) &alpha, (void*) acorr, 1, (void*) desc->a[i], 1);

      //  desc->a[i][aarx+(aatx*desc->nb_rx)].x = (sqrt(desc->forgetting_factor)*desc->a[i][aarx+(aatx*desc->nb_rx)].x) + sqrt(1-desc->forgetting_factor)*anew.x;
      //  desc->a[i][aarx+(aatx*desc->nb_rx)].y = (sqrt(desc->forgetting_factor)*desc->a[i][aarx+(aatx*desc->nb_rx)].y) + sqrt(1-desc->forgetting_factor)*anew.y;
    }

    /*
    for (aarx=0;aarx<desc->nb_rx;aarx++) {
      for (aatx=0;aatx<desc->nb_tx;aatx++) {
        //desc->a[i][aarx+(aatx*desc->nb_rx)].x = acorr[aarx+(aatx*desc->nb_rx)].x;
        //desc->a[i][aarx+(aatx*desc->nb_rx)].y = acorr[aarx+(aatx*desc->nb_rx)].y;
        printf("tap %d, a(%d,%d) = %f+1j*%f\n",i,aatx,aarx,desc->a[i][aarx+(aatx*desc->nb_rx)].x, desc->a[i][aarx+(aatx*desc->nb_rx)].y);
      }
    }
    */

  } //nb_taps
  stop_meas(&desc->random_channel);

  //memset((void *)desc->ch[aarx+(aatx*desc->nb_rx)],0,(int)(desc->channel_length)*sizeof(struct complex));

  if (abstraction_flag==0) {
    start_meas(&desc->interp_time);

    for (aarx=0; aarx<desc->nb_rx; aarx++) {
      for (aatx=0; aatx<desc->nb_tx; aatx++) {
        if (desc->channel_length == 1) {
          desc->ch[aarx+(aatx*desc->nb_rx)][0].x = desc->a[0][aarx+(aatx*desc->nb_rx)].x;
          desc->ch[aarx+(aatx*desc->nb_rx)][0].y = desc->a[0][aarx+(aatx*desc->nb_rx)].y;
        } else {

          for (k=0; k<(int)desc->channel_length; k++) {
            desc->ch[aarx+(aatx*desc->nb_rx)][k].x = 0.0;
            desc->ch[aarx+(aatx*desc->nb_rx)][k].y = 0.0;

            for (l=0; l<desc->nb_taps; l++) {
              if ((k - (desc->delays[l]*desc->sampling_rate) - NB_SAMPLES_CHANNEL_OFFSET) == 0)
                s = 1.0;
              else
                s = sin(M_PI*(k - (desc->delays[l]*desc->sampling_rate) - NB_SAMPLES_CHANNEL_OFFSET))/
                    (M_PI*(k - (desc->delays[l]*desc->sampling_rate) - NB_SAMPLES_CHANNEL_OFFSET));

              desc->ch[aarx+(aatx*desc->nb_rx)][k].x += s*desc->a[l][aarx+(aatx*desc->nb_rx)].x;
              desc->ch[aarx+(aatx*desc->nb_rx)][k].y += s*desc->a[l][aarx+(aatx*desc->nb_rx)].y;
              //    printf("l %d : desc->ch.x %f\n",l,desc->a[l][aarx+(aatx*desc->nb_rx)].x);

            } //nb_taps

#ifdef DEBUG_CH
        k=0;
        printf("(%d,%d,%d)->(%f,%f)\n",k,aarx,aatx,desc->ch[aarx+(aatx*desc->nb_rx)][k].x,desc->ch[aarx+(aatx*desc->nb_rx)][k].y);
#endif
        }
      } //channel_length
    } //aatx
  } //aarx
  stop_meas(&desc->interp_time);
  }

  if (desc->first_run==1)
    desc->first_run = 0;

  return (0);
}

double N_RB2sampling_rate(uint16_t N_RB)
{
  double sampling_rate;
  switch (N_RB) {
  case 6:
    sampling_rate = 1.92;
    break;

  case 25:
    sampling_rate = 7.68;
    break;

  case 50:
    sampling_rate = 15.36;
    break;

  case 100:
    sampling_rate = 30.72;
    break;

  default:
    AssertFatal(1==0,"Unknown N_PRB %d",N_RB);
  }

  return(sampling_rate);
}

double N_RB2channel_bandwidth(uint16_t N_RB)
{
  double channel_bandwidth;
  switch (N_RB) {
  case 6:
    channel_bandwidth = 1.25;
    break;

  case 25:
    channel_bandwidth = 5.00;
    break;

  case 50:
    channel_bandwidth = 10.00;
    break;

  case 100:
    channel_bandwidth = 20.00;
    break;

  default:
    LOG_E(PHY,"Unknown N_PRB\n");
    return(-1);
  }
  return(channel_bandwidth);
}

#ifdef RANDOM_CHANNEL_MAIN
#define sampling_rate 5.0
#define Td 2.0
main(int argc,char **argv) {

  double amps[8] = {.8,.2,.1,.04,.02,.01,.005};
  struct complex ch[(int)(1+2*sampling_rate*Td)],phase;
  int i;

  randominit();
  phase.x = 1.0;
  phase.y = 0;
  random_channel(amps,Td, 8,sampling_rate,ch,(double)1.0,&phase);
  /*
  for (i=0;i<(11+2*sampling_rate*Td);i++){
    printf("%f + sqrt(-1)*%f\n",ch[i].x,ch[i].y);
  }
  */
}

#endif

max_stars_repo_name: Imyne/opran

max_issues_repo_name: Imyne/opran

max_forks_repo_name: Imyne/opran