F10 generation

explain in depth params.scenario='O1_28'; % DeepMIMO Dataset scenario

params.active_BS=3; D_Lambda = 0.5; BW = 100e6; Ut_row = 850; Ut_element = 90; Ur_rows
= [1000 1200]; Validation_Size = 6200; K = 512; miniBatchSize = 500; Mx = 1; M = Mx * My
* Mz;
ChatGPT said:

Great — let's break this segment down line by line, focusing on the context of the simulation in
the paper “Enabling Large Intelligent Surfaces with Compressive Sensing and Deep Learning”.
These lines appear to be inside Main_fn.m and configure the DeepMIMO dataset parameters,
antenna configuration, and model input dimensions.

params.scenario = 'O1_28';

Meaning:

 This selects the channel propagation scenario for the DeepMIMO dataset.
 'O1_28' stands for:
o ‘O1’: A layout (Outdoor 1) from the DeepMIMO project.
o ‘28’: Refers to 28 GHz carrier frequency (mmWave).

Context:

 In DeepMIMO, each scenario represents a ray-traced 3D urban environment using

Remcom Wireless InSite.
 'O1_28' includes urban buildings and users placed in grid patterns, with ray-traced
channels precomputed at 28 GHz.

params.active_BS = 3;

Meaning:

 Selects Base Station 3 from the DeepMIMO O1 scenario.

 The DeepMIMO dataset includes multiple base stations — BS3 is one of them.

Context:

 The code will extract the channel data only from BS3 to a large set of users.
 Important: Each BS has a different position and orientation, which affects the channel
sparsity and user-BS visibility.
D_Lambda = 0.5;

Meaning:

 This sets the element spacing on the LIS in units of wavelength (λ).
 D_Lambda = 0.5 means: the elements are spaced at half-wavelength distance, which is
standard in antenna arrays.

Context:

 Ensures no grating lobes in uniform linear or rectangular arrays.

 In practice, if λ = 1 cm at 28 GHz, then each element is 0.5 cm apart.

BW = 100e6;

Meaning:

 Sets the system bandwidth to 100 MHz, a typical bandwidth for 5G mmWave channels.

Context:

 This is used for:

o Computing data rate using Rate = BW * log2(1 + SNR)
o Normalizing power across subcarriers
o Applying appropriate scaling in FFT/IFFT during signal processing (if needed)

Ut_row = 850; Ut_element = 90;

Meaning:

 These define which users in the DeepMIMO grid are used for training/testing.

In DeepMIMO:

 Users are laid out in a grid of rows and columns.

 Each user has a position and precomputed channels from each BS.
 Ut_row = 850 and Ut_element = 90 selects a specific user:
o Row 850, user index 90
o This may be the receiver for which the channel is used for beamforming or rate
prediction

Ur_rows = [1000 1200];

Meaning:
  Range of user rows used for training/validation/test data.
 i.e., this code will collect user data from rows 1000 to 1200, giving ~201 rows × some
users/row = many user samples.

Context:

 These users form the training dataset for supervised learning (input: partial CSI; output:
best beam index or reflection pattern).

Validation_Size = 6200;

Meaning:

 Number of users taken from Ur_rows used for validation or testing dataset.

Why?

 After training on a larger dataset, we evaluate DL model generalization on

Validation_Size unseen samples.
 Used to compute the rate metrics in Figure 10.

K = 512;

Meaning:

 Number of OFDM subcarriers in the system (used in DeepMIMO).

 The DeepMIMO channel matrices are frequency-selective with 512 subcarriers.
 Each subcarrier has a channel vector per user and BS antenna.

Note:

 The deep learning model will later only use a subset of K (like K_DL = 64) for
compressive input.

miniBatchSize = 500;

Meaning:

 Batch size for training the deep learning model.

 Standard DL hyperparameter that controls:
o How many training samples are processed before one gradient update.
o Affects convergence and memory usage.
Mx = 1;

Meaning:

 Mx is the number of elements along the x-dimension of the LIS.

 Since it’s a 2D surface (Y-Z plane), Mx = 1 implies the LIS is a flat 2D surface (not 3D).

