% estimate_entropies.m
clear; rng(0);

%% 1) Parameters
M        = 1e5;          % number of Monte-Carlo samples
EbNo_dB  = 0;           % SNR per bit in dB
EbNo     = 10^(EbNo_dB/10);
sigma2   = 1/(2*EbNo);   % noise variance per real dimension

% 4-QAM constellation (row vector)
X  = [1+1j, 1-1j, -1+1j, -1-1j];  
N  = numel(X);
PX = ones(1,N)/N;        % uniform PMF

%% 2) Generate transmit symbols and AWGN
idx = randi(N,1,M);             % 1×M random symbol indices
s   = X(idx);                   % 1×M transmitted symbols
n   = sqrt(sigma2)*(randn(1,M) + 1j*randn(1,M));
y   = s + n;                    % 1×M received samples

%% 3) Compute p_{Y|X}(y|x) for each constellation point
%    Create an N×M matrix where row n is |y - X(n)|^2
d2     = abs(bsxfun(@minus, X(:), y)).^2;        % N×M
pYgX   = (1/(pi*sigma2)) * exp(-d2 / sigma2);    % N×M

%% 4) Estimate H(Y) = -E[ log2 p_Y(Y) ]
%    Mixture density p_Y(y_m) = sum_n PX(n)*pYgX(n,m)
pY     = PX * pYgX;                 % 1×M
HY      = -mean(log2(pY));          % bits

%% 5) Estimate H(Y|X) = -E[ log2 p(Y|X) ]
%    For each m, pick the row corresponding to the true idx(m)
linearIdx = sub2ind([N, M], idx, 1:M);
pYgX_true = pYgX(linearIdx);        % 1×M
HYgX     = -mean(log2(pYgX_true));  % bits

%% 6) Mutual information
I = HY - HYgX;

%% 7) Display
fprintf('Estimated H(Y)      = %.4f bits\n', HY);
fprintf('Estimated H(Y|X)    = %.4f bits\n', HYgX);
fprintf('Estimated I(X;Y)    = %.4f bits\n', I);