INTRODUCTION
Let's supose we receive the following signal:
\begin{equation} y[n] = hx[n]+W[n] \end{equation}
where:
- $x[n] = Ae^{j2 \pi f_c t}$ is the transmitted signal
- $f_c$ is the carrier frequency
- $W[n]$ is AWGN complex noise, $W[n] \sim \mathcal{CN}(0,\sigma_w^2)$
- $h$ a complex normal random variable: $h \sim \mathcal{CN}(0,2\sigma_h^2)$ with real and imaginary parts i.i.d, which stands for Rayleigh Fading.
The aim here is to charactarize the distribution of the test statistic $T(Y) = \sum_{n=1}^N |y[n]|^2$. If we consider $h$ to be deterministic we obtain the following pdf:
\begin{equation} f_{T}^{fading}(x;|h|) = \frac{N}{\sigma_w^2} e^{-\frac{N}{\sigma_w^2} (x+|h|^2A^2)} \bigg(\frac{x}{|h|^2A^2}\bigg)^{\frac{1}{2}(N-1)} I_{N-1}\bigg(\frac{2N}{\sigma_w^2}|h|^2A\sqrt{x}\bigg) \end{equation}
Notice that $T$ follows a scaled non-central chi square distribution. This distribution has been proofed in MATLAB, you can run the following script to find out.
%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 = 2;
%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
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.
T_pdf_th = 2*N/sigma_w^2*pdf('Noncentral Chi-square',T_var*2*N/...
sigma_w^2,2*N,2*N*abs(h)^2*A^2/sigma_w^2);
plot(T_var,T_pdf,T_var,T_pdf_th);
title('Energy T(y) distribution');
legend('Experimental results','Theoretical results');
grid on;
PROBLEM
The problem occurs when I consider the distribution of $|h|$, since it is not a deterministic value: $|h| \sim Rayleigh(\sigma_h)$. Therefore, I suppose that the final expression should be something like:
\begin{equation} f_T^{Rayleigh}(x) = \int_{0}^{\infty} \frac{N}{\sigma_w^2} e^{-\frac{N}{\sigma_w^2} (x+|h|^2A^2)} \bigg(\frac{x}{|h|^2A^2}\bigg)^{\frac{1}{2}(N-1)} I_{N-1}\bigg(\frac{2N}{\sigma_w^2}|h|A\sqrt{x}\bigg) f_{|h|}(|h|)d|h| \tag{1} \end{equation}
where $f_{|h|}(x) = \frac{x}{\sigma_h^2} e^{-\frac{x^2}{2\sigma_h^2}}$ is the pdf of $|h|$.
Or maybe something like:
\begin{equation} f_T^{Rayleigh}(x) = \int_{0}^{\infty} \frac{N}{\sigma_w^2} e^{-\frac{N}{\sigma_w^2} (x+|h|^2A^2)} \bigg(\frac{x}{|h|^2A^2}\bigg)^{\frac{1}{2}(N-1)} I_{N-1}\bigg(\frac{2N}{\sigma_w^2}|h|A\sqrt{x}\bigg) f_{|h|^2}(|h|^2)d|h|^2 \tag{2} \end{equation}
where $f_{|h|^2}(x) = \frac{1}{2 \sigma_h^2} e^{-\frac{x}{2\sigma_h^2}}$ is the pdf of $|h|^2$.
However, I don't have any idea how can I compute these expressions with MATLAB, since I am not sure whether (1) or (2), or both, or none are correct.
Really appreciate your help, thanks.