0
$\begingroup$

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:

Autocorrelation Function of Jakes Model

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!

$\endgroup$
  • $\begingroup$ i would use confidence intervals as one kind of “validation”. $\endgroup$ – user28715 Jul 5 '19 at 15:59
  • $\begingroup$ 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. $\endgroup$ – JohnMarvin Jul 8 '19 at 0:11
0
$\begingroup$

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.

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.