# Validating Jakes model in Rayleigh Fading channels

I am trying to validate the Jakes model by comparing the autocorrelation function given by $$\rho(\Delta t ) = J_{0}(2 \pi f_{d} \Delta t)$$, where $$f_{d}$$ is the maximum Doppler shift, $$J_{0}$$ is the zero-th order Bessel function of the first kind, and $$\Delta t$$ is the time interval between channel samples. However, I am obtaining the following figure:

As it is possible to see, the simulation values does not exactly match the simulated values. I would like to know if this is expected or if I am doind anything wrong.

Also, this is the code I am using:

clear variables; %close all;

%% Simulation parameters
Ntrials = 1e2;  % Number of trials
M = 4;          % Modulation order (4 for QPSK)

%% Modulators
% Quadrature modulator and demodulator:
hMod = comm.RectangularQAMModulator( ...
'ModulationOrder', M, ...
'BitInput', true, ...
'NormalizationMethod', 'Average power');
hDemod = comm.RectangularQAMDemodulator( ...
'ModulationOrder', M, ...
'BitOutput', true, ...
'NormalizationMethod', 'Average power');
Nbpsym = log2(M); % Number of bits per symbol

%% Channel parameters
SNR = 100;      % SNR values for simulation
Nsnr = length(SNR);
symRate = 1e6; % Original value: 20e6
trms = 500e-9; % Original value: 25e-9
[pdp, tau] = IEEE802_11_exppdp(trms, 1/symRate); % Channel Power Profile
pdB = 10*log10(pdp);    % Channel Power Profile (dB)
fd = 6;                 % Maximum Doppler spread in Hz
L = length(pdB);        % Channel length in time

% Creates Rayleigh Channel Object
rayChan = comm.RayleighChannel( ...
'SampleRate', symRate, ...
'PathDelays', tau, ...
'AveragePathGains', pdB, ...
'MaximumDopplerShift', fd, ...
'PathGainsOutputPort', true);

%% Channel probing parameters
Tc = sqrt(9/(16*pi))/fd; % Channel coherence time from Clarke's model (seconds)
Tsym = 1/symRate; % Symbol time
Deltat = linspace(Tsym, 4*Tc, 10);  % Interval between subsequent
% measurements. DeltaT must be
% at minimum Tsym, since
% there wil be at least one
% OFDM symbol transmitted by
% Alice and other transmitted
% by Bob
NDeltat = length(Deltat);
Nsym = ceil(symRate*Deltat);        % Number of symbols necessary for
% each Delta t value.
Nb = Nbpsym*Nsym; % Number of bits to meet the required \Delta t

%% Initializations
rhohhat = zeros(Nsnr, NDeltat);
rhoHhat = zeros(Nsnr, NDeltat);
rhoHhatz = zeros(Nsnr, NDeltat);
rhoh = zeros(Nsnr, NDeltat);

% Wait bar:
cont = 0;
barStr = 'Simulation';
progbar = waitbar(0, 'Initializing waitbar...', ...
'Name', barStr);

%% Simulation
% parfor iDeltaT = 1:NDeltaT
for iDeltat = 1:NDeltat

fprintf(['Simulations for ' num2str(Deltat(iDeltat)*1000) ' ms.\n']);

for iSNR = 1:Nsnr

for trial = 1:Ntrials

%% Alice transmitter
bitsA = randsrc(Nb(iDeltat), 1, [0 1]);     % Bit stream in Alice
xA = step(hMod, bitsA);            % QAM Modulation

%% Channel
[yB, h] = rayChan(xA);

%% Channel correlation
if(trial ~= 1)
rhoh(iSNR, iDeltat) = rhoh(iSNR, iDeltat) + ...
corr(real(h(1, :).'), real(hant), 'type', 'Pearson');
end
hant = h(1, :).';

%% Progress:
cont = cont + 1;
% Update bar:
prog = cont/(Ntrials*Nsnr*NDeltat);
perc = 100*prog;
waitbar(prog, progbar, sprintf('%.2f%% Concluded', perc));

end

end

end

rhoh = rhoh/Ntrials;

close(progbar)

% Theoretical values:
Deltat_theo = 0:0.01:Deltat(end);
z = 2*pi*fd*Deltat_theo;
r = besselj(0, z);

%% Figures

figure
plot(Deltat, rhoh(end, :), 'o')
hold on;
plot(Deltat_theo, r);
legend('Simulation', 'Theoretical')
xlabel('\Deltat (s)'); ylabel('\rho(\Delta t)')


And also:

function [pdp, tau] = IEEE802_11_exppdp(sigma_t, Ts)
% Generate IEEE 802.11 power delay profile
% Input:
%   - sigma_t - RMS delay spread
%   - Ts - Sampling time
% Output:
%   - pdp - Power delay profile
% Example:
%   Ts = 50e-9;
%   trms = 25e-9
%

lmax = ceil(10*sigma_t/Ts); % Eq.(2.6)
sigma0 = (1-exp(-Ts/sigma_t))/(1-exp(-(lmax+1)*Ts/sigma_t)); % Eq.(2.9)
tau = (0:lmax)*Ts;
pdp = sigma0*exp(-tau/sigma_t); % Eq.(2.8)


I appreciate any insight or help about this.

Thanks!

• i would use confidence intervals as one kind of “validation”. – user28715 Jul 5 '19 at 15:59
• Hi @StanleyPawlukiewicz Do you think that would be enough? I thought about that, however, there seems to be some kind of stretching in the data I am computing, something that seems to go beyond only some error on the estimation of the correlation coefficient. Another thing is that I am using the actual channel taps to calculate the correlation coefficient, not some estimate of them, which makes me more confused. – JohnMarvin Jul 8 '19 at 0:11

## 1 Answer

Just to let anyone who might come across this post and is having the same issue.

I figured out the problem.

One cannot average the $$\rho$$ as the output of the corr function, since this estimate is not an ergodic random variable.

I checked this by using the following code:

% Test on averaging the result of corr function
clear variables; close all;

% Parameters on correlation:
rho = 0.75;
L = 6;

% XA = randn(N, 1);
% XB = rho(iRho)*XA + sqrt(1-rho(iRho)^2)*randn(N, 1);

%% Monte Carlo simulation:
Ntrials = 1e3;
rhoexp = 0;

for k = 1:Ntrials

% Generate data:
h1 = randn(L, 1);
h2 = rho*h1 + sqrt(1-rho^2)*randn(L, 1);

%% Channel correlation
rhoexp = rhoexp + corr(h1, h2, 'type', 'Pearson');

end

rhoexp = rhoexp/Ntrials;

fprintf('Theoretical: %f \nExperimental: %f\n', rho, rhoexp)

%% Test with no averaging:
h1 = zeros(L, Ntrials);
h2 = zeros(L, Ntrials);

for k = 1:Ntrials

% Generate data:
h1(:, k) = randn(L, 1);
h2(:, k) = rho*h1(:, k) + sqrt(1-rho^2)*randn(L, 1);

end
h1 = h1(:);
h2 = h2(:);
%% Channel correlation
rhoexp = corr(h1, h2, 'type', 'Pearson');

fprintf('Theoretical: %f \nExperimental: %f\n', rho, rhoexp)


Therefore, the right way of computing the $$\rho$$ value is to accumulate the channel values in a vector and, then, compute the correlation value after the trials loops.