close all; clear all; clc;

% Number of samples to generate
Nsamps = 100000;

% Number of signals
n = 4;

% Set the required second order statistics
P1 = 2.3;
P2 = 0.75;
P3 = 3.4;
P4 = 1.23;
rho12 = 0.2 - 1j*0.3;
rho13 = -0.6 + 1j*0.1;
rho14 = -1j*0.4;
rho23 = 0.1 + 1j*0.1;
rho24 = 0.5;
rho34 = -0.3 - 1j*0.1;

% 4x4 (theoretical) covariance matrix
Rxx = [P1,                      sqrt(P1*P2)*rho12,       sqrt(P1*P3)*rho13,       sqrt(P1*P4)*rho14;
       sqrt(P2*P1)*conj(rho12), P2,                      sqrt(P2*P3)*rho23,       sqrt(P2*P4)*rho24;
       sqrt(P3*P1)*conj(rho13), sqrt(P3*P2)*conj(rho23), P3,                      sqrt(P3*P4)*rho34;
       sqrt(P4*P1)*conj(rho14), sqrt(P4*P2)*conj(rho24), sqrt(P4*P3)*conj(rho34), P4]

% Generate zero-mean, unit-power, uncorrelated Gaussian random signals
s = randn(n, Nsamps);

% Create partially correlated signals
L = chol(Rxx,'lower');
X = L*randn(n,Nsamps);

% 4x4 (practical) sample covariance matrix
Rxx_samp = X*X' / Nsamps