Validating Jakes model in Rayleigh Fading channels

Question

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!

Answer 1

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.