function  avg_prob_det = PdSNR_RealData(m_tx, p_rx, K, lp_class, Prob_fa, snr_db_vec, rho, pow_gen)
% avg_prob_det = PdSNR_RealData(m_tx, p_rx, K, lp_class, Prob_fa, snr_db_vec, rho, pow_gen)
%
% This code can be used to generate the probability of detection vs SNR
% curve for the U_{BF} Eq. (24) detector in Fig. 1 of the paper listed
% below, where the data is real.
%
% Input
%   m_tx: scalar, number of Tx antennas.
%   p_rx: scalar, number of Rx antennas.
%   K: scalar, number of IID received data samples.
%   lp_class: scalar, objective prior iinducator:
%               1 for Jeffreys Prior,
%               2 for Independence Jeffreys Prior,
%               3 for Geisser and Cornfield Prior,
%               4 for Right Haar Measure,
%               5 for Reference prior.
%   Prob_fa: scalar, probability of false alarm.
%   snr_db_vec: (Q x 1), Q SNR values in log-scale (20 log_10(SNR)).
%   rho: scalar, correlation parameter to produce Eq. (40) in the reference paper.
%   pow_gen:  (p_rx x 1), diagonal values of Omega_0 in the log scale (20 log_10(pow_gen)).
%
% Output
%   avg_prob_det: (Q x 1), Probability of detection for the Q SNR values.
%
% Reference:
% M. H. Al-Ali and K.C. Ho, "Objective Bayesian approach for binary hypothesis
% testing of multivariate Gaussian observations," IEEE Trans. Inf. Theory,
% vol. 69, no. 2, pp. 1337-1354, Feb. 2023.
%
% M. H. Al-Ali and K. C. Ho, Dec. 12. 2022
%
%       Copyright (C) 2022
%       Compputer Engineering Department
%       University of Basrah
%       Basrah Province, Iraq
%       mohannd.khalaf@uobasrah.edu.iq
%
%       Copyright (C) 2022
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

% -- Objective prior parameters --
eps_val = .01;              % A small parameter to satisfy eq. (27) in the reference paper
m = 1:p_rx;
if lp_class == 1
    a_m = p_rx - m;         % Jeffreys Prior
    h = 1/2 + eps_val;
elseif lp_class == 2
    a_m = p_rx - m + 1;     % Independence Jeffreys Prior
    h = 0 + eps_val;

elseif lp_class == 3
    a_m = 2 - m;            % Geisser and Cornfield Prior
    h = (p_rx - 1)/2 + eps_val;
elseif lp_class == 4
    a_m = m;                % Right Haar Measure
    h = 0 + eps_val;
elseif lp_class == 5
    a_m = ones(1, p_rx);    % Reference Prior
    h = 0 + eps_val;
end

% -- Gamma distribution parameters --
h_tilde =  h + K/2;
delta_vec = 1e-4*ones(p_rx, 1);                 % Delta values for the gamma distribution

alf = (K - a_m + 1)/2;                          % The exponent of the denominator in the Bayesian detector
alf_max = max(max(alf), h_tilde);               % Used to scale exponents of the detector

% -- Simulation parameters --
num_rlzn = 1e5;                                   % Number of realizations
snr_length = length(snr_db_vec);
Prob_det = zeros(num_rlzn, snr_length);

% -- Covariance matrix for null hypothesis --
Cov_n = zeros(p_rx, p_rx);
for i = 1: p_rx
    for j = 1: p_rx
        if i <= j
            Cov_n(i, j) = rho^(j - i);
        else
            Cov_n(i, j) = rho^( i - j);
        end
    end
end
pow_noisedB = pow_gen(1: p_rx);
pow_noise = 10.^(pow_noisedB/10);
D = diag(pow_noise);
Omega_0 = (D.^(.5)*Cov_n*D.^(.5));

% -- Priors parameters for null hypothesis --
tau = 500;              % Large value for the parameter tau
F = zeros(p_rx, p_rx);

for k = 2:p_rx
    for b = 1: k-1
        F(b,k) = tau;
    end
end

F_cell = cell(p_rx - 1, 1);
D_cell = cell(p_rx - 1, 1);

for k = 2:p_rx
    F_cell{k} = diag(F(1:(k-1),k));
    D_cell{k} = eye(k-1);
end

% -- Detection threshold at specified Pfa --
[Th_val, ~, ~] = DtnTh_RealData(p_rx, delta_vec, Omega_0, F_cell, D_cell, Prob_fa, K, h, num_rlzn, lp_class);

% -- Simulation for obtaining Pd --
H = randn(p_rx, m_tx, num_rlzn);     % Generate channel matrices for all the realizations.

