For the below mentioned code the plots for the Figures:1,3,5,6, is not getting ploted for entire range of values x-axis. and I cannot see the Fig 7 at all.
Kindly help and guide.
Thank You.
clear; clc; close all;
%% --- Common Setup ---
set(0, 'defaultAxesFontName', 'Times New Roman');
set(0, 'defaultTextInterpreter', 'tex');
set(0, 'defaultAxesFontSize', 11);
%% --- Simulation Parameters ---
c = 3e8; % Speed of light [m/s]
fc = 5.9e9; % Carrier frequency [Hz]
lambda = c / fc; % Wavelength [m]
N_realizations = 100; % Monte Carlo realizations
N_paths = 20; % Number of NLOS subpaths per cluster
L = 5; % Number of clusters
sigma_CSI = 0.1; % Std. dev. of CSI phase error [rad]
G_RIS = 5; % RIS amplification gain
K = 10; % Rice factor
sigma_a = 0.1; % Acceleration variance [m/s^2]
sigma_tau_base = [40e-9, 60e-9, 80e-9, 100e-9, 120e-9]; % Base delay spread per scenario [s]
sigma_xi = 4; % Shadowing std. dev. [dB]
v_max = 200 / 3.6; % Max speed [m/s]
T_values = linspace(0.5e-3, 5e-3, 20); % Feedback latency [ms]
M_T = 4; % Number of transmit antennas
V_values = round(linspace(8, 80, 20)); % Number of RIS elements, integer values
P_max = 10^(43/10) / 1000; % Max BS power [W]
P_v = 0.01; % RIS element power [W]
P_circuit = 5; % Circuit power [W]
chi_max = 2; % Max RIS amplification
% Mobility scenarios
scenarios = {'Dynamic Acc.', 'Steady Acc.', 'Uniform Vel.', 'Mobile-to-Static', 'Static-to-Static'};
v0 = [10, 10, 10, 10, 0]; % Initial velocities [m/s]
a = [6, 3, 0, 0, 0]; % Accelerations [m/s^2]
kappa = [60, 50, 40, 30, 20]; % TCF decay rates [1/s]
times = linspace(0, 4, 20); % 20 points for full coverage
kappa_scaling = 0.2; % Scaling for Dynamic Acc.
chi_l = 0.8 * ones(1, L); % Cluster gains
tau_0 = [150e-9, 180e-9, 200e-9, 220e-9, 250e-9]; % PDP decay constants [s]
% Geometry parameters
d_tx_rx_base = [40, 60, 80, 100, 120]; % Base Tx-Rx distance per scenario [m]
d_tx_ris = 100; % Distance Tx-RIS [m]
d_ris_rx = 100; % Distance RIS-Rx [m]
theta_AoA = rand(N_paths, L) * pi; % Random AoA per cluster [rad]
theta_AoD = rand(N_paths, L) * pi; % Random AoD per cluster [rad]
phi_AoA = rand(N_paths, L) * 2 * pi; % Random azimuth AoA [rad]
t = linspace(0, 0.04, 100); % Time difference for TCF [s]
delta_f = linspace(0, 1e6, 100); % Frequency separation for FCF [Hz]
tau = linspace(0, 700e-9, 50); % Delay for PDP [s]
%% --- Monte Carlo Simulation ---
sigma_D_scenarios = nan(length(scenarios), length(times));
Tc_scenarios = nan(length(scenarios), length(times));
Bc_scenarios = nan(length(scenarios), length(times));
EE_scenarios = nan(length(scenarios), length(times));
PL_LOS_scenarios = nan(length(scenarios), length(times));
PL_NLOS_scenarios = nan(length(scenarios), length(times));
G_RIS_scenarios = nan(length(scenarios), length(times));
v_x = nan(length(scenarios), length(times)); % Velocity [m/s]
sigma_tau_x = nan(length(scenarios), length(times)); % Delay spread [ns]
d_tx_rx_x = nan(length(scenarios), length(times)); % Tx-Rx distance [m]
d_total_x = nan(length(scenarios), length(times)); % Total cluster path distance [m]
for s = 1:length(scenarios)
for ti = 1:length(times)
sigma_D_temp = 0;
Tc_temp = 0;
Bc_temp = 0;
EE_temp = 0;
PL_LOS_temp = 0;
PL_NLOS_temp = 0;
G_RIS_temp = 0;
d_total_temp = 0;
v_x(s, ti) = v0(s) + a(s) * times(ti) + s * times(ti) * 0.1; % Increased offset
sigma_tau_x(s, ti) = sigma_tau_base(s) * (1 + 0.6 * times(ti)); % Increased scaling
d_tx_rx_x(s, ti) = d_tx_rx_base(s) * (1 + 0.6 * times(ti)); % Increased scaling
valid_realizations = 0;
for r = 1:N_realizations
[sigma_D_r, Tc_r, Bc_r, EE_r, PL_LOS_r, PL_NLOS_r, G_RIS_r, d_total_r] = compute_metrics(t, delta_f, tau, v0(s), a(s), times(ti), kappa(s), tau_0(s), N_paths, L, d_tx_rx_x(s, ti), d_tx_ris, d_ris_rx, ...
theta_AoA, theta_AoD, phi_AoA, fc, c, G_RIS, sigma_CSI, sigma_a, lambda, s == 1, kappa_scaling, K, chi_l, sigma_tau_x(s, ti), sigma_xi, M_T, V_values(ti), P_max, P_v, P_circuit, chi_max, v_max, T_values(ti));
if all(isfinite([sigma_D_r, Tc_r, Bc_r, EE_r, PL_LOS_r, PL_NLOS_r, G_RIS_r, d_total_r]))
sigma_D_temp = sigma_D_temp + sigma_D_r;
Tc_temp = Tc_temp + Tc_r;
Bc_temp = Bc_temp + Bc_r;
EE_temp = EE_temp + EE_r;
PL_LOS_temp = PL_LOS_temp + PL_LOS_r;
PL_NLOS_temp = PL_NLOS_temp + PL_NLOS_r;
G_RIS_temp = G_RIS_temp + G_RIS_r;
d_total_temp = d_total_temp + d_total_r;
valid_realizations = valid_realizations + 1;
end
end
if valid_realizations > 0
sigma_D_scenarios(s, ti) = sigma_D_temp / valid_realizations + s * 20; % Offset
Tc_scenarios(s, ti) = Tc_temp / valid_realizations + s * 0.1e-3; % Offset in ms
Bc_scenarios(s, ti) = Bc_temp / valid_realizations + s * 1e3; % Offset in Hz
EE_scenarios(s, ti) = EE_temp / valid_realizations * (1 + 0.2 * s); % Scaling
PL_LOS_scenarios(s, ti) = PL_LOS_temp / valid_realizations + s * 3; % Offset in dB
PL_NLOS_scenarios(s, ti) = PL_NLOS_temp / valid_realizations + s * 3; % Offset in dB
G_RIS_scenarios(s, ti) = G_RIS_temp / valid_realizations + s * 0.8; % Offset in dB
d_total_x(s, ti) = d_total_temp / valid_realizations + s * times(ti) * 3; % Increased offset
end
end
end
%% --- Figure: Doppler Spread vs. Velocity ---
figure('Position', [100 100 900 500]);
colors = {'r-o', 'b--s', 'g:^', 'm-.d', 'k-*'};
for s = 1:length(scenarios)
plot(v_x(s, :), sigma_D_scenarios(s, :), colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Velocity [m/s]');
ylabel('Doppler Spread, \sigma_D [Hz]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('Doppler Spread vs. Velocity for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'XLim', [0 35]);
%% --- Figure: Coherence Time vs. Time ---
figure('Position', [1050 100 900 500]);
for s = 1:length(scenarios)
plot(times, Tc_scenarios(s, :) * 1000, colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Time, t [s]');
ylabel('Coherence Time, T_c [ms]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('Coherence Time vs. Time for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'XLim', [0 4]);
%% --- Figure: Coherence Bandwidth vs. Delay Spread ---
figure('Position', [2000 100 900 500]);
for s = 1:length(scenarios)
plot(sigma_tau_x(s, :) * 1e9, Bc_scenarios(s, :) / 1e3, colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Delay Spread, \sigma_\tau [ns]');
ylabel('Coherence Bandwidth, B_c [kHz]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('Coherence Bandwidth vs. Delay Spread for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'XLim', [0 300]);
%% --- Figure: Energy Efficiency vs. Feedback Latency ---
figure('Position', [100 650 900 500]);
for s = 1:length(scenarios)
plot(T_values * 1e3, EE_scenarios(s, :), colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Feedback Latency, T [ms]');
ylabel('Energy Efficiency, EE [Mbits/Joule]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('Energy Efficiency vs. Feedback Latency for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'YScale', 'log', 'XLim', [0.5 5]);
%% --- Figure: Path Loss (LOS) vs. Tx-Rx Distance ---
figure('Position', [1050 650 900 500]);
for s = 1:length(scenarios)
plot(d_tx_rx_x(s, :), PL_LOS_scenarios(s, :), colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Tx-Rx Distance, d_{tx-rx} [m]');
ylabel('LOS Path Loss [dB]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('LOS Path Loss vs. Tx-Rx Distance for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'XLim', [0 300]);
%% --- Figure: Path Loss (NLOS) vs. Total Cluster Path Distance ---
figure('Position', [2000 650 900 500]);
for s = 1:length(scenarios)
plot(d_total_x(s, :), PL_NLOS_scenarios(s, :), colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Total Cluster Path Distance, d_{total} [m]');
ylabel('NLOS Path Loss [dB]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('NLOS Path Loss vs. Total Cluster Path Distance for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'XLim', [0 500]);
%% --- Figure: RIS Gain vs. Number of RIS Elements ---
figure('Position', [100 1200 900 500]);
for s = 1:length(scenarios)
plot(V_values, G_RIS_scenarios(s, :), colors{s}, 'LineWidth', 2, 'MarkerSize', 8); hold on;
end
grid on;
xlabel('Number of RIS Elements, V');
ylabel('RIS Gain [dB]');
legend(scenarios, 'Location', 'northwest', 'FontSize', 10);
title('RIS Gain vs. Number of RIS Elements for Different Mobility Scenarios');
set(gca, 'Box', 'on', 'FontSize', 11, 'XLim', [8 80]);
%% --- Function Definition ---
function [sigma_D, Tc, Bc, EE, PL_LOS, PL_NLOS, G_RIS, d_total] = compute_metrics(t, delta_f, tau, v, a, t_eval, kappa_s, tau_0_s, N_paths, L, d_tx_rx, d_tx_ris, d_ris_rx, theta_AoA, theta_AoD, phi_AoA, fc, c, G_RIS, sigma_CSI, sigma_a, lambda, is_dynamic_acc, kappa_scaling, K, chi_l, sigma_tau, sigma_xi, M_T, V, P_max, P_v, P_circuit, chi_max, v_max, T)
% Initialize
h = zeros(N_paths * L + 1, length(tau)); % CIR
f_D = zeros(N_paths * L + 1, 1); % Doppler shifts
h_eff = zeros(M_T, 1); % Effective channel
h_direct = zeros(M_T, 1); % Direct channel
r_MT = [0, 0, 0]'; % MT position
r_MR = [max(d_tx_rx, 1e-3), 0, 0]'; % MR position, ensure non-zero
r_C = rand(3, L) * 150 + [d_tx_ris; 0; 0]; % Cluster positions, increased range
v_MT = [v + a * t_eval, 0, 0]'; % MT velocity
v_MR = [v + a * t_eval, 0, 0]'; % MR velocity (simplified)
v_C = randn(3, L) * sigma_a; % Cluster velocities
lambda_CSI = besselj(0, 2 * pi * (v_max * fc / c) * T); % CSI correlation
% LOS component
tau_LOS = d_tx_rx / c;
u_LOS = (r_MR - r_MT) / norm(r_MR - r_MT);
f_D_LOS = dot(v_MR - v_MT, u_LOS) / lambda;
h_LOS = sqrt(K / (K + 1) / (4 * pi * d_tx_rx / c)) * exp(1j * 2 * pi * f_D_LOS * t_eval) * exp(1j * normrnd(0, sigma_CSI));
h(1, :) = h_LOS * exp(-1j * 2 * pi * fc * (tau - tau_LOS));
f_D(1) = f_D_LOS;
h_direct = h_LOS * ones(M_T, 1); % Simplified MIMO
% NLOS components (RIS-reflected)
idx = 2;
H_r_b = sqrt(1 / (4 * pi * d_tx_ris / c)) * exp(1j * randn(M_T, V)); % BS-RIS channel
h_r = sqrt(1 / (4 * pi * d_ris_rx / c)) * exp(1j * randn(V, 1)); % RIS-MR channel
Theta = diag(chi_max * exp(1j * 2 * pi * rand(V, 1))); % RIS matrix
d_total = 0;
for ell = 1:L
tau_NLOS = (norm(r_C(:, ell) - r_MT) + norm(r_MR - r_C(:, ell))) / c;
d_total_ell = norm(r_C(:, ell) - r_MT) + norm(r_MR - r_C(:, ell));
if d_total_ell < 1e-3
d_total_ell = 1e-3; % Prevent zero distance
end
d_total = d_total + d_total_ell / L;
for n = 1:N_paths
u_T = [cos(phi_AoA(n, ell)) * cos(theta_AoA(n, ell)), sin(phi_AoA(n, ell)) * cos(theta_AoA(n, ell)), sin(theta_AoA(n, ell))]';
u_R = [cos(phi_AoA(n, ell)) * cos(theta_AoA(n, ell)), sin(phi_AoA(n, ell)) * cos(theta_AoA(n, ell)), sin(theta_AoA(n, ell))]';
f_D_NLOS = dot(v_MT - v_C(:, ell), u_T) / lambda + dot(v_MR - v_C(:, ell), u_R) / lambda;
h_NLOS = sqrt(chi_l(ell)^2 / (K + 1) * G_RIS / (4 * pi * d_total_ell / c)) * exp(1j * 2 * pi * f_D_NLOS * t_eval) * exp(1j * normrnd(0, sigma_CSI));
h(idx, :) = h_NLOS * exp(-1j * 2 * pi * fc * (tau - tau_NLOS));
f_D(idx) = f_D_NLOS;
idx = idx + 1;
end
end
h_eff = h_direct + H_r_b * Theta * h_r; % Effective channel
% Doppler Spread
sigma_D = sqrt((K / (K + 1)) * (f_D(1) - mean(f_D))^2 + sum((chi_l.^2 / (K + 1)) .* mean(reshape((f_D(2:end) - mean(f_D)).^2, N_paths, L), 1)));
% Coherence Time
if is_dynamic_acc
kappa_s = kappa_s * (1 + t_eval * kappa_scaling);
end
Tc = log(2) / kappa_s;
% Coherence Bandwidth
FCF = (K / (K + 1)) * exp(1j * 2 * pi * delta_f * tau_LOS);
for ell = 1:L
tau_NLOS = (norm(r_C(:, ell) - r_MT) + norm(r_MR - r_C(:, ell))) / c;
FCF = FCF + (chi_l(ell)^2 / (K + 1)) * exp(1j * 2 * pi * delta_f * tau_NLOS);
end
FCF = abs(FCF) / max(abs(FCF));
Bc = delta_f(find(FCF <= 0.5, 1, 'first'));
if isempty(Bc)
Bc = delta_f(end);
end
% Energy Efficiency (simplified)
w_c = sqrt(P_max / (M_T + 1)) * ones(M_T, 1); % Common stream precoding
w_k = sqrt(P_max / (M_T + 1)) * ones(M_T, 1); % Private stream precoding
sigma_r = 1e-3; % RIS noise variance
sigma_k = 1e-3; % AWGN variance
gamma_c = (abs(w_c' * (lambda_CSI * h_eff + sqrt(1 - lambda_CSI^2) * randn(M_T, 1)))^2) / ...
(abs(w_k' * (lambda_CSI * h_eff + sqrt(1 - lambda_CSI^2) * randn(M_T, 1)))^2 + sigma_r * sum(abs(diag(Theta)).^2) + sigma_k);
R_c = log2(1 + gamma_c);
R_k = log2(1 + gamma_c); % Simplified for one CV
P_total = norm(w_c)^2 + norm(w_k)^2 + sum(abs(diag(Theta)).^2) * P_v + P_circuit;
EE = (R_c + R_k) / P_total;
% Path Loss
PL_LOS = -20 * log10(sqrt(K / (K + 1)) * lambda / (4 * pi * d_tx_rx));
PL_NLOS = 0;
for ell = 1:L
d_total_ell = norm(r_C(:, ell) - r_MT) + norm(r_MR - r_C(:, ell));
if d_total_ell < 1e-3
d_total_ell = 1e-3; % Prevent zero distance
end
PL_NLOS = PL_NLOS + (chi_l(ell)^2 / L) * (-20 * log10(sqrt(1 / (K + 1)) * lambda / (4 * pi * d_total_ell)) + normrnd(0, sigma_xi));
end
% RIS Gain
G_RIS = 10 * log10(mean(abs(h_eff).^2) / mean(abs(h_direct).^2));
end
