Probability of True Positive of a random variable defined by an integral expression

Question

$\newcommand{\Prob}{\operatorname{Prob}}$Let's assume that we have a random variable with the following pdf:

\begin{equation} f_T(x) = \int_0^\infty f_T(x,g) \cdot f_{g}(g) \, dg = \int_0^\infty \frac{1}{\sqrt{2 \pi \big[(2gA^2\sigma_w^2+\sigma_w^4)/N\big]}} e^{-\frac{(x-(gA^2+\sigma_w^2))^2}{2 (2gA^2+\sigma_w^4)/N}} \cdot \frac{1}{2\sigma_h^2} e^{-\frac{g}{2\sigma_h^2}} \, dg \end{equation}

where:

• $A$, $N$, $\sigma_w$ and $\sigma_h$ are known deterministic values.

I want to obtain the mathematical expression for $\Prob(T>\lambda)$, where $\lambda$ is an arbitrary value. As far as I know I would compute it like following:

\begin{equation} \Prob(T>\lambda) = \int_\lambda^\infty \bigg[ \int_0^\infty \frac{1}{\sqrt{2 \pi \big[(2gA^2\sigma_w^2+\sigma_w^4)/N\big]}} e^{-\frac{(x-(gA^2+\sigma_w^2))^2}{2 (2gA^2\sigma_w^2+\sigma_w^4)/N}} \cdot \frac{1}{2\sigma_h^2} e^{-\frac{g}{2\sigma_h^2}} \, dg \bigg] \, dx = \int_0^\infty \bigg[ \int_{\lambda}^{\infty} \frac{1}{\sqrt{2 \pi \big[(2gA^2\sigma_w^2+\sigma_w^4)/N\big]}} e^{-\frac{(x-(gA^2+\sigma_w^2))^2}{2 (2gA^2\sigma_w^2+\sigma_w^4)/N}} \, dx \bigg] \cdot \frac{1}{2\sigma_h^2} e^{-\frac{g}{2\sigma_h^2}} \, dg \\= \int_{0}^{\infty} Q\Bigg(\frac{\lambda-(gA^2+\sigma_w^2)}{\sqrt{(2gA^2\sigma_w^2+\sigma_w^4)/N}}\Bigg) \cdot \frac{1}{2\sigma_h^2} e^{-\frac{g}{2\sigma_h^2}} \, dg \end{equation}

where $Q\big(\frac{x-\mu}{\sigma}\big) = \int_{x}^{\infty} f_x(x) dx$ where $f_x(x)$ is the pdf of a normal distributed random variable $X \sim \mathcal{N}(\mu,\sigma^2)$.

However this result seems not to be correct, since I have been trying to run it with MATLAB and it does not fit the experimental results. The pdf $f_T(x)$ does fit the experimental results.

Matlab code:

%RAYFAD_CHANNEL_SIM


S = 2e7;
N = 2e3;

% PU signal
fs = 150e3; % sampling frequency
Ts = 1/fs; % sampling period
T = (1:1:S)*Ts; % time vector
fc = 50e3; % carrier frequency
A = sqrt(7e-10); % PU signal amplitude
x = A*exp(2i*pi*fc*t); % PU signal

% Noise
sigma_w = sqrt(7e-10); % Noise variance
w = sigma_w/sqrt(2)*(randn(1,2e7)+1i*randn(1,2e7)); % Noise vector

% Fading coefficient
sigma_h = 2; % Rayleigh parameter

h = sigma_h*(randn(1,S/N)+1i*randn(1,S/N)); % h coefficient
h = repelem(h,N);

% Signal
y = h.*x + w; % h value changes every 2000 samples, being 2e7 the total 
              % number of samples, we have 1e4 different h values 

Y = reshape(abs(y),N,S/N); % Reshape the received signal into a 2e3x1e4 matrix


% This is the random variable !!!
T = mean(Y.^2,1); % Compute the avarage energy for each 2e3 subvector -> 
                      % We have 1e4 realizations of the energy vector T(y)

[T_pdf, T_var] = var2pdf(T,200); % Obtain the pdf of the energy vector T(y)
% var2pdf works the same way as MATLAB native function hist, but the values
% are normalized.


Mx = max(abs(h).^2); % Since infinit is not a valid number, we use the maximum achieved value.
cont = 0;
Pd_exp = zeros(size(T_var)); % This is probability we are computing (experimental)
Pd_th = Pd_exp; % This is probability we are computing (theoretical)


% Compute the probability, both experimental and analytical
for Tk=T_var
        cont = cont + 1;

            Pd_exp(cont) = sum(T > Tk)/length(T);

            % Theoretical (exact)


        ff = @(g) qfunc((Tk-g.*A.^2-sigma_w^2)/sqrt((2.*g.*A^2 ...
        .*sigma_w^2+sigma_w^4)/N))...
        .*1./(2.*sigma_h^2).*exp(-g./(2*sigma_h^2));
            Pd_th(cont) = integral(ff, 0, Mx);

end

% PLOT
plot(T_var,Pd_exp,'*'); hold on; plot(T_var,Pd_th,'g*'); title('Pd'); hold off;

Really appreciate your help

Answer 1

Ok, so the theoretical results were correct. The error was with the matlab code.

The expression of ff should be like:

ff = @(g) qfunc((Tk-(g*A^2+sigma_w^2))./(sqrt((2.*g*A^2*sigma_w^2+sigma_w^4)/N)))...
 .*1./(2.*sigma_h^2).*exp(-g./(2*sigma_h^2));

instead of

ff = @(g) qfunc((Tk-g.*A.^2-sigma_w^2)/sqrt((2.*g.*A^2.*sigma_w^2+sigma_w^4)/N))...
.*1./(2.*sigma_h^2).*exp(-g./(2*sigma_h^2));

I don't exactly get how the character . works here. Anyway, problem solved!