for lp_real = 1:num_rlzn,

    Cov_H = H(:,:,lp_real)*H(:,:,lp_real)';
    normalize_cov_H = Cov_H/trace(Cov_H);

    for lp_snr = 1:snr_length,

        % Generate model data
        snr_linear = 10^(snr_db_vec(lp_snr)/10);
        sig_pow = snr_linear*trace(Omega_0);
        Omega_1 = sig_pow*normalize_cov_H + Omega_0;
        Y = mvnrnd(zeros(1, p_rx), Omega_1, K);     % Generates an K-by-p_rx matrix of random vectors chosen

        % Sample covariance matrix
        S = Y'*Y;

        % Calculate the Bayes factor eq. (24)
        Bayes_factor = BayesFactor(S,F_cell,D_cell,delta_vec,h_tilde,alf);

        % detection decision
        Prob_det(lp_real, lp_snr) = Bayes_factor > Th_val;

    end

end

% -- Pd by averaging over realizations --
avg_prob_det = mean(Prob_det)';


% =================================================================================================
function [Th_val, Out_Pfa, BF] = DtnTh_RealData(p_rx, delta_vec, Omega_0, F_cell, D_cell, Pfa, K, h, num_ensembles, lp_class)
% [Th_val, Out_Pfa, BF] = DtnTh_RealData(p_rx, delta_vec, Omega_0, F_cell, D_cell, Pfa, K, h, num_ensembles, lp_class)
%
% The code provides the threshold value for the U_{BF} detector at a given
% prob of false alarm, where the data is real.
%
% Input
%   p_rx: scalar, number of Rx antennas.
%   delta_vec: (p_rx x 1), parameters for the Gamma distribution.
%   Omega_0: (p_rx x p_rx): Covariance matrix in the null hypothesis,
%   F_cell, D_cell: cell array, parameters of the multivariate normal prior of the
%     off-diagonal elements of H0, where the number of elements in F_cell or
%      D_cell is (p_rx - 1).
%   Pfa: scalar, desired prob of false alarm.
%   K: scalar, number of received samples.
%   h: parameter for the prior of H1.
%   num_ensembles: scalar, number of realizations.
%   lp_class: scalar, index for the objective prior, possible values are
%       1,..,5.
%
% Output
%   Th_val: scalar, threshold value for the given Pfa.
%   Out_Pfa: scalar, simulated prob. of false alarm to match with the desired Pfa.
%   BF (num_ensembles x 1):  Bayes factor values.
%
% Reference:
% M. H. Al-Ali and K.C. Ho, "Objective Bayesian approach for binary hypothesis
% testing of multivariate Gaussian observations," IEEE Trans. Inf. Theory,
% vol. 69, no. 2, pp. 1337-1354, Feb. 2023.
%
% M. H. Al-Ali and K. C. Ho, Dec. 12. 2022
%
%       Copyright (C) 2022
%       Compputer Engineering Department
%       University of Basrah
%       Basrah Province, Iraq
%       mohannd.khalaf@uobasrah.edu.iq
%
%       Copyright (C) 2022
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

% -- Objective prior parameters --
m = 1:p_rx;

if lp_class == 1            % Possible values for a_m: su_rx - k, su_rx + 1 - k, 2 - k, k, or 1
    a_m = p_rx - m;         % Jeffreys Prior
elseif lp_class == 2
    a_m = p_rx - m + 1;     % Independence jeffreys Prior
elseif lp_class == 3
    a_m = 2 - m;            % Geisser and Cornfield Prior
elseif lp_class == 4
    a_m = m;                % Right Haar Measure
elseif lp_class == 5
    a_m = ones(1, p_rx);    % Reference Prior
end

% -- Gamma distribution parameters --
h_tilde = h + K/2;

alf = (K - a_m + 1)/2;              % The exponent of the denominator in the Bayesian detector
alf_max = max(max(alf), h_tilde);   % Used to scale exponents of the detector

% -- Null model generation --
BF = zeros(num_ensembles, 1);
A = chol(Omega_0);

for lp_real = 1:num_ensembles
    Y = mvnrnd(zeros(1, p_rx), Omega_0, K);     % Generates an K-by-su_rx matrix of random vectors

    S = Y'*Y ;      % sample covariance matrix

    % Bayes factor calculation
    Bayes_factor = BayesFactor(S,F_cell,D_cell,delta_vec,h_tilde,alf);
    BF(lp_real, 1) = Bayes_factor;

end

% -- Determining threshold --
nbins = 1e4;    % step resolution to determine threshold
thresh_low = min(BF);
thresh_high  = max(BF);
thresh_steps = linspace(thresh_low,thresh_high,nbins);
pfa_count = zeros(1,nbins);

for k = 1:nbins
    thresh = thresh_steps(k);
    pfa_count(k) = sum(BF >= thresh);
end

pfa_sim = pfa_count/num_ensembles;
[err, indx] = min(abs(pfa_sim - Pfa));
Th_val = thresh_steps(indx);
Out_Pfa = pfa_sim(indx);

