LTENR_Fading.c

/************************************************************************************
* Copyright (C) 2023                                                                *
* TETCOS, Bangalore. India                                                          *
* This source code is licensed per the NetSim license agreement.                    *
* No portion of this source code may be used as the basis for a derivative work.    *
*                                                                                   *
* Author: Suraj(Rician) && Charu(Rayleigh)                                          *
************************************************************************************/
 
#include <stdlib.h>
#include <math.h>
#include "NetSim_utility.h"
#include "LTENR_Fading.h"
#include "lapacke.h"
#include "main.h"
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
#define MAX_RETRIES 10
#define MIN(a, b) ((a) < (b) ? (a) : (b))
 
/* Box–Muller Gaussian */
static double gaussian_random(void) {
    double U1 = NETSIM_RAND_01();
    double U2 = NETSIM_RAND_01();
    return sqrt(-2.0 * log(U1)) * cos(2.0 * M_PI * U2);
}
 
/* Descending sort comparator */
static int cmp_desc(const void* a, const void* b) {
    double da = *(double*)a, db = *(double*)b;
    if (da < db) return 1;
    if (da > db) return -1;
    return 0;
}
 
/* Rayleigh H entries */
static void generate_rayleigh_H(lapack_complex_double* H, UINT NR, UINT NT) {
    for (UINT i = 0; i < NR * NT; i++) {
        H[i].real = gaussian_random() / sqrt(2.0);
        H[i].imag = gaussian_random() / sqrt(2.0);
    }
}
 
/* Rician H entries */
static void generate_rician_H(lapack_complex_double* H, UINT NR, UINT NT, double K) {
    double mu = sqrt(K / (2.0 * (K + 1.0)));
    double sigma = sqrt(1.0 / (2.0 * (K + 1.0)));
    for (UINT i = 0; i < NR * NT; i++) {
        H[i].real = mu + sigma * gaussian_random();
        H[i].imag = mu + sigma * gaussian_random();
    }
}
 
/* Unified Rayleigh/Rician eigen-solver */
void LTENR_Beamforming_GetValue(
    UINT NT, UINT NR,
    UINT fastFadingModel,
    double K_factor,
    double** eigenvalues)
{
    UINT layers = MIN(NT, NR);
    double* results = malloc(layers * sizeof(*results));
    if (!results) { *eigenvalues = NULL; return; }
 
    lapack_complex_double* H = malloc(NR * NT * sizeof(*H));
    lapack_complex_double* G = malloc(NR * NR * sizeof(*G));
    double* w = malloc(NR * sizeof(*w));
    if (!H || !G || !w) {
        free(results); free(H); free(G); free(w);
        *eigenvalues = NULL; return;
    }
 
    lapack_int info = -1;
    int retry = 0;
    while (retry < MAX_RETRIES && info != 0) {
        if (fastFadingModel == FADING_MODEL_RICIAN) {
            generate_rician_H(H, NR, NT, K_factor);
        }
        else if (fastFadingModel == FADING_MODEL_RAYLEIGH) {
            // Default to Rayleigh for any other model
            generate_rayleigh_H(H, NR, NT);
        }
 
        /* Build G = H * H^H */
        for (UINT i = 0; i < NR; i++) {
            for (UINT j = 0; j < NR; j++) {
                lapack_complex_double sum = { 0.0, 0.0 };
                for (UINT k = 0; k < NT; k++) {
                    double ar = H[i * NT + k].real, ai = H[i * NT + k].imag;
                    double br = H[j * NT + k].real, bi = -H[j * NT + k].imag;
                    sum.real += ar * br - ai * bi;
                    sum.imag += ar * bi + ai * br;
                }
                G[i * NR + j] = sum;
            }
        }
 
        info = LAPACKE_zheevd(LAPACK_ROW_MAJOR, 'N', 'U', NR, G, NR, w);
        if (info != 0) retry++;
    }
 
    if (info != 0) {
        for (UINT i = 0; i < layers; i++) results[i] = 0.0;
    }
    else {
        qsort(w, NR, sizeof(*w), cmp_desc);
        for (UINT i = 0; i < layers; i++) results[i] = w[i];
    }
 
    for (int i = 0; i < layers; i++) {
        double v = results[i];
        if (v < 1e-12) v = 1e-12;      // avoid log(0)
        results[i] = 10.0 * log10(v);  // convert to dB
    }
 
    for (UINT i = 0; i < layers; i++) {
        (*eigenvalues)[i] = results[i];
    }
    free(results);
    free(H); free(G); free(w);
}