Cholesky factorization
%Matrix A here is needed for factorization (A = L' * D * L)
%clear console; Remove all variables from memory; Close all figures
clc;clearvars;close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%All_Possible_Routes_On_Tree%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TxRx = 4; %Number(Equal) of transmitters and Recievers %Number_Of_Levels_On_Tree
iter = 1; %Number of iterations/data signals
PTX = 10;
%Complex Gaussian Noise Signal, Noise
Cn = [4 1 0 2; 1 6 1 1; 0 1 5 2; 2 1 2 10];
[V,D] = eig(Cn);
squarerootCn = V*sqrt(D)*V';
Noise = 1/sqrt(2) * squarerootCn * (randn(TxRx,iter) + 1i*randn(TxRx,iter));
%Normal Distribution Matrix, H
mean = 0; variance = 1;
a = randn(TxRx) + 1i*randn(TxRx);
H = mean + sqrt(variance)*a;
%With PTX in block diagram
Hptx = sqrt(PTX/TxRx) * H;
%Random Signal (WO influence of Noise)
Integers = [0 1 2 3];
Constellations = [exp(1i*(pi/4)) exp(1i*(3.*pi/4)) exp(1i*(5.*pi/4)) exp(1i*(7.*pi/4))];
RandomIntegers = randi([0 3],iter,TxRx);
RandomConstellations = interp1(Integers, Constellations, RandomIntegers);
Sq = transpose(RandomConstellations);
%A = H' * Cn^-1 * H
A = H' * Cn^-1 * H
%A is non-negative definite if eigenvalues of A are non-negative.
eigenvalues_of_A = eig(A); %Checked: A is non-negative definite
%Cholesky factorization A = L' * D * L
[L, D] = ldl_golub(A); % but achieved A = L * D * L'
function [L, D] = ldl_golub(A)
n = length(A);
v = zeros(n,1);
for j = 1:n
for i = 1:j - 1
v(i) = A(j, i)*A(i, i);
end
A(j, j) = A(j, j) - A(j, 1:j - 1) * v(1:j - 1);
A(j + 1:n, j) = (A(j + 1:n,j) - A(j + 1:n, 1:j - 1) * v(1:j - 1))/A(j, j);
end
L = tril(A,-1)+eye(n);
D = diag(diag(A));
end