【通信】基于 ADMM 的大规模 MIMO 无穷范数检测附matlab代码

1 内容介绍

在本文中，我们为大规模多用户 (MU) 多输入多输出 (MIMO) 无线系统提出了一种新颖的数据检测算法和相应的 VLSI 设计。我们的算法使用基于交替方向乘法器 (ADMM) 的无限范数约束均衡，称为 ADMIN。ADMIN 是一种迭代算法，当基站 (BS) 和用户天线的数量之比较小时，它的性能大大优于线性检测器。在第一次迭代中，ADMIN 计算线性最小均方误差 (MMSE) 解，这在 BS 和用户天线数量之比较大时就足够了。我们为支持 16 和 32 用户系统的基于 LDL 分解的软输出 ADMIN 开发分时和迭代 VLSI 架构。我们提出了用于 28-nm CMOS 中的 16-64 个天线基站的专用集成电路 (ASIC) 设计，支持多达 64 个正交幅度调制 (QAM)。16 用户 ADMIN ASIC 达到 303 Mb/s，同时耗散 85 mW。

2 部分代码

% -----------------------------------------------------

% -- Simulation of state-of-the-art massive MIMO detection algorithms

% -- Author: Shahriar Shahabuddin

% -- Email: shahriar.cwc@gmail.com

% This simulator is based on Christoph Studer's simple MIMO simulator. This

% simulator contains the following algorithms:

% (1) Conventional detection schemes: matched filtering, MMSE, SIMO

% (2) Approximate Inversion Based Detection: neumann-series approximation,

% gauss-seidel detection, conjugate-gradient detection

% (3) BOX detection based methods: ADMIN, OCD

% Please cite the following paper if you use this simulator,

% S. Shahabddin, M. Juntti and C. Studer, "ADMM-based infinity-norm

% detector for large-scale MIMO", IEEE International symposium of circuits

% and systems, Maryland, USA, May 2017.

% -----------------------------------------------------

function massiveMIMOdetectors(varargin)

% -- set up default/custom parameters

close all

if isempty(varargin)

disp('using default simulation settings and parameters...')

% set default simulation parameters

par.suffix = 'exp'; % simulation name suffix: 'exp' experimental

par.runId = 0; % simulation ID (used to reproduce results)

par.MR = 64; % receive antennas

par.MT = 16; % user terminals (set not larger than MR!)

par.mod = '64QAM'; % modulation type: 'BPSK','QPSK','16QAM','64QAM'

par.simName = ['ERR_' num2str(par.MR) 'x' num2str(par.MT) '_' par.mod '_' par.suffix] ; % simulation name (used for saving results)

par.trials = 100; % number of Monte-Carlo trials (transmissions)

par.SNRdB_list = 10:2:20; % list of SNR [dB] values to be simulated

par.detector = {'Conjugate-Gradient','Neumann','Gauss-Seidel','OCDBOX','ADMIN'}; % define detector(s) to be simulated

% algorithm specific

par.alg.maxiter = 3;

else

disp('use custom simulation settings and parameters...')

par = varargin{1}; % only argument is par structure

end

% -- initialization

% use runId random seed (enables reproducibility)

% rng(par.runId);

% set up Gray-mapped constellation alphabet (according to IEEE 802.11)

switch (par.mod)

case 'BPSK'

par.symbols = [ -1 1 ];

case 'QPSK'

par.symbols = [ -1-1i,-1+1i, ...

+1-1i,+1+1i ];

case '16QAM'

par.symbols = [ -3-3i,-3-1i,-3+3i,-3+1i, ...

-1-3i,-1-1i,-1+3i,-1+1i, ...

+3-3i,+3-1i,+3+3i,+3+1i, ...

+1-3i,+1-1i,+1+3i,+1+1i ];

case '64QAM'

par.symbols = [ -7-7i,-7-5i,-7-1i,-7-3i,-7+7i,-7+5i,-7+1i,-7+3i, ...

-5-7i,-5-5i,-5-1i,-5-3i,-5+7i,-5+5i,-5+1i,-5+3i, ...

-1-7i,-1-5i,-1-1i,-1-3i,-1+7i,-1+5i,-1+1i,-1+3i, ...

-3-7i,-3-5i,-3-1i,-3-3i,-3+7i,-3+5i,-3+1i,-3+3i, ...

+7-7i,+7-5i,+7-1i,+7-3i,+7+7i,+7+5i,+7+1i,+7+3i, ...

+5-7i,+5-5i,+5-1i,+5-3i,+5+7i,+5+5i,+5+1i,+5+3i, ...

+1-7i,+1-5i,+1-1i,+1-3i,+1+7i,+1+5i,+1+1i,+1+3i, ...

+3-7i,+3-5i,+3-1i,+3-3i,+3+7i,+3+5i,+3+1i,+3+3i ];

end

% extract average symbol energy

par.Es = mean(abs(par.symbols).^2);

% precompute bit labels

