MATLAB Code
close all; clear; rng('shuffle'); clc;
% Simulation parameters
t = 32; % Number of Tx/Rx Antennas
r = 32; % Number of Tx/Rx Antennas
numRF = 6; % Number of RF Chains
G = 64; % Grid Size
L = 8; % Grid Size
Ns = 6; % Sparsity level
ITER = 100; % Number of iterations
% Initializations
H = zeros(r,t); % Channel matrix
SNRdB = -5:1.5:15; % SNR in dB
C_HYB = zeros(length(SNRdB), 1); % Capacity of Hybrid MIMO
C_MIMO = zeros(length(SNRdB), 1); % Capacity of Conventional MIMO
A_T = zeros(t, G); % Transmit response matrix
A_R = zeros(r, G); % Receive response matrix
% Generate transmit array response matrix
for l = 1:G
dirCos = 2/G*(l-1) - 1; % Direction cosine for each grid point
for K = 1:t
A_T(K, l) = 1/sqrt(t) * exp(-j*pi*(K-1)*dirCos);
end
end
% Assume A_R = A_T for simplicity
A_R = A_T;
% Main simulation loop
for iter1 = 1:ITER
% ===== Channel Generation =====
A_T_genie = zeros(t, L);
A_R_genie = zeros(r, L);
% Generate AoD/AoA uniformly in grid
AoDlist = randperm(G, L);
AoAlist = randperm(G, L);
% Tx/Rx array response matrices
A_T_genie = A_T(:, AoDlist);
A_R_genie = A_R(:, AoAlist);
% Channel gain
chGain = 1/sqrt(2)*(randn(L,1) + j*randn(L,1));
% Channel matrix
H = sqrt(t*r/L) * A_R_genie * diag(chGain) * A_T_genie';
% ===== SVD of Channel =====
[U,S,V] = svd(H);
Fopt = V(:,1:Ns); % Optimal unconstrained precoder
% ===== OMP-Based Transmit Hybrid Precoder Design =====
[FBB, FRF] = SOMP_mmW_precoder(Fopt, A_T, eye(t), numRF);
FBB_NORM = sqrt(Ns)/(norm(FRF*FBB,'fro')) * FBB; % Normalized baseband precoder
% ===== Capacity Evaluation over SNR =====
for i_snr = 1:length(SNRdB)
np = 10^(-SNRdB(i_snr)/10); % Noise power for signal power = 1
% --- Optimal unconstrained MMSE combiner for full MIMO ---
Wmmse_opt = H*Fopt*inv(Fopt'*H'*H*Fopt + np*Ns*eye(Ns));
% Capacity of unconstrained MIMO system
disp('Size of Wmmse_opt''*H*Fopt:');
disp(size(Wmmse_opt'*H*Fopt)); % Should output [Ns, Ns]
disp('Size of 1/Ns*eye(Ns):');
disp(size(1/Ns*eye(Ns))); % Should output [Ns, Ns]
% Call to mimo_capacity
C_MIMO(i_snr) = C_MIMO(i_snr) + ...
mimo_capacity(Wmmse_opt'*H*Fopt, 1/Ns*eye(Ns), np); % Corrected: Use Wmmse_opt
% --- Optimal unconstrained MMSE combiner for hybrid precoder ---
Ryy = 1/Ns * H*FRF*FBB_NORM*FBB_NORM'*FRF'*H' + np*eye(r);
Wmmse_Hyb = H*FRF*FBB_NORM*inv(FBB_NORM'*FRF'*H'*H*FRF*FBB_NORM + np*Ns*eye(Ns));
% --- OMP-based receive hybrid combiner design ---
[WBB, WRF] = SOMP_mmW_precoder(Wmmse_Hyb, A_R, Ryy, numRF); % Ensure WBB and WRF are computed here
% --- Capacity of hybrid precoder/combiner ---
% np is the noise power, so use it directly as the noise power scalar
noisePower = np;
C_HYB(i_snr) = C_HYB(i_snr) + ...
mimo_capacity(WBB'*WRF'*H*FRF*FBB_NORM, 1/Ns*eye(Ns), noisePower);
end
end
C_MIMO = C_MIMO / ITER;
C_HYB = C_HYB / ITER;
plot(SNRdB, C_MIMO, 'b', 'linewidth', 3.0);
hold on;
plot(SNRdB, C_HYB, 'm-.', 'linewidth', 3.0);
grid on;
axis tight;
xlabel('SNR (dB)');
ylabel('Capacity (b/s/Hz)');
legend('Conventional MIMO', 'Hybrid MIMO');
title('Capacity vs SNR');
% Function: SOMP_mmW_precoder
function [FBB, FRF] = SOMP_mmW_precoder(Fopt, A_T, Ryy, numRF)
% SOMP_mmW_precoder implements Sparse Orthogonal Matching Pursuit (SOMP)
% Inputs:
% Fopt - Optimal unconstrained precoder (full dimension)
% A_T - Transmit array response matrix
% Ryy - Covariance matrix (used in receive hybrid precoder case)
% numRF - Number of RF chains (desired rank of hybrid precoder)
% Outputs:
% FBB - Baseband precoder (lower-dimensional)
% FRF - RF precoder (upper-dimensional)
[t, G] = size(A_T); % t: Number of transmit antennas, G: Grid size
L = numRF; % Number of RF chains (rank of the hybrid precoder)
FBB = zeros(L, G); % Baseband precoder (output)
FRF = zeros(t, L); % RF precoder (output)
% Step 1: Initialize the residual as the full precoder (Fopt)
res = Fopt;
selectedColumns = [];
% Step 2: Iterative column selection
for l = 1:L
% Compute the correlation of each column of A_T with the residual
corr = abs(A_T' * res); % Correlation between the residual and array response
[~, idx] = max(corr); % Find index of the best matching column
selectedColumns = [selectedColumns, idx]; % Store the selected index
% Step 3: Update the RF precoder matrix with the selected column
FRF(:, l) = A_T(:, selectedColumns(end)); % Assign the selected column to RF precoder
% Update the residual by removing the contribution of the selected column
res = res - A_T(:, selectedColumns(end)) * (A_T(:, selectedColumns(end))' * res);
end
% Step 4: Solve for the baseband precoder using least squares
% Now, FBB should be computed such that FRF * FBB approximates Fopt.
% Since we have L columns in FRF and L rows in FBB, we use FRF' to solve for FBB.
FBB = pinv(FRF' * FRF) * FRF' * Fopt; % Pseudo-inverse solution for least squares
end
% Function: mimo_capacity
function capacity = mimo_capacity(H, W, noisePower)
% mimo_capacity calculates the capacity of a MIMO system.
% Inputs:
% H - Channel matrix (r x t)
% W - Receiver combiner matrix (r x Ns)
% noisePower - Noise power (scalar, typically 1 or computed from SNR)
% Outputs:
% capacity - MIMO capacity in b/s/Hz
% Effective received signal matrix after applying combiner W
effectiveH = H * W; % (r x Ns), where r is number of receive antennas, Ns is number of streams
% Compute the SNR matrix
snrMatrix = effectiveH' * effectiveH; % (Ns x Ns)
% Check the dimensions of snrMatrix and ensure it is square
[Ns, Ns_check] = size(snrMatrix);
if Ns ~= Ns_check
error('SNR matrix must be square: Ns x Ns');
end
% Compute the capacity using the Shannon capacity formula:
% C = log2(det(I + (1/noisePower) * snrMatrix))
% Identity matrix of the same size as snrMatrix
I = eye(Ns);
% Check if noisePower is scalar and valid
if numel(noisePower) ~= 1
error('noisePower must be a scalar');
end
% Compute the capacity
capacity = log2(det(I + (1 / noisePower) * snrMatrix));
end
% Simulation parameters
t = 32; % Number of Tx/Rx Antennas
r = 32; % Number of Tx/Rx Antennas
numRF = 6; % Number of RF Chains
G = 64; % Grid Size
L = 8; % Grid Size
Ns = 6; % Sparsity level
ITER = 100; % Number of iterations
% Initializations
H = zeros(r,t); % Channel matrix
SNRdB = -5:1.5:15; % SNR in dB
C_HYB = zeros(length(SNRdB), 1); % Capacity of Hybrid MIMO
C_MIMO = zeros(length(SNRdB), 1); % Capacity of Conventional MIMO
A_T = zeros(t, G); % Transmit response matrix
A_R = zeros(r, G); % Receive response matrix
% Generate transmit array response matrix
for l = 1:G
dirCos = 2/G*(l-1) - 1; % Direction cosine for each grid point
for K = 1:t
A_T(K, l) = 1/sqrt(t) * exp(-j*pi*(K-1)*dirCos);
end
end
% Assume A_R = A_T for simplicity
A_R = A_T;
% Main simulation loop
for iter1 = 1:ITER
% ===== Channel Generation =====
A_T_genie = zeros(t, L);
A_R_genie = zeros(r, L);
% Generate AoD/AoA uniformly in grid
AoDlist = randperm(G, L);
AoAlist = randperm(G, L);
% Tx/Rx array response matrices
A_T_genie = A_T(:, AoDlist);
A_R_genie = A_R(:, AoAlist);
% Channel gain
chGain = 1/sqrt(2)*(randn(L,1) + j*randn(L,1));
% Channel matrix
H = sqrt(t*r/L) * A_R_genie * diag(chGain) * A_T_genie';
% ===== SVD of Channel =====
[U,S,V] = svd(H);
Fopt = V(:,1:Ns); % Optimal unconstrained precoder
% ===== OMP-Based Transmit Hybrid Precoder Design =====
[FBB, FRF] = SOMP_mmW_precoder(Fopt, A_T, eye(t), numRF);
FBB_NORM = sqrt(Ns)/(norm(FRF*FBB,'fro')) * FBB; % Normalized baseband precoder
% ===== Capacity Evaluation over SNR =====
for i_snr = 1:length(SNRdB)
np = 10^(-SNRdB(i_snr)/10); % Noise power for signal power = 1
% --- Optimal unconstrained MMSE combiner for full MIMO ---
Wmmse_opt = H*Fopt*inv(Fopt'*H'*H*Fopt + np*Ns*eye(Ns));
% Capacity of unconstrained MIMO system
disp('Size of Wmmse_opt''*H*Fopt:');
disp(size(Wmmse_opt'*H*Fopt)); % Should output [Ns, Ns]
disp('Size of 1/Ns*eye(Ns):');
disp(size(1/Ns*eye(Ns))); % Should output [Ns, Ns]
% Call to mimo_capacity
C_MIMO(i_snr) = C_MIMO(i_snr) + ...
mimo_capacity(Wmmse_opt'*H*Fopt, 1/Ns*eye(Ns), np); % Corrected: Use Wmmse_opt
% --- Optimal unconstrained MMSE combiner for hybrid precoder ---
Ryy = 1/Ns * H*FRF*FBB_NORM*FBB_NORM'*FRF'*H' + np*eye(r);
Wmmse_Hyb = H*FRF*FBB_NORM*inv(FBB_NORM'*FRF'*H'*H*FRF*FBB_NORM + np*Ns*eye(Ns));
% --- OMP-based receive hybrid combiner design ---
[WBB, WRF] = SOMP_mmW_precoder(Wmmse_Hyb, A_R, Ryy, numRF); % Ensure WBB and WRF are computed here
% --- Capacity of hybrid precoder/combiner ---
% np is the noise power, so use it directly as the noise power scalar
noisePower = np;
C_HYB(i_snr) = C_HYB(i_snr) + ...
mimo_capacity(WBB'*WRF'*H*FRF*FBB_NORM, 1/Ns*eye(Ns), noisePower);
end
end
C_MIMO = C_MIMO / ITER;
C_HYB = C_HYB / ITER;
plot(SNRdB, C_MIMO, 'b', 'linewidth', 3.0);
hold on;
plot(SNRdB, C_HYB, 'm-.', 'linewidth', 3.0);
grid on;
axis tight;
xlabel('SNR (dB)');
ylabel('Capacity (b/s/Hz)');
legend('Conventional MIMO', 'Hybrid MIMO');
title('Capacity vs SNR');
% Function: SOMP_mmW_precoder
function [FBB, FRF] = SOMP_mmW_precoder(Fopt, A_T, Ryy, numRF)
% SOMP_mmW_precoder implements Sparse Orthogonal Matching Pursuit (SOMP)
% Inputs:
% Fopt - Optimal unconstrained precoder (full dimension)
% A_T - Transmit array response matrix
% Ryy - Covariance matrix (used in receive hybrid precoder case)
% numRF - Number of RF chains (desired rank of hybrid precoder)
% Outputs:
% FBB - Baseband precoder (lower-dimensional)
% FRF - RF precoder (upper-dimensional)
[t, G] = size(A_T); % t: Number of transmit antennas, G: Grid size
L = numRF; % Number of RF chains (rank of the hybrid precoder)
FBB = zeros(L, G); % Baseband precoder (output)
FRF = zeros(t, L); % RF precoder (output)
% Step 1: Initialize the residual as the full precoder (Fopt)
res = Fopt;
selectedColumns = [];
% Step 2: Iterative column selection
for l = 1:L
% Compute the correlation of each column of A_T with the residual
corr = abs(A_T' * res); % Correlation between the residual and array response
[~, idx] = max(corr); % Find index of the best matching column
selectedColumns = [selectedColumns, idx]; % Store the selected index
% Step 3: Update the RF precoder matrix with the selected column
FRF(:, l) = A_T(:, selectedColumns(end)); % Assign the selected column to RF precoder
% Update the residual by removing the contribution of the selected column
res = res - A_T(:, selectedColumns(end)) * (A_T(:, selectedColumns(end))' * res);
end
% Step 4: Solve for the baseband precoder using least squares
% Now, FBB should be computed such that FRF * FBB approximates Fopt.
% Since we have L columns in FRF and L rows in FBB, we use FRF' to solve for FBB.
FBB = pinv(FRF' * FRF) * FRF' * Fopt; % Pseudo-inverse solution for least squares
end
% Function: mimo_capacity
function capacity = mimo_capacity(H, W, noisePower)
% mimo_capacity calculates the capacity of a MIMO system.
% Inputs:
% H - Channel matrix (r x t)
% W - Receiver combiner matrix (r x Ns)
% noisePower - Noise power (scalar, typically 1 or computed from SNR)
% Outputs:
% capacity - MIMO capacity in b/s/Hz
% Effective received signal matrix after applying combiner W
effectiveH = H * W; % (r x Ns), where r is number of receive antennas, Ns is number of streams
% Compute the SNR matrix
snrMatrix = effectiveH' * effectiveH; % (Ns x Ns)
% Check the dimensions of snrMatrix and ensure it is square
[Ns, Ns_check] = size(snrMatrix);
if Ns ~= Ns_check
error('SNR matrix must be square: Ns x Ns');
end
% Compute the capacity using the Shannon capacity formula:
% C = log2(det(I + (1/noisePower) * snrMatrix))
% Identity matrix of the same size as snrMatrix
I = eye(Ns);
% Check if noisePower is scalar and valid
if numel(noisePower) ~= 1
error('noisePower must be a scalar');
end
% Compute the capacity
capacity = log2(det(I + (1 / noisePower) * snrMatrix));
end