初识MIMO(六):MU-MIMO的仿真

初识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
image-20220504223857351

image-20220504223916370

二. 信道传播方式一——信道反转

1.理论

image-20220504224649784

image-20220504224706302

信道反转就是所谓的预编码均衡,没有任何区别,但注意,此时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.理论

image-20220504230645054

image-20220504230827535

image-20220504230841616

image-20220504230853055

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.一些思考

  1. 代码我做了推广,可以适用于任意的NT和NR,但都要满足NT=NR*N_user 是严格按照书本公式来的,不过随着NT增大,效果越差这,这也很好理解,因为干扰越来越难以消除
  2. 我们可以看到,块对角化最大的特点是(13.17),他考虑了不属于当前MS的信号对当前MS的影响,并尝试消除他,就是说专门的NT服务专门的MS。

四. 信道传播方式三——DPC和TH

image-20220504233017488

image-20220504233041869

image-20220504233052099

image-20220504233106267

image-20220504233137585

image-20220504233147014

image-20220504233156184

image-20220504233205786

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.一些思考

  1. % 选出四个最好的信道
    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,:)’😅;

这里在做什么呢?

这里是在选出最好的四个信号

排序后最小的最大的就是四个最大的啦

  1. % Matlab没有LQ分解,要这样曲线救国
    [Q_temp,R_temp] = qr(H_selected’);
    L=R_temp’;
    Q=Q_temp’;

image-20220504234528003

  1. frame_prepre(ipre,:) = frame_prepre(ipre,:) - L(ipre,1:ipre-1)/L(ipre,ipre)*frame_prepre(1:ipre-1,:);

这是

image-20220504234753058

  1. inv(diag(diag(L)))

这里是matlab的特性,diag(L)得到的是一个数组,diag(diag(L))则提取出对角元素

  1. 为什么叫脏纸编码?

image-20220504235039976

这是非常好的比喻,因为编码正是挨个解决问题的

image-20220504235121757

总结

其实MU-MIMO最终都是希望得到一个易于分解的对角矩阵,为此有各种方法,其实每种方法都对应着一种分解,对角化对应着SVD,DPC则对应着LQ分解,如果还有分解方法,一定还有更多的分解方法

MIMO-OFDM的matlab仿真到此就结束了,后续如果有需要会更新之前没有看懂的内容,非常感谢这本书以及陈老湿的博客带我入门了通信!我后续也会继续学习,更新!

  • 6
    点赞
  • 36
    收藏
    觉得还不错? 一键收藏
  • 8
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 8
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值