M = Mx * My * Mz;

Meaning:

 Total number of LIS elements (both active and passive).

 My, Mz are previously defined (e.g., 32 and 64)
 Mx = 1 → the array is 1 × 32 × 64 = 2048 elements.

💡 Summary

Variable Meaning Use

scenario = Urban 28 GHz mmWave Used by DeepMIMO to load ray-traced
'O1_28' scenario channels
active_BS = 3 Select BS#3 Only channels from this BS are used
D_Lambda = 0.5 Element spacing = 0.5λ Avoids grating lobes
BW = 100e6 100 MHz bandwidth Used for rate calculation
Ut_row, Ut_element Target user index For evaluation
Ur_rows Range of user rows for dataset Source of training/testing users
Validation_Size 6200 users For evaluating DL performance
K = 512 # of OFDM subcarriers Full frequency-domain channel
miniBatchSize DL training batch size Optimizer parameter
Mx, My, Mz LIS dimensions Total # of elements: M = Mx * My * Mz

Rate_DLt=zeros(numel(My_ar),numel(Training_Size));

Rate_OPTt=zeros(numel(My_ar),numel(Training_Size));

for rr = 1:1:numel(My_ar)

save Fig10_data.mat L My_ar Mz_ar M_bar Training_Size K_DL Rate_DLt Rate_OPTt

[Rate_DL,Rate_OPT]=Main_fn(L,My_ar(rr),Mz_ar(rr),M_bar,K_DL,Pt,kbeams,Training_Size);

Rate_DLt(rr,:)=Rate_DL; Rate_OPTt(rr,:)=Rate_OPT;
End

save Fig10_data.mat L My_ar Mz_ar M_bar Training_Size K_DL Rate_DLt Rate_OPTt

This block of code is responsible for generating data for Figure 10, which compares the
achievable rates of:

 DL-based beamforming
 Genie-aided (optimal) beamforming

across various LIS sizes and training dataset sizes.

Rate_DLt = zeros(numel(My_ar), numel(Training_Size));

Purpose: Pre-allocate a 2D matrix for storing Deep Learning-based achievable rates.

Details:

My_ar is an array of possible values for the number of LIS elements along the y-axis.

Training_Size is a vector of different dataset sizes used to train the deep learning model.

Rate_DLt(i,j) will contain the achievable rate using DL for the:

i-th LIS configuration (from My_ar)

j-th training dataset size

Rate_OPTt = zeros(numel(My_ar), numel(Training_Size));

Purpose:

Pre-allocate a matrix for optimal (Genie-aided) beamforming rates for comparison.

Note:

The "Genie-aided" method uses perfect channel knowledge to select the best reflection beam. It’s an
upper bound for DL performance.

for rr = 1:1:numel(My_ar)
Purpose: Start a loop over the LIS sizes defined in My_ar.

Loop Behavior:

For each value in My_ar, the system simulates the training and evaluation process using that LIS
configuration.

save Fig10_data.mat L My_ar Mz_ar M_bar Training_Size K_DL Rate_DLt Rate_OPTt

Purpose: Save the current workspace variables to a .mat file (in case of crash or checkpointing).

Useful for long simulations, to avoid data loss.

File Includes:

Channel and system params: L, My_ar, Mz_ar, M_bar, Training_Size, K_DL

Current results: Rate_DLt, Rate_OPTt

[Rate_DL, Rate_OPT] = Main_fn(L, My_ar(rr), Mz_ar(rr), M_bar, K_DL, Pt, kbeams, Training_Size);

Purpose: Call the main simulation function Main_fn.m for a specific LIS configuration.

Arguments Explained:

Param Meaning

L Number of channel paths (usually 1)

My_ar(rr) Y-axis elements for current LIS size

Mz_ar(rr) Z-axis elements for current LIS size

M_bar Number of active elements used (others are passive)

K_DL Number of subcarriers used as input to DL model

Pt Transmit power in dB

kbeams Number of best beams to consider

