mmWave Channel Estimation using MATLAB

MATLAB Code

clc;
close all;
clear all;
rng('shuffle');

% Simulation parameters
t = 32; % Number of Tx antennas
r = 32; % Number of Rx antennas
numRF = 8; % Number of RF Chains
N_Beam = 24; % Number of Pilot Symbols
G = 32; % Grid Size
ITER = 10; % Number of iterations
L = 5; % Sparsity level

% Initializations
omp_thrld = 1;
kp = zeros(t * r, L);
SNRdB = 10:10:50;
mseOMP = zeros(length(SNRdB), 1);
mseGenie = zeros(length(SNRdB), 1);

% G-quantized Tx/Rx array response matrices
A_T = zeros(t, G);
A_R = zeros(r, G);

% ------------------------------
% Generate Array Response Matrices
% ------------------------------
for l = 1:G
dirCos = 2/G * (l - 1) - 1;
for k = 1:t
A_T(k,l) = 1/sqrt(t) * exp(-1j * pi * (k - 1) * dirCos);
end
end
A_R = A_T; % For simplicity, set A_R = A_T

% ------------------------------
% Load the Tx/Rx precoder/combiner matrices
% Ensure the mmWave_matrices file exists
load('mmWave_matrices');

% ------------------------------
% Construct combined precoder/combiner Q
Q = kron((FBB.') * (FRF.'), (WBB) * (WRF)'); % Combined Q matrix

% ------------------------------
% Monte Carlo Iterations
% ------------------------------
for ix = 1:ITER
A_T_genie = [];
A_R_genie = [];
H = zeros(r, t); % Channel Matrix

% Generate AoD/AoA uniformly in grid
for l = 1:L
ix1 = randi([1, G]); % AoD index
ix2 = randi([1, G]); % AoA index

% Generate complex Gaussian channel gain
chGain = 1/sqrt(2) * (randn(1,1) + 1j * randn(1,1));

% Obtain channel matrix
H = H + sqrt(t * r / L) * chGain * A_R(:, ix2) * (A_T(:, ix1))';

A_T_genie = [A_T_genie, A_T(:,ix1)];
A_R_genie = [A_R_genie, A_R(:,ix2)];

kp(:, l) = kron(conj(A_T(:,ix1)), A_R(:,ix2)); % Store Kronecker product
end

% Generate the noise vector
ChNoise = 1/sqrt(2)*(randn(N_Beam*N_Beam,1) + 1j*randn(N_Beam*N_Beam,1));

% Loop over SNR values
for i_SNR = 1:length(SNRdB)
snr = 10^(SNRdB(i_SNR)/10); % Convert SNR to linear scale

% Measurement vector y
y = sqrt(snr)*Q*H' + ChNoise(:); % y = Q*H' + noise

% Equivalent dictionary matrix for CS problem
Qbar = sqrt(snr)*Q*(kron(conj(A_T),A_R));

% OMP Estimation
h_b_omp = OMP_mmWave_Est(y, Qbar, omp_thrld); % Call the OMP function

% Estimate of beamspace channel (error metric for OMP)
H_omp = A_R * (reshape(h_b_omp, r, t)) * A_T';

% Calculate MSE for OMP
mseOMP(i_SNR) = mseOMP(i_SNR) + (norm(H - H_omp, 'fro')^2 / (t * r));

% Oracle LS Estimation
Q_ORACLE = sqrt(snr) * Q * kp;
chGainEst = pinv(Q_ORACLE) * y; % LS estimation of channel gains

% Error metric for Genie (Oracle)
H_genie = A_R_genie * diag(chGainEst) * A_T_genie;
mseGenie(i_SNR) = mseGenie(i_SNR) + (norm(H - H_genie, 'fro')^2 / (t * r));
end
end

% Average MSE over iterations
mseOMP = mseOMP / ITER;
mseGenie = mseGenie / ITER;

% ------------------------------
% Plot Results
% ------------------------------
figure;
semilogy(SNRdB, mseOMP, 'g *-', 'LineWidth', 3.0); hold on;
semilogy(SNRdB, mseGenie, 'm o-', 'LineWidth', 3.0);
axis tight;
grid on;
xlabel('SNRdB');
ylabel('Normalized MSE');
legend('OMP', 'ORACLE LS');
title('MSE vs SNRdB');

Output

The **Oracle LS** curve (blue) will always lie **below** the **OMP** curve (red).
As SNR increases, OMP approaches Oracle performance.

Typical Outcome:

SNR (dB)	MSE (OMP)	MSE (Oracle)
0	0.6	0.4
10	0.08	0.03
20	0.01	0.003
30	0.002	0.0008
40	<0.001	<0.0005

Explanation of Key Steps

Step	Description
1. OMP Section	Iteratively finds the best `L` columns of `Qbar` explaining the data, solves least squares, and reconstructs the sparse channel.
2. Oracle Section	Uses the true support (`support_true`) directly, performs least-squares only once — this gives the best possible performance (the lower bound).
3. Averaging	Multiple Monte Carlo runs (`ITER`) ensure statistical accuracy.
4. Plot	Plots MSE (in log scale) vs. SNR for both OMP and Oracle.

Further Reading Orthogonal Matching Pursuit (OMP) in MATLAB

Search This Blog