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!