Training_Size Vector of dataset sizes

What It Returns:

Rate_DL: Achievable rate using the Deep Learning-based beamformer

Rate_OPT: Achievable rate using perfect (Genie-aided) beam selection

These are arrays of the same size as Training_Size.

Rate_DLt(rr,:) = Rate_DL;

Rate_OPTt(rr,:) = Rate_OPT;

Purpose: Store the outputs for this LIS configuration (My_ar(rr) and Mz_ar(rr)) into their respective
matrices.

Impact: These results will later be plotted as curves showing rate vs. training dataset size for different
LIS sizes.

end

Purpose: End of the loop over different LIS sizes.

save Fig10_data.mat L My_ar Mz_ar M_bar Training_Size K_DL Rate_DLt Rate_OPTt

Purpose:

Save final results into a .mat file for plotting Figure 10.

This file will be loaded by the next block to generate the plot.

What’s the Overall Flow?

Loop over two LIS configurations: e.g., My_ar = [32, 64]

For each:

Run DL training/evaluation over multiple training dataset sizes

Compare it with optimal beamforming

Store and save the achievable rate results

Dimensions:

Let’s assume:

My_ar = [32, 64]; % → numel = 2

Training_Size = [2, 10000, 14000, 18000, 22000, 26000, 30000]; % → numel = 7

Then:

size(Rate_DLt) = [2, 7];

Each row → a LIS size

Each column → a dataset size

Each cell → an achievable rate for that configuration

Summary Table:

Line Function

zeros(...) Allocate memory for results

for loopIterate over LIS configurations

Main_fn(...) Simulate training and test of DL vs. optimal

save ... Periodic checkpointing

Rate_DLt(rr,:) = ... Store results per LIS size

end Close loop

Final save Export data for plotting

Colour = 'brgmcky';

f10 = figure('Name', 'Figure10', 'units','pixels');

hold on; grid on; box on;

