初识MIMO(六):MU-MIMO的仿真
零 代码地址
https://github.com/liu-zongxi/MIMO_simulation
请大家看完觉得有用别忘了点赞收藏,github项目给star哦
一.基本概念
多用户MIMO主要是希望实现一对多,这要求对于每一个MS,要去除掉发送给其他MS的信号,这也就是这一章的核心了,在接下来所有的仿真中,我们都更关注下行广播,而我感觉到所有仿真里都隐含了一个条件就是
N
T
=
N
R
×
N
u
s
e
r
NT = NR \times N_{user}
NT=NR×Nuser
二. 信道传播方式一——信道反转
1.理论
信道反转就是所谓的预编码均衡,没有任何区别,但注意,此时NR=1,直接去掉信道的影响来看代码吧
2. 代码
%------------------MU-MIMO的信道反转方式---------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月3日15点50分-----------------%
%% 参数设置
clear;clf;
L_frame = 100; % 一个帧里多少符号
N_packet = 2000; % 一个包里多少帧
Nmod = 2; % 调制阶数
NT = 4; % 发射天线数
N_user = 20; % 用户数
N_user_active = NT; % 当前使用的用户数
I = eye(NT,NT);
Ncase = 4;
SNRs_dB = 0:2:20;
SNRs = 10.^(SNRs_dB/10);
N_SNR = length(SNRs);
BERs = zeros(Ncase, N_SNR);
gss = ["-kx" "-b^" "-ro" "-m+"]; % 画图图像,注意使用双引号
%% 主函数
for icase = 1:Ncase
gs = gss(icase);
if icase ==1
W_formula = @(H, sigma, I) H'*inv(H*H');
flag_select = 0;
elseif icase ==2
W_formula = @(H, sigma, I) H'*inv(H*H'+sigma*I);
flag_select = 0;
elseif icase ==3
W_formula = @(H, sigma, I) H'*inv(H*H');
flag_select = 1;
elseif icase ==4
W_formula = @(H, sigma, I) H'*inv(H*H'+sigma*I);
flag_select = 1;
end
for isnr = 1:N_SNR
SNR = SNRs(isnr);
n_biterror = 0;
sigma = sqrt(NT/(2*SNR));
for ipacket = 1:N_packet
% 生成信号
frame_origin = randi([0,1],L_frame,N_user_active*Nmod);
frame_mod = QPSKMod(frame_origin,L_frame, N_user_active);
frame_reshape = reshape(frame_mod, N_user_active, L_frame);
% 生成信道
if flag_select == 0
H_DL = (randn(N_user_active,NT)+1j*randn(N_user_active,NT))/sqrt(2);
elseif flag_select ==1
H = zeros(N_user, NT);
Channel_norm = zeros(1, N_user);
for iuser = 1:N_user
H(iuser,:) = (randn(1,NT)+1j*randn(1,NT))/sqrt(2);
Channel_norm(iuser)=norm(H(iuser,:));
end
[Ch_norm,Index]=sort(Channel_norm,'descend');
H_DL = H(Index(1:N_user_active),:);
end
W = W_formula(H_DL, sigma, I);
beta = sqrt(NT/trace(W*W'));
W_pre = beta*W;
frame_trnasmit = W_pre*frame_reshape;
noise = sigma*(randn(N_user_active,L_frame)+1j*randn(N_user_active,L_frame));
y = H_DL * frame_trnasmit + noise;
% 接收
frame_re = y/beta;
frame_pre_demod = reshape(frame_re, L_frame, N_user_active);
frame_demod = QPSKDemod(frame_pre_demod,L_frame,N_user_active);
% 计算误码率
n_biterror_tmp = sum(sum(abs(frame_demod - frame_origin)))
n_biterror = n_biterror + n_biterror_tmp;
end
BERs(icase, isnr) = n_biterror/(N_packet*L_frame*Nmod);
end
semilogy(SNRs_dB,BERs(icase,:),gs);
hold on;
end
%% 画图
title('MU-MIMO BER inv');
xlabel('SNR[dB]');
ylabel('BER')
grid on;
set(gca,'fontsize',9)
legend('ZF','MMSE','ZF with select', 'MMSE with select')
3. 一些思考
啊,没有什么好思考的,上一张都思考过了,看不懂这个代码的可以看我的前一个博客~
三. 信道传播方式二——块对角化
1.理论
2.代码
%------------------MU-MIMO的对角块化---------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月3日15点50分-----------------%
%% 设置
clear;
clf;
L_frame = 100;
N_packet = 2000;
Nmod = 2;
NT = 4;
NR = 2;
N_user = 2;
gss = ["-kx" "-^" "-ro" "-b>" "-g<" "-m+"];
SNRs_dB = 0:2:20;
SNRs = 10.^(SNRs_dB/10);
N_SNR = length(SNRs_dB);
BERs = zeros(5, N_SNR);
%% 主函数
for isnr = 1:N_SNR
SNR = SNRs(isnr);
n_biterror = 0;
sigma = sqrt(NT/(2*SNR));
for ipacket = 1:N_packet
% 生成信号
frame_origin = randi([0,1],L_frame,NT*Nmod);
frame_mod = QPSKMod(frame_origin,L_frame, NT);
frame_reshape = reshape(frame_mod, NT, L_frame);
% 信道
noise = sigma* (randn(NR,L_frame)+1j*randn(NR,L_frame));
H1 = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);
H2 = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);
H = (randn(NR,NT, N_user)+1j*randn(NR,NT, N_user))/sqrt(2);
% 对角化处理
for iuser = 1:N_user
if iuser == 1
region_tilde = [2:N_user];
elseif iuser == N_user
region_tilde = [1:N_user - 1];
else
region_tilde = [1:iuser-1 iuser+1:N_user];
end
H_tilde_temp = [];
for ii = region_tilde
H_tilde_temp = [H_tilde_temp H(:,:,ii)'];
end
H_tilde = H_tilde_temp';
[U,S,V] = svd(H_tilde);
index_zero = size(S,2)-rank(S);
W(:,:,iuser) = V(:,end-index_zero+1:end);
end
frame_transmit = 0;
for iuser = 1:N_user
frame_transmit = frame_transmit+W(:,:,iuser)*frame_reshape(NR*(iuser-1)+1:NR*(iuser-1)+NR,:);
end
% [U1,S1,V1] = svd(H1);
% W2 = V1(:,3:4);
% [U2,S2,V2] = svd(H2);
% W1 = V2(:,3:4);
% frame_transmit = W1*frame_reshape(1:2,:) + W2*frame_reshape(3:4,:);
for iuser = 1:N_user
y(:,:,iuser) = H(:,:,iuser)*frame_transmit+noise;
end
% y1 = H1*frame_transmit+noise;
% y2 = H2*frame_transmit+noise;
% HV1 = H1*W1;
% EQ1 = HV1'*inv(HV1*HV1'); % Equalizer for the 1st user
% HV2 = H2*W2;
% EQ2 = HV2'*inv(HV2*HV2'); % Equalizer for the 2nd user
y_re=[];
for iuser = 1:N_user
HV = H(:,:,iuser)*W(:,:,iuser);
EQ = HV'*inv(HV*HV');
y_re_temp = EQ*y(:,:,iuser);
y_re = [y_re;y_re_temp];
end
% y = [EQ1*y1; EQ2*y2];
% frame_pre_demod = reshape(y,L_frame,NT);
frame_pre_demod = reshape(y_re,L_frame,NT);
frame_demod = QPSKDemod(frame_pre_demod,L_frame,NT);
% 计算误码率
n_biterror_tmp = sum(sum(abs(frame_demod - frame_origin)))
n_biterror = n_biterror + n_biterror_tmp;
end
BERs(1, isnr) = n_biterror/(N_packet*L_frame*Nmod);
end
%% 画图
BERs(2,:) = [1.40608250000000 1.27856000000000 1.12480750000000 0.941037500000000 0.738272500000000 0.582960000000000 0.417507500000000 0.287120000000000 0.197902500000000 0.115562500000000 0.0775400000000000];
BERs(3,:) = [0.892967500000000 0.736142500000000 0.576275000000000 0.428812500000000 0.306602500000000 0.212430000000000 0.143212500000000 0.0964100000000000 0.0644250000000000 0.0476550000000000 0.0363025000000000];
BERs(4,:) = [1.16728000000000 1.02494250000000 0.830452500000000 0.641750000000000 0.481362500000000 0.320395000000000 0.211695000000000 0.144177500000000 0.0884050000000000 0.0634100000000000 0.0345100000000000];
BERs(5,:) = [0.679290000000000 0.521507500000000 0.376435000000000 0.253155000000000 0.160840000000000 0.102062500000000 0.0656150000000000 0.0394150000000000 0.0272925000000000 0.0190700000000000 0.0148800000000000];
for i = 1:5
semilogy(SNRs_dB,BERs(i,:), gss(i));
hold on;
end
3.一些思考
- 代码我做了推广,可以适用于任意的NT和NR,但都要满足NT=NR*N_user 是严格按照书本公式来的,不过随着NT增大,效果越差这,这也很好理解,因为干扰越来越难以消除
- 我们可以看到,块对角化最大的特点是(13.17),他考虑了不属于当前MS的信号对当前MS的影响,并尝试消除他,就是说专门的NT服务专门的MS。
四. 信道传播方式三——DPC和TH
2.代码
%------------------------DPC和TH编码-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月4日10点38分-----------------%
%% 参数设置
clear;clf;
L_frame =100;
N_packet = 2000;
Nmod = 2;
NT = 4;
N_user = 10;
N_user_ative = 4;
I=eye(NT,NT);
SNRs_dB = 0:2:20;
SNRs = 10.^(SNRs_dB/10);
N_SNR = length(SNRs);
BERs = zeros(2, N_SNR);
Ncase = 2;
gss = ["-kx" "-b^" "-ro" "-m+"]; % 画图图像,注意使用双引号
%% 主函数
for icase = 1:Ncase
gs = gss(icase);
for isnr = 1:N_SNR
SNR = SNRs(isnr);
n_biterror = 0;
sigma = sqrt(NT/(2*SNR));
for ipacket = 1:N_packet
% 生成数据
frame_origin = randi([0,1],L_frame,NT*Nmod);
% QPSK调制
frame_mod=QPSKMod(frame_origin,L_frame, NT);
frame_reshape = reshape(frame_mod,NT,L_frame);
% 信道
H = (randn(N_user,NT)+1j*randn(N_user,NT))/sqrt(2);
% 选出四个最好的信道
Types = nchoosek([1:N_user],N_user_ative);
for itype = 1:size(Types, 1)
H_sel = H(Types(itype,:)',:);
[Q_H, R_H] = qr(H_sel);
L_min(itype) = min(diag(R_H));
end
% 最大的最小值就是最好的四个
[val_max_L_min,index_max_L_min] = max(L_min);
H_selected = H(Types(index_max_L_min,:)',:);
% Matlab没有LQ分解,要这样曲线救国
[Q_temp,R_temp] = qr(H_selected');
L=R_temp';
Q=Q_temp';
% 预编码
frame_prepre = frame_reshape;
if icase == 1
for ipre = 2:4
frame_prepre(ipre,:) = frame_prepre(ipre,:) - L(ipre,1:ipre-1)/L(ipre,ipre)*frame_prepre(1:ipre-1,:);
end
else
for ipre = 2:4
frame_prepre(ipre,:) = ModTH(frame_prepre(ipre,:) - L(ipre,1:ipre-1)/L(ipre,ipre)*frame_prepre(1:ipre-1,:), 2);
end
end
frame_pre = Q'*frame_prepre;
% 信道
noise = sigma*(randn(N_user_ative,L_frame)+1j*randn(N_user_ative,L_frame));
y = H_selected*frame_pre+noise;
% 接收
frame_re = inv(diag(diag(L)))*y;
frame_recieved = reshape(frame_re,NT*L_frame,1);
if icase==2 % in the case of TH precoding
frame_recieved = ModTH(frame_recieved,2);
end
frame_pre_demod = reshape(frame_recieved,L_frame,NT);
frame_demod = QPSKDemod(frame_pre_demod,L_frame,NT);
% 计算误码率
n_biterror_tmp = sum(sum(abs(frame_demod - frame_origin)))
n_biterror = n_biterror + n_biterror_tmp;
end
BERs(icase, isnr) = n_biterror/(N_packet*L_frame*Nmod);
end
semilogy(SNRs_dB,BERs(icase,:),gs);
hold on;
end
3.一些思考
- % 选出四个最好的信道
Types = nchoosek([1:N_user],N_user_ative);
for itype = 1:size(Types, 1)
H_sel = H(Types(itype,:)‘😅;
[Q_H, R_H] = qr(H_sel);
L_min(itype) = min(diag(R_H));
end
% 最大的最小值就是最好的四个
[val_max_L_min,index_max_L_min] = max(L_min);
H_selected = H(Types(index_max_L_min,:)’😅;
这里在做什么呢?
这里是在选出最好的四个信号
排序后最小的最大的就是四个最大的啦
- % Matlab没有LQ分解,要这样曲线救国
[Q_temp,R_temp] = qr(H_selected’);
L=R_temp’;
Q=Q_temp’;
- frame_prepre(ipre,:) = frame_prepre(ipre,:) - L(ipre,1:ipre-1)/L(ipre,ipre)*frame_prepre(1:ipre-1,:);
这是
- inv(diag(diag(L)))
这里是matlab的特性,diag(L)得到的是一个数组,diag(diag(L))则提取出对角元素
- 为什么叫脏纸编码?
这是非常好的比喻,因为编码正是挨个解决问题的
总结
其实MU-MIMO最终都是希望得到一个易于分解的对角矩阵,为此有各种方法,其实每种方法都对应着一种分解,对角化对应着SVD,DPC则对应着LQ分解,如果还有分解方法,一定还有更多的分解方法
MIMO-OFDM的matlab仿真到此就结束了,后续如果有需要会更新之前没有看懂的内容,非常感谢这本书以及陈老湿的博客带我入门了通信!我后续也会继续学习,更新!