par.Q = log2(length(par.symbols)); % number of bits per symbol

par.bits = de2bi(0:length(par.symbols)-1,par.Q,'left-msb');

% track simulation time

time_elapsed = 0;

% -- start simulation

% initialize result arrays (detector x SNR)

res.VER = zeros(length(par.detector),length(par.SNRdB_list)); % vector error rate

res.SER = zeros(length(par.detector),length(par.SNRdB_list)); % symbol error rate

res.BER = zeros(length(par.detector),length(par.SNRdB_list)); % bit error rate

% generate random bit stream (antenna x bit x trial)

bits = randi([0 1],par.MT,par.Q,par.trials);

% trials loop

tic

for t=1:par.trials

% generate transmit symbol

idx = bi2de(bits(:,:,t),'left-msb')+1;

s = par.symbols(idx).';

% generate iid Gaussian channel matrix & noise vector

n = sqrt(0.5)*(randn(par.MR,1)+1i*randn(par.MR,1));

H = sqrt(0.5)*(randn(par.MR,par.MT)+1i*randn(par.MR,par.MT));

% transmit over noiseless channel (will be used later)

x = H*s;

% SNR loop

for k=1:length(par.SNRdB_list)

% Current SNR point in dBs

SNR_dB = par.SNRdB_list(k);

% Linear SNR

SNR_lin = 10.^(SNR_dB./10);

% Variance of complex noise per receive antenna

N0 = par.Es*par.MT/SNR_lin;

% transmit data over noisy channel

y = x+sqrt(N0)*n;

% algorithm loop

for d=1:length(par.detector)

switch (par.detector{d}) % select algorithms

case 'MF' % Matched Filter

[idxhat,bithat] = MF(par,H,y,N0);

case 'MMSE' % MMSE detector

[idxhat,bithat] = MMSE(par,H,y,N0);

case 'SIMO' % SIMO lower bound

[idxhat,bithat] = SIMO(par,H,y,N0,s);

case 'ADMIN' % ADMM-based Infinity Norm detector

[idxhat,bithat] = ADMIN(par,H,y,N0);

case 'OCDBOX' % co-ordinate descent (optimized) detector

[idxhat,bithat] = OCDBOX(par,H,y);

case 'Neumann' % coordinate descent

[idxhat,bithat] = Neumann(par,H,y,N0);

case 'Gauss-Seidel' % Gauss-Seidel detector

[idxhat,bithat] = Gauss_Seidel(par,H,y,N0);

case 'Conjugate-Gradient' % conjugate gradient detector

[idxhat,bithat] = CG(par,H,y,N0);

otherwise

error('par.detector type not defined.')

end

% -- compute error metrics

err = (idx~=idxhat);

res.VER(d,k) = res.VER(d,k) + any(err);

res.SER(d,k) = res.SER(d,k) + sum(err)/par.MT;

res.BER(d,k) = res.BER(d,k) + sum(sum(bits(:,:,t)~=bithat))/(par.MT*par.Q);

end % algorithm loop

end % SNR loop

% keep track of simulation time

if toc>10

time=toc;

time_elapsed = time_elapsed + time;

fprintf('estimated remaining simulation time: %3.0f min.\n',time_elapsed*(par.trials/t-1)/60);

tic

end

end % trials loop

% normalize results

res.VER = res.VER/par.trials;

res.SER = res.SER/par.trials;

res.BER = res.BER/par.trials;

res.time_elapsed = time_elapsed;

% -- save final results (par and res structure)

% save([ par.simName '_' num2str(par.runId) ],'par','res');

% -- show results (generates fairly nice Matlab plot)

marker_style = {'bo-','rs--','mv-.','kp:','g*-','c>--','yx:'};

figure(1)

for d=1:length(par.detector)

if d==1

semilogy(par.SNRdB_list,res.BER(d,:),marker_style{d},'LineWidth',2)

hold on

else

semilogy(par.SNRdB_list,res.BER(d,:),marker_style{d},'LineWidth',2)

end

hold off

grid on

xlabel('average SNR per receive antenna [dB]','FontSize',12)

ylabel('bit error rate (BER)','FontSize',12)

axis([min(par.SNRdB_list) max(par.SNRdB_list) 1e-4 1])

legend(par.detector,'FontSize',12)

set(gca,'FontSize',12)

end

% -- set of detector functions

%% Matched filter

function [idxhat,bithat] = MF(par,H,y)

xhat = H' * y / norm(H(:));