title(['Achievable Rate for different dataset sizes using only ' num2str(M_bar) ' active
elements'],'fontsize',12)

xlabel('Deep Learning Training Dataset Size (Thousands of Samples)','fontsize',14)

ylabel('Achievable Rate (bps/Hz)','fontsize',14)

set(gca,'FontSize',13)

if ishandle(f10)

set(0, 'CurrentFigure', f10)

hold on; grid on;

for rr=1:1:numel(My_ar)

plot((Training_Size*1e-3),Rate_OPTt(rr,:),[Colour(rr) '*--'],'markersize',8,'linewidth',2,
'DisplayName',['Genie-Aided Reflection Beamforming, M = ' num2str(My_ar(rr)) '*' num2str(Mz_ar(rr))])

plot((Training_Size*1e-3),Rate_DLt(rr,:),[Colour(rr) 's-'],'markersize',8,'linewidth',2, 'DisplayName',

['DL Reflection Beamforming, M = ' num2str(My_ar(rr)) '*' num2str(Mz_ar(rr))])

end

legend('Location','SouthEast')

legend show

end

drawnow

hold off
Let’s deeply analyze this MATLAB plotting code line by line, understanding both:

 Code functionality, and

 The mathematical/graphical meaning in the context of the paper “Enabling Large
Intelligent Surfaces with Compressive Sensing and Deep Learning”.

GOAL OF THIS BLOCK

To generate Figure 10, which compares:

 Genie-aided (optimal) reflection beamforming vs.

 Deep Learning-based reflection beamforming
...over different training dataset sizes and LIS configurations (different M = My × Mz
active elements).

🔢 DATA INPUTS USED

 Training_Size: array of DL training dataset sizes

 Rate_OPTt: achievable rates using optimal reflection beam
 Rate_DLt: achievable rates using DL-predicted reflection beam
 My_ar, Mz_ar: LIS sizes in Y and Z direction
 M_bar: number of active elements used
 Colour: character string for plotting color markers

🔍 LINE-BY-LINE EXPLANATION

Line 1
Colour = 'brgmcky';

🎨 Purpose:

 Define a sequence of 7 plot colors/markers:

'b' = blue, 'r' = red, 'g' = green, etc.
 Each entry corresponds to a different LIS configuration (My × Mz)

🧠 Used later to differentiate curves visually.

Line 2
f10 = figure('Name', 'Figure10', 'units','pixels');

📊 Purpose:

 Create a new figure window, stored in f10.

 'Name' sets the title in the window.
 'units','pixels' defines figure dimension units.

Line 3
hold on; grid on; box on;

📐 Purpose:

 hold on: Allows multiple plots to appear on the same figure.

 grid on: Displays background grid for clarity.
 box on: Draws a box around the plot area.

Line 4
title(['Achievable Rate for different dataset sizes using only '
num2str(M_bar) ' active elements'],'fontsize',12)

📝 Purpose:

 Adds a title to the figure.

 Describes what the plot shows:
o Achievable Rate vs. DL training size
o For M_bar active elements

📢 Output title:

Achievable Rate for different dataset sizes using only 16 active elements (for example)

Lines 5–6
xlabel('Deep Learning Training Dataset Size (Thousands of
Samples)','fontsize',14)
ylabel('Achievable Rate (bps/Hz)','fontsize',14)

Purpose:

 Label the x-axis and y-axis:

o X = training data size (in 1,000s of samples)
o Y = achieved spectral efficiency (rate in bps/Hz)

Line 7
set(gca,'FontSize',13)

🔠 Purpose:

 Set the font size of the current axes labels and ticks for better readability.

Line 8
if ishandle(f10)

🧠 Purpose:

 Check that f10 is still a valid figure handle.

 Ensures robustness: if the figure was accidentally closed, the rest won't run.

Line 9
set(0, 'CurrentFigure', f10)

Purpose:

 Set focus back to the figure f10 in case another figure opened.

Line 10
hold on; grid on;

✅ Redundant but safe: ensures grid is on and figure is in multi-plot mode.

Loop Block: Lines 11–14
matlab
CopyEdit
for rr = 1:1:numel(My_ar)

🔁 Purpose:

 Loop over each LIS configuration (each value of My_ar)

 For each, plot both optimal and DL-based performance

Line 12
matlab
CopyEdit
plot((Training_Size*1e-3), Rate_OPTt(rr,:), [Colour(rr) '*--'], ...
'markersize', 8, 'linewidth', 2, ...
'DisplayName', ['Genie-Aided Reflection Beamforming, M = '
num2str(My_ar(rr)) '*' num2str(Mz_ar(rr))])

🔹 Purpose:

 Plot optimal beamforming rate vs. training size

📌 Key Parts:

Part Description

Training_Size*1e-3 Convert size to thousands

Rate_OPTt(rr,:) Genie-aided rate for LIS config rr

[Colour(rr) '*--'] Marker = asterisk; Line = dashed

'DisplayName' Label used for legend

'M = Y×Z' Size of the LIS for this config

Line 13
plot((Training_Size*1e-3), Rate_DLt(rr,:), [Colour(rr) 's-'], ...
'markersize', 8, 'linewidth', 2, ...
 'DisplayName', ['DL Reflection Beamforming, M = '
num2str(My_ar(rr)) '*' num2str(Mz_ar(rr))])

🔹 Purpose:

 Plot Deep Learning-based rate for same LIS config

📌 Key Parts:

 's-':Square marker, solid line

 Otherwise same format as Line 12

🔍 Now we can compare DL vs. optimal on same plot for each LIS size.

Line 14–15
legend('Location','SouthEast')
legend show

Purpose:

 Add a legend to the plot showing which curve corresponds to which method and LIS size.
 'SouthEast' = bottom right of plot.

Line 16
end

✅ Ends the if ishandle(f10) block.

Line 17
drawnow

Purpose:

 Force MATLAB to immediately render the figure, rather than waiting.

 Useful for real-time visualization in long loops.
Line 18
hold off

🔚 Purpose:

 Releases the figure from holding more plots.

📊 Final Output: What Figure 10 Shows

A set of curves:

 X-axis: DL training size (in thousands)

 Y-axis: Achievable rate (bps/Hz)

For each LIS configuration:

 One dashed curve with asterisks: Genie-aided beamforming

 One solid curve with squares: Deep Learning beamforming

Helps visualize how DL approaches optimal with more training data, and how LIS size affects
performance.

Main function

Function Definition and Description Block

function
[Rate_DL,Rate_OPT]=Main_fn(L,My,Mz,M_bar,K_DL,Pt,kbeams,Training_Size)

 Inputs:
o L: Number of channel paths.
o My, Mz: LIS dimensions along y and z axes.
o M_bar: Number of active LIS elements (compressed measurement).
o K_DL: Number of OFDM subcarriers for DL input.
o Pt: Transmit power.
o kbeams: Number of top beams to consider.
o Training_Size: Array of training dataset sizes to evaluate.

params.scenario='O1_28';
params.active_BS=3;

 Load DeepMIMO Scenario O1_28 and specify BS index 3 as the active one.

D_Lambda = 0.5;
BW = 100e6;

 Antenna spacing is 0.5λ (standard to avoid grating lobes).

 Bandwidth: 100 MHz.

Ut_row = 850; Ut_element = 90; Ur_rows = [1000 1200];

 Specify fixed Tx user (Ut) and a range of Rx users (Ur) for validation.

Validation_Size = 6200; K = 512; miniBatchSize = 500;

 Dataset size for validation, number of subcarriers (OFDM), batch size.

Mx = 1;
M = Mx.*My.*Mz;

 LIS is assumed to be on y-z plane (Mx=1). Total LIS elements: M.

Output Initialization
Rate_DL = zeros(1,length(Training_Size));
Rate_OPT = Rate_DL;
LastValidationRMSE = Rate_DL;

 Preallocate storage for output rates and RMSE for different training sizes.

Noise and SNR Modeling

Gt=3; Gr=3; NF=5; Process_Gain=10;

 Tx/Rx antenna gains, noise figure, and processing gain for estimation.

noise_power_dB = -204 + 10*log10(BW/K) + NF - Process_Gain;

SNR = 10^(.1*(-noise_power_dB))*(10^(.1*(Gt+Gr+Pt)))^2;

 Noise power and corresponding SNR formula derived from Friis’ equation.

noise_power_bar = 10^(.1*(noise_power_dB))/(10^(.1*(Gt+Gr+Pt)));
  Used to compute standard deviation of additive Gaussian noise for channels.

User Pair Configuration

No_user_pairs = (Ur_rows(2)-Ur_rows(1))*181;
RandP_all = randperm(No_user_pairs).';

 Generate total user pairs and shuffle the indices randomly.

Display and Beamforming Codebook

disp([' Calculating for M = ' num2str(M)]);
Rand_M_bar_all = randperm(M);

 Randomly select a subset of M_bar elements from full LIS.

[BF_codebook]=sqrt(Mx*My*Mz)*UPA_codebook_generator(Mx,My,Mz,...);

 Generate UPA beamforming codebook scaled by √M to maintain power.

DeepMIMO Dataset Parameters

params.num_ant_x= Mx;
params.num_ant_y= My;
params.num_ant_z= Mz;
...
params.num_paths=L;

 Define LIS geometry, channel resolution (OFDM), and number of rays.

Dataset Generation: Transmitting User

params.active_user_first=Ut_row;
params.active_user_last=Ut_row;
DeepMIMO_dataset=DeepMIMO_generator(params);
Ht = single(DeepMIMO_dataset{1}.user{Ut_element}.channel);

 Load fixed TX user channel Ht (size: M x K).

Dataset Generation: Receiving Users
Validation_Ind = RandP_all(end-Validation_Size+1:end);
Ur_rows_step = 100;
Ur_rows_grid=Ur_rows(1):Ur_rows_step:Ur_rows(2);
Delta_H_max = single(0);

 Grid Ur_rows to load subsets of users. Will loop in batches of 100 rows.

for pp = 1:1:numel(Ur_rows_grid)-1

 Loop to find the max |Ht .* Hr| (used for normalization later).

DL Input/Output Initialization
Rand_M_bar =unique(Rand_M_bar_all(1:M_bar));
Ht_bar = reshape(Ht(Rand_M_bar,:),M_bar*K_DL,1);
...
DL_input = single(zeros(M_bar*K_DL*2,No_user_pairs));
DL_output = single(zeros(No_user_pairs,codebook_size));
DL_output_un= single(zeros(numel(Validation_Ind),codebook_size));

 Create real-valued input vectors and codebook-size output labels for each user pair.

Input Generation + Rate Labels (Training + Validation)

for pp = 1:1:numel(Ur_rows_grid)-1

 For each user batch:

o Load users.
o For each user, compute sampled channel:
Hbar=(Ht+n1)∘(Hr+n2)H_{\text{bar}} = (H_t + n_1) \circ (H_r + n_2)Hbar=(Ht
+n1)∘(Hr+n2)

DL_input(:,u+uu-1+((pp-1)*params.num_user))= reshape([real(H_bar)
imag(H_bar)].',[],1);

 Construct real-valued input: stack real + imaginary parts vertically.

R = single(log2(1+(SNR_sqrt_var/Delta_H_max).^2));
Rn = diag(1./Delta_Out_max)*R;
DL_output(...) = Rn;

 Normalize achievable rate vector per user and store as DL label.

Final Normalization
DL_input= 1*(DL_input/Delta_H_bar_max);

 Normalize entire DL input range based on worst-case magnitude.

Dataset Reshaping
DL_output_reshaped = reshape(...);
DL_input_reshaped= reshape(...);

 Reshape to match imageInputLayer expectations: [features, 1, 1, N].

Loop over Training Sizes

for dd=1:1:numel(Training_Size)

 For each dataset size, split the shuffled dataset into training and validation sets.

layers = [
imageInputLayer(...)
fullyConnectedLayer(...)
...
regressionLayer('Name','outReg')];

 Define a feedforward regression network with ReLU activations and dropout.

options = trainingOptions('sgdm', ...)

 Use Stochastic Gradient Descent with Momentum for training.

Model Training
[trainedNet,traininfo] = trainNetwork(XTrain,YTrain,layers,options);
YPredicted = predict(trainedNet,XValidation);

 Train the model and predict on validation set.

Achievable Rate Calculation
[~,Indmax_DL] = maxk(YPredicted,kbeams,2);

 For each user, find top-k predicted beams.

MaxR_DL(b) = max(squeeze(YValidation_un(1,1,Indmax_DL(b,:),b)));

 Evaluate true achievable rate of those predicted beams.

Rate_OPT(dd) = mean(MaxR_OPT);
Rate_DL(dd) = mean(MaxR_DL);

 Average achievable rates across validation users.

Final Output
end

 Return:
o Rate_DL: Achievable rates via DL-based beam prediction.
o Rate_OPT: Upper bound via ideal beam selection.

UPA CODEBOOK GENERATION

Function Definition
function [F_CB, all_beams] = UPA_codebook_generator(Mx, My, Mz,
over_sampling_x, over_sampling_y, over_sampling_z, ant_spacing)

 Inputs:
o Mx, My, Mz: Number of antenna elements in x, y, z dimensions.
o over_sampling_x/y/z: How finely to sample steering angles in each direction.
o ant_spacing: Antenna spacing in wavelengths (usually 0.5).
 Outputs:
o F_CB: Final 3D beamforming codebook (M × total_beams).
o all_beams: Indices of selected beams in each direction (total_beams × 3).

Constants and Antenna Indexing

kd = 2*pi*ant_spacing;
  kd is the wave number times antenna spacing: kd=2πdkd = 2\pi dkd=2πd.
It determines the phase shift per element due to angle-of-arrival/departure.

antx_index = 0:1:Mx-1;
anty_index = 0:1:My-1;
antz_index = 0:1:Mz-1;
M = Mx * My * Mz;

 Indices for antenna elements along each axis.

 Total number of antennas: M=Mx⋅My⋅MzM = M_x \cdot M_y \cdot M_zM=Mx⋅My⋅Mz.

Codebook Size per Dimension

codebook_size_x = over_sampling_x * Mx;
codebook_size_y = over_sampling_y * My;
codebook_size_z = over_sampling_z * Mz;

 Codebook oversamples the DFT basis in each direction for finer resolution.

X-Direction Beamforming Vectors

theta_qx = 0:pi/codebook_size_x:pi-1e-6;

 Azimuth beamsteering angles uniformly quantized from 0 to π.

F_CBx = zeros(Mx, codebook_size_x);

for i = 1:1:length(theta_qx)
F_CBx(:, i) = sqrt(1/Mx) * exp(-1j * kd * antx_index' * cos(theta_qx(i)));
end

 For each angle:

o Generate a steering vector along x.
o Normalized by √Mx for unit power.

Y-Direction Beamforming Vectors

range_y = (20 + 307) * pi / 180;
theta_qy = 20*pi/180 : -range_y/codebook_size_y : -307*pi/180 + 1e-6;

 Y-direction has a non-uniform angular range: from +20° to -307° in radians.

 theta_qy: Elevation or azimuth angles (paper-dependent choice).

F_CBy = zeros(My, codebook_size_y);

for i = 1:1:length(theta_qy)
F_CBy(:, i) = sqrt(1/My) * exp(-1j * anty_index' * theta_qy(i));
end

 Create normalized steering vectors for y-direction.

Z-Direction Beamforming Vectors

theta_qz = 0:pi/codebook_size_z:pi-1e-6;
F_CBz = zeros(Mz, codebook_size_z);
for i = 1:1:length(theta_qz)
F_CBz(:, i) = sqrt(1/Mz) * exp(-1j * kd * antz_index' * cos(theta_qz(i)));
end

 Z-direction beamsteering (e.g., elevation) from 0 to π.

 Phase shift depends on cos(θ) as in array response.

Combine 3D Codebooks
F_CBxy = kron(F_CBy, F_CBx);

 Kronecker product to create 2D UPA (x–y plane) codebook:

Mx⋅My×Nx⋅NyM_x \cdot M_y \times N_x \cdot N_yMx⋅My×Nx⋅Ny

F_CB = kron(F_CBz, F_CBxy);

 Final 3D codebook by extending to z-axis:

⋅My⋅Mz×Nx⋅Ny⋅Nz
Mx⋅My⋅Mz×Nx⋅Ny⋅NzM_x \cdot M_y \cdot M_z \times N_x \cdot N_y \cdot N_zMx

Beam Index Tracking

beams_x = 1:1:codebook_size_x;
beams_y = 1:1:codebook_size_y;
beams_z = 1:1:codebook_size_z;

 These are just index vectors for each beam angle.

Mxx_Ind = repmat(beams_x, 1, codebook_size_y * codebook_size_z)';

 Repeat beams_x to match all combinations with y and z.

o Shape: Nx⋅Ny⋅Nz×1N_x \cdot N_y \cdot N_z \times 1Nx⋅Ny⋅Nz×1

Myy_Ind = repmat(reshape(repmat(beams_y, codebook_size_x, 1), 1, ...), 1,

codebook_size_z)';
  Stretch beams_y across x and z combinations.

Mzz_Ind = reshape(repmat(beams_z, codebook_size_x * codebook_size_y,

1), ...)';

 Stretch beams_z across x and y combinations.

Pack Beam Indices

Tx = cat(3, Mxx_Ind', Myy_Ind', Mzz_Ind');
all_beams = reshape(Tx, [], 3);

 Stack all 3 beam index vectors together and reshape into:

o all_beams: [total number of beams, 3], where each row gives the beam indices
along x, y, z directions.

Summary

Output:

 F_CB: Complex-valued codebook matrix of size M × (codebook_size_x ×

codebook_size_y × codebook_size_z).
 all_beams: Index matrix of shape [total beams × 3] describing beam combinations.