MATLAB: Validating Jakes model in Rayleigh Fading channels

MATLABrayleigh fadingwireless communications

Hi,
I am trying to validate the Jakes model by comparing the autocorrelation function given by , where is the maximum Doppler shift, is the zero-th order Bessel function of the first kind, and 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 = 1e3;  % 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. It is
                                    % remove the number of sampls
                                    % in 2 OFDM symbols, one
                                    % transmitted by Alice and
                                    % other by Bob, so the total
                                    % transmission time results in
                                    % Delta t.
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')
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!

Best 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 ρ 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 ρ value is to accumulate the channel values in a vector and, then, compute the correlation value after the trials loops.
Related Question