[~,idxhat] = min(abs(xhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% MMSE detector (MMSE)

function [idxhat,bithat] = MMSE(par,H,y,N0)

xhat = (H'*H+(N0/par.Es)*eye(par.MT))\(H'*y);

[~,idxhat] = min(abs(xhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% SIMO bound

function [idxhat,bithat] = SIMO(par,H,y,s)

z = y-H*s;

xhat = zeros(par.MT,1);

for m=1:par.MT

hm = H(:,m);

yhat = z+hm*s(m,1);

xhat(m,1) = hm'*yhat/norm(hm,2)^2;

end

[~,idxhat] = min(abs(xhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% Neumann-Series Approximation based massive MIMO detection

function [idxhat,bithat] = Neumann(par,H,y,N0)

A = H'*H+(N0/par.Es)*eye(par.MT);

MF = H'*y;

D = diag(diag(A));

E = triu(A,1)+tril(A,-1);

Ainv = 0;

for i = 0:par.alg.maxiter

Ainv = Ainv+((-inv(D)*E)^i)*inv(D);

end

xhat = Ainv*MF;

[~,idxhat] = min(abs(xhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% Gauss-Seidel massive MIMO detection

function [idxhat,bithat] = Gauss_Seidel(par,H,y,N0)

A = H'*H+(N0/par.Es)*eye(par.MT);

MF = H'*y;

D = diag(diag(A));

E = -triu(A,1);

F = -tril(A,-1);

xhat = diag(inv(D));% inv(D)*MF; %%% Check Gauss Seidel detection paper

for i = 0:par.alg.maxiter

xhat = inv(D-E)*(F*xhat+MF);

end

[~,idxhat] = min(abs(xhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% Conjugate Gradient massive MIMO detection

function [idxhat,bithat] = CG(par,H,y,N0)

A = H'*H+(N0/par.Es)*eye(par.MT);

MF = H'*y;

r = MF;

p = r;

v = zeros(par.MT,1);

for k = 1:par.alg.maxiter

e = A*p;

alpha = norm(r)^2/(p'*e);

v = v+alpha*p;

new_r = r-alpha*e;

beta = norm(new_r)^2/norm(r)^2;

p = new_r+beta*p;

r = new_r;

end

xhat = v;

[~,idxhat] = min(abs(xhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% ADMM-based infinity norm (ADMIN) detector

function [idxhat,bithat] = ADMIN(par,H,y,N0)

% -- preprocessing

% by setting beta to N0/par.Es we get the MMSE estimator in the first iteration

% this is pretty neat as this is a very good detector already

beta = N0/par.Es;%*3; % tweaking this one by 3 improved performance significantly

A = H'*H + beta*eye(par.MT);

L = chol(A,'lower');

yMF = H'*y;

% -- initialization

gamma = (1+sqrt(5))/2;%*2; %% tweaked with 2 to improve performance

alpha = max(real(par.symbols)); % symbol box

zhat = zeros(par.MT,1);

lambda = zeros(par.MT,1);

% -- ADMM loop

for iter=1:par.alg.maxiter

xhat = (L')\(L\(yMF+beta*(zhat-lambda))); % step 1

zhat = projinf(par,xhat+lambda,alpha); % step 2

lambda = lambda-real(gamma*(zhat-xhat)); % step 3

lambda = real(lambda);

end

% -- hard output detection

[~,idxhat] = min(abs(zhat*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

%% Optimized Coordinate Descent (OCD) BOX version

function [idxhat,bithat] = OCDBOX(par,H,y)

% -- initialization

[row, col] = size(H);

alpha = 0; % no regularization for BOX detector

beta = max(real(par.symbols));

% -- preprocessing

dinv = zeros(col,1);

p = zeros(col,1);

for uu=1:col

normH2 = norm(H(:,uu),2)^2;

dinv(uu,1) = 1/(normH2+alpha);

p(uu,1) = dinv(uu)*normH2;

end

r = y;

zold = zeros(col,1);

znew = zeros(col,1);

deltaz = zeros(col,1);

% -- OCD loop

for iters=1:par.alg.maxiter

for uu=1:col

tmp = dinv(uu)*(H(:,uu)'*r)+p(uu)*zold(uu);

znew(uu) = projinf(par,tmp,beta);

deltaz(uu) = znew(uu)-zold(uu);

r = r - H(:,uu)*deltaz(uu);

zold(uu) = znew(uu);

end

[~,idxhat] = min(abs(znew*ones(1,length(par.symbols))-ones(par.MT,1)*par.symbols).^2,[],2);

bithat = par.bits(idxhat,:);

end

% project onto alpha infinity-tilde-norm ball

function sproj = projinf(par,s,alpha)

switch par.mod

case 'BPSK'

v = real(s);

sproj = min(abs(v),alpha).*v;

otherwise

sr = real(s);

idxr = abs(sr)>alpha;

sr(idxr) = sign(sr(idxr))*alpha;

si = imag(s);

idxi = abs(si)>alpha;

si(idxi) = sign(si(idxi))*alpha;

sproj = sr +1i*si;

end

3 运行结果

4 参考文献

[1] Shahabuddin S , Hautala I , Juntti M , et al. ADMM-Based Infinity-Norm Detection for Massive MIMO: Algorithm and VLSI Architecture[J]. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 2021, PP(99):1-13.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机、雷达通信、无线传感器等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。